Using publicly available GPS solutions for fast estimations of first-order source details from coseismic deformations

We here explore the potential use of publicly available GPS solutions to obtain first-order constraints on a source model immediately following an earthquake, within the limits of GPS solution timeliness and near-field coverage. We use GPS solutions from the Scripps Orbit and Permanent Array Center to carry out simple inversions of the coseismic displacement field induced by the 2010 Maule earthquake (Chile), by inferring the seismic moment and the rake angle of a fixed-geometry seismic source. The rake angle obtained from the inversion (m=117.8 ̊) is consistent with seismological estimates. The seismic moment, which corresponds to a moment magnitude MW =8.9, is about 1.6 times greater than seismological estimates. This suggests that as in other recent megathrust events, a consistent fraction of the energy was released aseismically. In this respect, the additional information obtained from GPS can help to provide a better estimate of the weight of the aseismic contribution to the energy release.


Introduction
GPS solutions have very important roles in highprecision geodetic and geophysical measurements, particularly for the study of earthquake-induced permanent deformations, tectonic plate motions, plate boundary deformations, and meteorological processes.High-rate GPS data that belong to regional and global networks are often made available in near real-time through open-source projects, where the intention is to produce and disseminate these data to the largest possible community.
The aim of the present study is to demonstrate how some first-order information about a source model can be retrieved almost immediately after an earthquake, through a simple inversion of coseismic deformations obtained from publicly available near-real-time GPS solutions.More accurate estimates of source details can clearly be obtained using a collection of different datasets, e.g.teleseismic records, interferometric synthetic aperture radar (InSAR), local GPS and seismological data, rather than continuous GPS data alone.On the other hand, such datasets are not generally available in near real-time (and such datasets are often not publicly available at all).
The Scripps Orbit and Permanent Array Center (SOPAC) is an open-access archive of high-precision GPS data that provides, among other resources, daily solutions from about 800 continuous GPS sites that belong to more than 20 scientific networks around the World.In the present study, we used GPS daily solutions downloaded from SOPAC to infer first-order details about the source model of the 2010 Maule earthquake (Chile), through nonlinear inversion of the coseismic deformation field.Although the resulting dataset contains only two near-field offsets (from the CONZ and SANT sites), it can be shown to improve the characterization of the seismic source with respect to preliminary seismological models.

The 2010 Maule earthquake
On February 27, 2010, a devastating megathrust earthquake occurred offshore of Maule, in central Chile.On the moment magnitude scale (M W ), this 2010 Maule earthquake was estimated at 8.8 (United States Geological Survey, USGS; http://earthquake.usgs.gov/earthquakes/eqinthenews/), corresponding to a seismic moment release of 1.8 × 10 22 Nm; this therefore remains amongst the largest events of the instrumental age, and the second-largest recorded event in South America, after the great 1960 Chile earthquake.
The 2010 Maule earthquake was generated by the convergence of the Nazca plate downwards to the South American plate (Figure 1).The two plates are converging at a rate of 66 mm per year.The fault rupture extent has been estimated to be more than 100 km in width and to stretch nearly 500 km parallel to the coast of Chile (USGS, UNAVCO, http://www.unavco.org/support/event-response/eventresponse.html and Lay et al. [2010]).The rupture spread westward, northward and southward, and occurred

Subject classification:
Geodesy, Tectonophysics, GPS, Coseismic deformations, Source inversion, Maule earthquake.mostly offshore, generating a Pacific-wide tsunami that had devastating effects on the nearby coasts (National Oceanic Atmospheric Administration, http://nctr.pmel.noaa.gov/chile20100227/ and Saito et al. [2010]).The mid-latitude localization of the event maximized its effects on the dynamics of the Earth rotation, which induced perturbations that should have shortened the length of the day by 1.26 µs and to have moved the Earth figure axis by ~8 cm, according to current estimates (National Aeronautics and Space Administration, http://www.nasa.gov/topics/earth/features/earth-20100301.html).
After the 2004 M W = 9.3 Sumatra event, several studies found evidence of detectable coseismic offsets at GPS sites also on an extremely large scale (comparable with plate dimensions) [Banerjee et al. 2005, Vigny et al. 2005, Gahalaut et al. 2006].These far-field offsets have been successfully used to provide additional constraints to the source geometry and Earth rheological layering [Pollitz et al. 2006[Pollitz et al. , 2008]].Therefore the 2010 Maule earthquake is also likely to have left a detectable signature on GPS recordings, even on a wide regional scale.
Moreover, analysis of the GPS-derived coseismic deformation field induced by the Sumatra earthquake revealed that a consistent fraction of the energy was released aseismically [Banerjee et al. 2005].For the 2010 Maule earthquake, the present evidence is controversial: an analysis of normal-mode amplitudes has implied that 26% of the Global Centroid Moment Tensor solution (GCMT, http://www.globalcmt.org/)seismic moment was released aseismically [Tanimoto and Ji 2010], whilst using GPS and InSAR data, Pollitz et al. [2011] inferred a seismic moment that exceeds the GCMT estimate by 10%.On the other hand, joint inversion of high-rate GPS, InSAR and teleseismic data has been used to model the seismic source without invoking any aseismic energy contributions [Delouis et al. 2010].For this reason, we analyzed the SOPAC daily solutions of permanent GPS sites located in a vast region centered on the 2010 Maule earthquake epicenter, with a suitably designed MATLAB code that transforms the daily solutions into the International Terrestrial Reference Frame (ITRF2005) first, and then estimates the coseismic displacement field; this procedure was partly described by Devoti et al. [2010].Afterwards, our best estimate of the coseismic displacement field was used to infer the scalar seismic moment and the rake direction of the seismic source, through a nonlinear inversion procedure.

GPS analysis and displacement field
To extract the deformation field associated with the Maule earthquake, we analyzed the daily coordinates of 26 publicly available permanent GPS sites that belong to different networks, including: the International Global Navigation Satellite Systems Service (IGS), the Continuously Operating Reference Stations (CORS), and the Pacific GPS Facility (PGF).The GPS sites are located in a geographical area with longitudes between 0˚W and 160˚W, and latitudes between 75˚S and 30˚N.The coordinate time series were obtained from the daily quasi-observation files that are available from the SOPAC data center (ftp://garner.ucsd.edu/).The quasi-observation solutions contain coordinate estimates of nearly 200 sites of the IGS network, and the reference frame is only loosely constrained at the 1m level.We applied minimal inner constraints to the daily looseconstrained solutions, constraining translations, scale and rotations, to 1 mm.Each daily network is then transformed into the ITRF2005 reference frame, which is defined by 19 globally distributed sites that were extracted from the IGS cumulative solution [Dow et al. 2009] (ftp://igscb.jpl.nasa.gov/igscb/station/coord/IGS05.snx).These were chosen both on the basis of their data quality and their distance from the epicenter, which assured "noninfluence" from coseismic displacements associated with the Maule event.
The above procedure was applied to a dataset that included solutions from February 14, 2010 (GPS week 1571, day 0) to March 6, 2010 (GPS week 1573, day 6).More precisely, for each site, we stopped the analysis after the first two available daily solutions following the event in the GPS week 1573, to avoid any bias towards the instantaneous offset computation, due to a remarkable post-seismic drift that many stations presented, even in the first week following the earthquake (Figure 2).The sites located in the selected geographical region but without data in the proximity of the event (BRAZ, GLPS, VESL), were discarded from this analysis.
To measure the amount of coseismic displacement due to the seismic event, we estimated a common three-FIRST-ORDER SOURCE DETAILS FROM GPS dimensional offset at the epoch of February 27, 2010, in a leastsquares sense.This used the whole daily covariance matrices, and we evaluated the mean site positions in the 14 days before the earthquake and in the first two available days after the earthquake.The formal standard deviations associated with the estimates are likely to be underestimated, depending on the deviation from normality of the de-trended residuals [Williams et al. 2004].Figure 3 shows the offsets (blue vectors) and their associated errors (confidence ellipses) from the Maule earthquake for the North and East components.We only focused on the horizontal displacements, because the uncertainties of the vertical components are too large to allow any reliable interpretation.
Table 1 shows the preliminary estimated coseismic offsets and their associated errors for the North and East components.The maximum horizontal coseismic surface displacements were at CONZ (305.0 ±0.4 cm) and SANT (29.9 ±0.6 cm), the sites that were closest to the epicenter, both directed westward.Appreciable deformations that were still of the order of millimeters were evident even for the farfield sites that were included in the analysis (KOKB, HNLC, LPAL, GMAS).Except for the near-field sites (CONZ, SANT, LPGS), the relative uncertainties on the horizontal displacements ranged from 200% to 300%.Even though characterized by a lack of any public GPS station near the earthquake rupture zone, the coseismic deformation field obtained is consistent with the Nazca-South American plate dynamics [Haberland et al. 2009] and can be used to infer at least large-scale information on the source characteristics.

Inversion and source modeling
We inverted the residual deformation field associated with the Maule earthquake, with the aim of obtaining information about the seismic source from the geodetic data.The semi-analytical model used to compute the residual deformation field was widely discussed by [Piersanti et al. 1997, Soldati et al. 1998, Boschi et al. 2000], and it generally allows the computation of both coseismic and postseismic responses to a seismic dislocation for a spherical, incompressible, self-gravitating Earth with arbitrary rheological layerings [Melini et al. 2008, and references therein].In this case, we used a 45-layer structure with Preliminary Reference Earth Model (PREM)-averaged [Dziewonski and Anderson 1981] elastic constants; the corresponding density and rigidity profiles are shown in Supplementary Figure S1.Due to the poor near-field data coverage, the source geometry was fixed according to GCMT and USGS preliminary seismological estimates, and the only free parameters were the seismic moment and the rake angle (the position and extent of the fault plane is shown in the inset of Figure 3).Since our forward model assumes point-like dislocations, the rupture plane was discretized with 1,080 evenly discrete sources.For each elementary source, we computed the corresponding Green's Functions for the displacement field at the station locations for a unitary seismic moment, assuming both pure thrust and pure strike mechanisms and a uniform slip over the rupture plane, so the same slip parameters are assigned to all of the elementary sources.The deformation field for an arbitrary configuration of seismic moment and rake can therefore be easily obtained by a suitable combination of the elementary Green's Functions.
The optimal values of the seismic moment and rake are then found by minimizing the weighted chi-square between the model prediction and the observed displacement.The minimization was performed with the MINUIT package [James 1998], which implements a variable-metric method with inexact line search, a stable metric updating scheme, and a positive-definiteness check [Fletcher 1970].Radial displacements were not used in the inversion because of the lower resolution of the GPS system for vertical movements.The modeled field that results from the inversion is shown in Figure 3, along with our best estimate of the observed field, as described in the previous section.The normalized FIRST-ORDER SOURCE DETAILS FROM GPS chi-square between the data and the model is 84.2, and the wrms of residuals is 31 mm.For the South American sites, we obtain a good agreement on the direction of the deformation vectors (within the error ellipses) except at UNSA, which might be located along a nodal line of our model.The vector lengths at some sites appear to be overestimated by our model, by up to a factor of two (UNSA, LPGS and SANT).While the overestimation on the SANT site might be explained by edge effects, on the other sites, the incompressibility approximation of our analytical model might have a role, as in compressible models a larger fraction of the energy is stored in elastic loading, which yields a faster decay of the deformation with distance [Nostro et al. 1999].The agreement in the Central American region and the Antarctic peninsula is less satisfactory; however, we note that the modeled displacements in these regions are coherent, whilst the observed data are characterized by large variability, both in amplitude and direction, on small spatial scales.On the very far-field sites (Pacific Ocean and Northwestern Africa), the agreement is good, even if these sites are less significant, due to the large relative errors of the modeled vectors.
The best-fitting source model is characterized by a seismic moment M 0 =2.82×10 22 Nm and a rake angle m= 118˚.The rake angle that results from the inversion corresponds to a thrust mechanism with a small right-lateral component, which is consistent with the GCMT solution (m = 112˚) and the W-phase inversion (m = 109˚) [Lay et al. 2010].The seismic moment, which corresponds to a moment magnitude M W = 8.9, is about 60% greater than its (GCMT, USGS) seismological estimate M 0 = 1.8 × 10 22 , which is equivalent to the (well-known) moment magnitude M W = 8.8.
The resulting excess of energy released from geodetic inversion with respect to purely seismological methods has already been seen for the 2004 Sumatra earthquake [Banerjee et al. 2005, Lay et al. 2005], and it was interpreted as an indication of a large, very low frequency energy release in giant subduction events, even if Hearn and Burgmann [2005] suggested that the modelistic artifacts connected with different rheological layerings used in seismological and geodetic inversions might have a role too.For the 2010 Maule earthquake, the seismic moment estimates vary from 1.8 to 2.6 × 10 22 Nm, depending on the data and modelistic approaches used (see Supplementary Table S1 of Pollitz et al. [2011] for a summary), and they indicate the possibility of a large aseismic contribution to the energy release also for the 2010 Maule earthquake.Our estimate of M 0 is at the upper bounds of the range of published values; while it may be overestimated due to the lack of detailed near-field coverage in our dataset, it is in qualitative agreement with indications coming from the literature.However, we also note that our forward model describes the seismic source in terms of a double couple, and therefore it provides a natural estimate of the seismic moment with no need for any slipmoment conversion that might be biased by local rigidity assumptions.
To assess the resolution of the model parameters and to ascertain the presence of any trade-offs between them, we carried out a Montecarlo analysis.Here, we generated 100,000 synthetic datasets, each of which was obtained by adding random realizations of data noise to the original displacement field.Afterwards, for each synthetic dataset, we computed the best-fit model parameters.From the distribution of these synthetic models, we can establish the effects of experimental uncertainties on the inverted model parameters.Supplementary Figure S2 shows the a-posteriori distribution of the seismic moment and the rake angle, and the correlation between them.As can be seen from Supplementary Figure S2, the distribution of parameter models is well approximated by a Gaussian function, and there is no correlation.By fitting the distributions with a Gaussian curve, we obtained an estimate of the error of the parameters, which turns out to be M 0 = (2.818±0.034) × 10 22 Nm and m = (117.85±0.77)˚.

Conclusions
In the present study, we show how it is possible to use publicly available real-time data from worldwide networks to quickly infer first-order information that can help to better constrain the source parameters of a large earthquake.
Here, we used publicly available records of permanent GPS stations belonging to the IGS, CORS and PGF networks, to define the static deformation field associated with the coseismic effects of the M W = 8.8 2010 Maule earthquake, by analyzing short time series, including two weeks of daily solutions.The resulting estimate of the coseismic field shows a maximum horizontal deformation of over 3 m at the CONZ site, which is located on the hanging wall of the fault plane.Appreciable displacements of the order of millimeters are even evident for the far-field sites located at over 1,000 km from the epicenter.
The deformation field estimated from these GPS data was used to infer information about the seismic source with a nonlinear inversion.Due to the poor near-field data coverage, we used a fixed-geometry source model, using information from GCMT and USGS seismological analyses, and inverted only for seismic moment and rake angle.The modeled field is in agreement (within the associated errors) with the observations, with a normalized chi-square of 84.2 and a wrms of residuals of 31 mm.For a few sites, the model overestimates the data, most likely because of the incompressibility approximation assumed by the model.The best-fit rake angle is consistent with GCMT and USGS seismological estimates and with Lay et al. [2010], while the scalar moment (M 0 = 2.82 × 10 22 ) is about 1.6 times the seismological estimate, which corresponds to a moment magnitude M W = 8.9.The overestimation of the seismic moment using geodetic inversions has already been seen for the Sumatra earthquake [Banerjee et al. 2005, Lay et al. 2005], and this has been interpreted in terms of a large aseismic energy release in giant subduction events.For the 2010 Maule earthquake, a similar indication comes from the analysis of normal-mode amplitudes, even if the situation is not yet clear [Delouis et al. 2010, Tanimoto and Ji 2010, Pollitz et al. 2011].
The timeliness and public availability of high-rate GPS data make them an ideal tool for quick post-earthquake analysis.Of course, in the present study, we discussed an approach that can be used only if significant offsets have been recorded by at least a few local or regional-scale sites.As a result, this kind of analysis will be possible only for earthquakes above a magnitude threshold that depends on the (spatially variable) density of open-access global GPS networks.At the moment, we reasonably estimate that this analysis can be performed for earthquakes above a magnitude of 8.0.Over the past decade, more and more regional and global networks of continuously operating GPS ground stations have opened up the access to their data; if this trend is maintained in the future, with an increasing production and dissemination of open-access data, near realtime GPS analysis might become a viable tool for rapid characterization of even smaller earthquakes.

Figure 1 .
Figure 1.Regional seismotectonic framework of the 2010 Maule earthquake.The gray line marks the major plate boundaries [Bird 2003].The yellow star and the red beach ball show the epicenter location and the focal mechanism of the Harvard CMT solution for the main shock.The focal mechanisms from the CMT catalog (since 1976) are also reported.The red box marks the surface projection of the USGS rupture geometry estimate (Gavin Hayes, USGS).The lower-right inset shows the CMT seismicity cross-section along profile A-A´ (each tick mark corresponds to 50 km).

Figure 2 .
Figure 2. Time series for the North, East and vertical components of four GPS sites, as indicated.Error bars correspond to the 1-v formal errors.The vertical red line represents t = 2010.02.27 (February 27, 2010).

Figure 3 .
Figure3.Observed GPS coseismic displacements (blue vectors) with 68% confidence ellipses and modeled (red vector) displacements of the Maule earthquake, with the exception of the SANT and CONZ sites.The yellow star indicates the earthquake epicenter.The BRAZ, GLPS and VESL sites are included on the map, even though their data were not included in the analysis.The lower left inset shows an enlarged view of the near-field area with the horizontal displacements detected by the SANT and CONZ sites, and the position and extent of the modeled fault plane.The displacement at SANT is magnified by a factor of 200.