# Retrieving earthquake signature in grace gravity solutions

O. de Viron, I. Panet, V. Mikhailov, M. van Camp, M. Diament

► To cite this version:

O. de Viron, I. Panet, V. Mikhailov, M. van Camp, M. Diament. Retrieving earthquake signature in grace gravity solutions. Geophysical Journal International, 2008, 174, pp.14-20. 10.1111/j.1365-246X.2008.03807.x. insu-03603718

HAL Id: insu-03603718

https://insu.hal.science/insu-03603718

Submitted on 10 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Retrieving earthquake signature in grace gravity solutions

O. de Viron,¹,² I. Panet,¹ V. Mikhailov,¹,³ M. Van Camp⁴ and M. Diament¹

¹Institut de Physique du Globe de Paris (assoc. CNRS), France. E-mail: deviron@ipgp.jussieu.fr

²Université Denis Diderot (assoc. CNRS), Paris, France

³Institute of physics of the Earth RAS, Moscow 123995, Russia

⁴Royal Observatory, Belgium

Accepted 2008 March 27. Received 2008 March 26; in original form 2007 December 6

## SUMMARY

The GRACE satellites have been orbiting the Earth since 2002, monitoring the time variable gravity field. Some of the observed fluctuations are due to geodynamic causes, but they are often hidden in the complex signal, composed of hydrology, ocean, atmosphere, and geodynamics, the signal of geodynamic origin being usually the smallest. In addition, dealiasing residuals and noise make the separation of the signal from the different causes more difficult. We proposed a method based on the Empirical Orthogonal Function decomposition to extract the signal of physical origin, under the hypothesis that the physical signal is spatially more consistent than the noise and aliasing incomplete correction. We used synthetic geoid variations associated with earthquakes located at nearly 2000 positions at the Earth surface, based on several examples of large actual subduction events. We show that, with the present day accuracy, we can retrieve the geoid variations associated with more than 98 per cent of the earthquakes of magnitude 9 or above, around 60 per cent for magnitude 8.8, 40 per cent for magnitude 8.6 and 33 per cent for magnitude 8.3. Some events, with the right properties and location, can be detected with magnitude as low as 8. We then applied the method to the GRACE solutions, and retrieved the Hokkaido event (2003) and the Sumatra event (2004), which is in agreement with the retrieval rates mentioned here above.

**Key words:** Time series analysis; Time variable gravity; Earthquake ground motions.

## 1 INTRODUCTION

Since 2002, the GRACE satellites have been orbiting the Earth (see Tapley et al. 2004), monitoring the gravity field of the solid Earth, atmosphere, ocean, and hydrology. As the gravity is an integrated quantity, the recorded signal includes effects caused by mass redistribution inside and outside the solid Earth. The most challenging part of the GRACE data analysis is to attribute the observation to the right causative sources. This task is made even more difficult by the larger than expected noise level in the data, at least partly caused by errors in the dealiasing products.

The variation of mass distribution inside the Earth is dominated by the climate signal, mostly linked to the water exchange between the different parts of the world. The signal of geodynamic origin is a small contribution to the total mass redistribution. If we want to use the gravity data to improve our knowledge of the geodynamic processes, for instance to complement seismic data about earthquake events, it is necessary to separate this signal from the climate one, without entering too much *a priori* knowledge of the event. Previous studies (Han et al. 2006; Chen et al. 2007; Ogawa & Heki 2007; Panet et al. 2007) were based on the hypothesis that the spectral signature (spatial or temporal) of the earthquake was different enough from the climate signal to allow separation. Mikhailov et al. (2004)

proposed a statistical method, requiring an *a priori* knowledge of the earthquake gravity signature. The method proposed in this study allows to use the full spatio-temporal signature, under the assumption that the large earthquake signal will have a spatially coherent signature, with a Heaviside-like temporal behaviour.

As shown by Barnes (1997), the observable gravity potential signature of an earthquake at the surface would be mainly associated with the vertical displacement of the Earth's surface. For the very large earthquakes, in relatively short series of spherical harmonic expansion, effect of rock density changes during deformation becomes more important. As a result, effect of dilatation should also be taken into account (Gross & Chao 2001; Han et al. 2006; Panet et al. 2007). Mikhailov et al. (2004) have shown that the signal from major subduction earthquakes should be above the GRACE observation precision level, and that a statistical method should allow to test several focal mechanisms against each other and against pure noise in order to see if the signature significantly differs from noise, and if a model is significantly more probable than the others. This method could then be used to discriminate between possible fault plane models as well.

Retrieving earthquake signal from global gravity data is important from three points of view: (1) coverage of those data is globally uniform, allowing to study earthquakes whose epicentral area is
situated undersea, which is not so easy by geodetic and other surface techniques; (2) the vertical displacement which can be retrieved from GRACE gravity data is the least precise component of the GPS observed motion and (3) the gravity data provide information on the density variations at different depths.

The signature of the Sumatra megathrust event has already been extracted from the GRACE data by several workers (Han et al. 2006; Chen et al. 2007; Ogawa & Heki 2007; Panet et al. 2007).

de Viron et al. (2006) proposed a method based on the Empirical Orthogonal Function (EOF) decomposition to extract meaningful signal from the GRACE data. In this study, we apply this method to investigate earthquake signatures, under the hypothesis that the gravity signals associated with them are one of the dominant EOF modes in the area. The main advantage of this method is that it does require less a priori information on the event, which makes the detection more objective. Indeed, the other methods used up to now consider that we know when and where the event occurred, and use that information to separate the earthquake signature. The EOF-based method only use the information of location and timing as a test of the robustness of the result.

In this study, we first use synthetic data to systematically test which are the earthquakes producing a retrievable gravity signal.

In Section 2, the different GRACE time-series are presented, discussing how they were pre-processed for this study. Synthetic and real earthquake gravity data are also described. In Section 3, the precision of the GRACE data with respect to the expected signature of the earthquakes is evaluated. In Section 4, the proposed detection procedure is described and synthetic data are used to systematically test which earthquakes can produce a retrievable gravity signal. Then, in Section 5, the detection procedure is applied on real earthquakes contemporary with the GRACE mission and the results are compared to those obtained with synthetic data. Section 6 is devoted to discussions and conclusions.

## 2 DATA USED AND THEIR PRE-PROCESSING

### 2.1 Gravity data

In this study, we used the CNES release of the GRACE data (see Biancale et al. 2005). The data are given as either a grid (1° × 1°) or a spherical harmonic expansion (up to degree 50) of the geoid height. The data are given every 10 d, as a monthly mean centred on those 10 d. The last release includes data from July 2002 up to June 2007. They have been corrected for the atmospheric effect on the geoid, as predicted by the European Center for Medium-range Weather Forecast (ECMWF) analysis model, for the ocean effect as predicted by MOG2D (Carrère & Lyard (2003)). The ocean tides effect has been corrected for, using outputs from the FES-2004 model

(Lefèvre et al. 2000). Consequently, only hydrology, geodynamics, and the difference between the barotropic and the baroclinic ocean signal should remain in the data. Nevertheless, according to for example, Schmidt et al. (2006), the time variation in the GRACE geoid data exhibits meridionally oriented stripes of notable amplitude, presently attributed to aliasing of the miscorrected high-frequency signal originating in the Earth fluid layers (atmosphere and ocean).

We estimated a seasonal cycle (annual + semi-annual) by least-square fit at every gridpoint, and removed it.

### 2.2 Simulation data

To systematically investigate the parameters of the earthquakes that can be retrieved from GRACE data, we generated synthetic earthquake gravity signals that we added to the original GRACE geoids. In order to limit our investigation to realistic seismic events affecting the geoid, we used the characteristics of five past earthquakes of large magnitude, whose vertical deformation reaches the surface. Information about them is provided in Table 1. The geoid variations caused by an earthquake result from the displacements of density interfaces, the main one being the Earth's surface, and from the density variations associated with the Earth's deformation. This last contribution is particularly important at GRACE spatial scales (Han et al. 2006). We thus computed the geoid change associated with an earthquake as follows. First, we estimated the vertical coseis-mic surface displacement using an elastic fault plane model for the investigated earthquake. Even if more precise models can be found, displacements linked to the earthquake motions are indeed efficiently described using the analytical solution for finite, rectangular sources in an elastic half-space (Savage & Hastie 1966; Okada 1985). The obtained vertical displacement is then projected onto a sphere. This is a reasonable approximation, since the effect of the Earth's sphericity is small in the area of main deformation (Sun & Okubo 1993), and negligible in the far field (see supporting on-line material to Banerjee et al. 2005). We then computed the geoid effect of crustal and mantle density changes using an analytical solution for the internal deformation in an elastic half-space as a trace of the strain tensor (Okada 1992). Whereas superficial layers are affected by rock dilatation, compression occurs in the deeper layers and partially compensate the superficial contributions. Geoid effect of density changes depends on how compliant the rocks are. In this study, we assumed a Poisson ratio equal to 0.3, a realistic value for the lithosphere and upper mantle rocks according to the PREM model (Dziewonski & Anderson 1981). The surface displacement and density variations can finally be converted into loads at different depths, from which we compute the geoid variation according to Wahr et al. (1998). The approach is also explained in Mikhailov et al. (2004).

<table><caption>**Table 1.** Typical earthquakes used for the synthetic study, and the references where the characteristics where found. In addition, the web site of the NEIC (National Earthquake Information Center) has been used.</caption><thead><tr><th>Date</th><th>Latitude</th><th>Longitude</th><th>Dip</th><th>Slip</th><th>M<sub>w</sub></th><th>Location</th><th>Reference</th></tr></thead><tbody><tr><td>1960 May 22</td><td>-74.5°</td><td>-39.5°</td><td>45°</td><td>90°</td><td>9.5</td><td>Chile</td><td>(Plafker & Savage 1970)</td></tr><tr><td>2004 December 26</td><td>95.9°</td><td>3.3°</td><td>Shallow: 11-18°<br>Deep: 35°</td><td>90-140°<br>90-140°</td><td>9.3</td><td>Sumatra</td><td>(Banerjee et al. 2005)</td></tr><tr><td>2005 March 28</td><td>97.0°</td><td>2.1°</td><td>15°</td><td>90-120°</td><td>8.6</td><td>Nias</td><td>(Banerjee et al. 2007)</td></tr><tr><td>1906 January 31</td><td>-81.5°</td><td>1.0°</td><td>25°</td><td>129°</td><td>8.8</td><td>Ecuador</td><td>(Gutenberg & Richter 1959),<br>(Kanamori & McNally 1982)</td></tr><tr><td>2006 November 15</td><td>153.2°</td><td>46.6°</td><td>15°</td><td>92°</td><td>8.3</td><td>Kuril</td><td>(Tanioka & Hasegawa 2007)</td></tr></tbody></table>
## 2.3 Expected data precision

Before testing every single possible configuration, it is wise to start by an order of magnitude estimation. We can expect to retrieve an earthquake if the associated geoid signal is characterized by a not too low signal to noise ratio. As a first step, we tried to get estimation of the data precision. According to R. Biancale (personal communication, 2007), the CNES GRACE solution are precise at the level of a few $10^{-4}$ m of geoid topography, for length scale of about 1000 km.

It is difficult to assess the precision of high quality data, when there is nothing else to compare with. At present, several solution using the same GRACE data have been provided by several groups, and we can try to get precision estimate from their intercomparison. In general, when three independent measurements of the same quantity are available, it is possible to estimate the noise level of each measurement. Let us assume that each data $[x_i(t), i = 1, 2, 3]$ is composed of a true signal $[X(t)]$ and some errors, which are different for each measurement $[\epsilon_i(t)]$. If those errors are independent, we can link the variance of the difference between two series and the variances of their individual errors:

$$ \text{var}(x_i - x_j) = \text{var}(\epsilon_i) + \text{var}(\epsilon_j). \quad (1) $$

From three measurements, it is possible to compute three differences of variances, and the obtained equation system can then be solved for the estimation of the variance of the noise in each measurement.

When the measurements are not independent, as it is the case in our study, this method does not hold, as a covariance term appears in the equation:

$$ \text{var}(x_i - x_j) = \text{var}(\epsilon_i) + \text{var}(\epsilon_j) - 2\text{cov}(\epsilon_i, \epsilon_j). \quad (2) $$

In the three solutions used, the GFZ and CSR solutions are very similar, and the CNES solution differs, because a regularisation procedure has been applied by CNES, and not by the other two institutions (Bettadpur 2006; Flechtner 2006). As explained above,

the GRACE solutions present strike-like noise in the meridional direction; the regularisation allows to reduce them in the solutions. In addition, data from the LAGEOS satellite have been used by CNES, but this should not affect much the final Earth gravity models, except for the degree two of the geoid spherical harmonic expansion. As the GFZ and CSR solution are very similar to each other, their covariance will be large, and it will decrease the variance of their difference. The similarities between any of these solutions with the CNES data is much smaller, and the variances of their difference will be larger. Consequently, when solving the three-cornered hat equation system, we will get a larger number for the CNES solution, simply because the solution is more independent. On the other hand, the results for the other two solutions gives us information on how close the solutions are to each other. The conclusion we can reach from this study are the following: (1) All the solutions cannot be all more precise that the smallest of the three error estimate, and we can use this as a lower boundary for the precision and (2) the results of the three-cornered hat system provide us with information on the level of the independent signal from each solution.

In Fig. 1, we show the results of the three-cornered hat applied to each data point smoothed on a $7.5^\circ \times 7.5^\circ$ window to emphasize the large scale structures. It can be observed that the differences between the solutions do not coincide with some geographic structure or tectonic units, but are stronger at high latitude. From these results, we estimate the noise level to be between 0.1 and 1.5 mm of geoid: it is unclear which signal is right using the three-cornered hat method only, but all of the solutions cannot be more accurate than 1.5 mm, as they differ at that level. Nevertheless, considering that part of the difference comes from non-geophysical sources which are considerably reduced in the CNES solutions, this estimate is probably pessimistic for the CNES solutions. The estimate 0.1 mm is probably too optimistic, and we expect the true value to lie in between. Note that the estimated precision provided by R. Biancale is close to and included in the interval between the three-cornered hat estimation for the two other solutions. These results

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Three-cornered hat standard deviation of the independent signal evaluation for the CSR, GFZ and CNES time variable geoid, as a function of the position. The results have been smoothed on a $7.5^\circ \times 7.5^\circ$ window to emphasize the large-scale structures.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Standard deviation of the geoid deformation as a function of the magnitude of the simulated earthquake. The horizontal lines indicate the different estimations of the data precision (TCH stands for three-cornered hat).</figcaption></figure>

indicate that a 0.3 mm estimation of the data precision is probably reasonable. As we are using EOF analysis, it is not unreasonable to hope that there will be some error cancellation by averaging. If we assume that an EOF is a mean like operation over the time-series, we can hope for an error reduction at most of the order of the square root of the number of independent data points (about 60 sets of monthly independent data), which means an error level around 0.04–0.05 mm on the EOF. For that reason, it is useless to look for earthquake that would generate a geoid variation smaller than that amount.

Using those estimates of the data precision, we can get a first idea of the retrievable and non retrievable events. The size of the geoid signature of an earthquake in the GRACE data depends mostly on the normal component of displacement on a fault plant, that is, on vertical displacement at the surface associated with the event, and this displacement can be estimated from the earthquake parameters. We generated a large set of synthetic shallow subduction earthquakes of magnitude ranging between 6.0 and 9.5, and we estimated the geoid deformation. This estimates are independent on the location and time of the events. The minimum effect is close to zero, for strike-slip events, as expected. In this study, we focus on the megathrust events, which are the most likely to generate an observable signal (see Fig. 2). For near-surface events, the geoid change is a monotonic function of the magnitude. Fig. 2 shows the interval of expected effect from the synthetic earthquakes, as a function of the magnitude, estimated as the standard deviation of the associated geoid change on an 20° × 20° area. Consequently, it is very unlikely to detect any earthquake of magnitude smaller than 8.0, simply because the signal is not large enough. For magnitudes 8.5 and larger, we should be able to retrieve most of the earthquakes associated with large megathrust events.

## 2.4 Real earthquakes

After showing the efficiency of the method on synthetic data in Section 4, we apply it on real earthquakes in Section 5, in order to investigate which ones can be retrieved. The times and locations of the tested earthquakes are from the USGS earthquake database, NEIC, available through the website http://earthquake.usgs.gov/.

## 3 TEST USING SYNTHETIC DATA

Using the focal mechanisms of real large subduction earthquakes, we computed the geoid change associated with the event. We then performed a systematic test of detection. To this end, we changed the position, the epoch and the strike of the studied earthquake, and added the rotated and translated signal to the observed GRACE geoids. We tried 1848 different positions at the Earth surface (regularly spaced at 5° intervals, of latitude from −82° to 82°). For each position, we tested each earthquake with 16 different strike angles ($k\pi/8$) and four different occurrence times covering the whole period, except for the four first and four last months of the time-series. As the EOF method is based on the variance explained, an event occurring very close to the beginning or to the end of the data set will not be detected, except for a very large event.

The EOF decomposition (see Preisendorfer 1988) allows to represent a set of time-series as a set of empirical orthogonal functions (which are spatial patterns) associated with time-series. The main interest of this method is to represent most of the variance of the set of time-series with only a few of those EOF. In this study, we focus on the first five EOFs, which represent about 90 per cent of the signal variance, and analyse them, in order to determine if one of them is likely to represent the geoid change associated with the earthquake. This method is explained in details for instance in (Preisendorfer 1988). The EOF decomposition represents a space–time data set $[x_i(t_j), i = 1, ..., n, j = 1, ..., m]$ in terms of a given number $N$ of variability modes, each of which is a time-series $A_k(t)$ and a geographical distribution $X_k(i)$, with

$$x_i(t_j) = \sum_{k=1}^{N} A_k(t_j)X_k(i). \quad (3)$$

These modes are obtained from the eigenvalue decomposition of the covariance matrix $R = F^T F$, with the matrix $F_{ij} = x_i(t_j)$. The eigenvalues represent the variance explained by the mode of variability, and the eigenvectors, often called EOF as they are orthogonal, represent the space distribution of the modes. The modes are sorted by decreasing eigenvalue, so that the first mode explains more variance than any other. The time variability associated with each EOF can be retrieved by projecting the matrix $F$ on the eigenvector.

We propose a very simple detection method, which have been chosen after a set of test, to optimize an automatic detection, limiting false alarm and missed detection in presence of noise. We define a region around the event of twice the magnitude in longitude and in latitude (for instance, a magnitude 8 event will be looked for in a 16° × 16° box). For magnitude 8 events, the area size is about four times-the typical rupture length. We make an EOF decomposition of the geoid time-series inside the box, and analyse the EOF spatial structure and time-series. We fit a Heaviside function on the EOF time-series at every possible time in order to find at which time it explains the most variance. We compute the correlation between the spatial pattern of the EOF and the signal we added to the geoid. We consider that we retrieved the earthquake if (1) the best time for a jump is at less than 3 months from the ‘true time’, (2) the Heaviside function, fitted on the observation, explains at least 20 per cent of the variance of the EOF time-series and (3) the synthetic and the retrieved patterns have a correlation coefficient larger than 0.5. If these three conditions are verified together for at least one of the EOF mode, we consider that we detected the earthquake.

Of course, for synthetic event, all the information is available, allowing to use the whole set of criteria. If we want to estimate the vertical displacement from gravity data without a priori information on the ground displacement, we cannot use the condition (3).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Results of the synthetic test, in per cent of retrieved earthquake, as a function of the azimuth of the perpendicular to the largest slope direction in the added signal. The mean retrieval rate is given in the legend.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** EOF time-series and spatial pattern associated with the Hokkaido event. The time of the event is shown by a grey line on the left-hand panel.</figcaption></figure>

## 3.1 Null hypothesis test

We need to make sure that this test is efficient enough to ensure that what we detect is not simply random fluctuations of the geoid which mimics the earthquake signal. To this end, we tried to detect each event, as described here above, but with no synthetic signal added. The results are fairly good. The detection rate is smaller than 3.5 per cent for the eight tested synthetic 'events', and more than 5 out of them were falsely detected in less than 0.5 per cent of the tested locations. We accordingly put an error bar of 3.5 per cent on our detection rate, in order to be on the safe side.

## 3.2 Results of the synthetic test

The results of the synthetic test are summarized in Fig. 3 as a function of magnitude and strike angle. We observed a low decrease of the detection rate when the fault plane is oriented north-south, may be due to the dealiasing errors having the same orientation.

We then use the Ecuador event (with a total detection rate close to 55 per cent) to study which other parameters can affect the detection probability. The synthetic tests show a geographical effect: Australian, Eurasian and South American event are not detected as often as the other ones, North American and African events being the best detected. We have presently no explanation for this effect. Results for oceans areas are in between, but close to the best continent area. This might be due to the absence of signal from hydrological origin that could mask the earthquake signal.

Averaging on the whole Earth surface, we get a slight raise of detection rate for events occurring in the middle of the data interval, but the difference is not significative according to a $\chi^2$ test (at least

for events not too close to the 4-month limits of the observation epoch).

# 4 RESULTS ON ACTUAL EARTHQUAKES

The characteristics of the largest earthquakes that occurred during the GRACE era are summarized in Table 2. The first and third events (Denali and Macquarie Island earthquakes) were not detected according to our tests. Indeed, for those strike-slip events, the vertical motion is very small and, consequently, so is the gravity signal.

When applying the method to real earthquakes detection, we do not keep the last detection criterium. Indeed, using the correlation between the retrieved signal and a predicted one would introduce a strong *a priori* in the detection, whereas we want to keep the method as objective as possible. On the other hand, it is always possible to compute *a posteriori* the correlations between the EOF retrieved signal (based on criteria 1 and 2), and candidates rupture models, in order to determine which one is the most realist.

The detection of the Hokkaido subduction event succeeded the detection tests. We do have a first EOF mode with a jump at the right time explaining the most possible variance (close to 54 per cent). Considering the degree of smoothing of the GRACE data, the spatial pattern is consistent with the coseismic deformation determined from GPS data and waveform inversion (see Miura *et al.* 2004, and references there in). The results are shown on Fig. 4.

We get a very clean signal for the Sumatra-Andaman earthquakes (about 77 per cent of variance explained). It is shown on Fig. 5. The

<table><caption>**Table 2.** Large events occurred during the GRACE era. Hokkaido according to Tanioka and Hasegawa (2007); Alaska 2002 according to Wright *et al.* (2004); Macquarie according to Robinson & Sandron (2006).</caption><thead><tr><th>Date</th><th>Latitude</th><th>Longitude</th><th>Dip</th><th>Slip</th><th>Depth</th><th>M<sub>w</sub></th><th>Location</th></tr></thead><tbody><tr><td>2002 March 11</td><td>63.52°</td><td>-147.44°</td><td>90°</td><td>180°</td><td>4</td><td>8.50</td><td>Denali Fault, Alaska</td></tr><tr><td>2003 September 25</td><td>41.81°</td><td>143.91°</td><td>20°</td><td>109°</td><td>27</td><td>8.30</td><td>Hokkaido, Japan</td></tr><tr><td>2004 December 23</td><td>-49.31°</td><td>161.35°</td><td>61°</td><td>1°</td><td>10</td><td>8.10</td><td>Macquarie Island</td></tr><tr><td>2004 December 26</td><td>3.30°</td><td>95.98°</td><td>Shallow: 11-18°<br>Deep: 35°</td><td>90-140°</td><td>30</td><td>9.00</td><td>Sumatra-Andaman</td></tr><tr><td>2005 March 28</td><td>2.09°</td><td>97.11°</td><td>15°</td><td>90-120°</td><td>30</td><td>8.60</td><td>Nias (Sumatra)</td></tr><tr><td>2006 November 15</td><td>46.59°</td><td>153.27°</td><td>15°</td><td>92°</td><td>10</td><td>8.30</td><td>Kuril Island</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** EOF time-series and spatial pattern associated with the Sumatra event. The time of the event is shown by a grey line on the left-hand panel.</figcaption></figure>

geographical pattern is consistent with those of Han *et al.* (2006) or Panet *et al.* (2007), for instance.

We do not get a detection for the Nias event. It may be due to the mixing with the Sumatra event signature: the two events occurred in the same area, and it is possible that the EOF has difficulties to separate the Nias event.

With the data set ending at the end of 2006, the Kuril event is too close to the end of the time-series we used to generate enough variance to be detected.

In conclusion, we detected the Sumatra-Andaman earthquake, we probably detected the Hokkaido event, and it is possible that, as Panet *et al.* (2007), we get signal from the Nias event mixed in the Sumatra earthquake EOF-mode. Considering the probability of detection and the characteristics of the earthquakes, the results are fairly good, as we detected what was likely (or even possible) to be detected.

## 5 DISCUSSIONS AND CONCLUSIONS

In this paper, we study the possibility to detect and characterize large earthquakes using the GRACE data. Only earthquakes associated with large vertical displacements at long wavelength generate observable signal in the GRACE data. Considering the present data quality, we cannot detect any event of magnitude smaller than 8, as they would generate signature in the geoid smaller than the noise level. The associated vertical motion needs to be at the metre level, which is obviously within the detection range of other techniques (seimology, geology, geodesy, etc.). Nevertheless, (1) we can hope that time variable gravity will improve its quality and resolution in the next few years and (2) gravity from space allows to study in areas where field work is difficult or other techniques inoperative, as undersea for example.

Considering that an earthquake will generate a jump in the geopotential variations, that would be very different from any other geophysical signature, we propose to use the EOF decomposition to separate the geopotential variations associated with the earthquake from the noise and geopotential from other sources. The method is very simple: we perform an EOF decomposition of the geopotential variation around the event location. Then, in the EOF modes, we

look for a mode with a jump at the time of the event, explaining a large part of the variance of the EOF time-series.

We tested this method on synthetic cases. We showed that it was promising on events associated with large vertical displacements. The rate of false alarm is very small, and the detection is above 30 per cent for events of magnitude 8.1, above 50 per cent for events of magnitude 8.8 events, and more than 95 per cent for magnitude 9 events.

Using the locations and occurrence times of events in the GRACE era, we tested the method, and we were able to detect two events, the Sumatra 2004 and the Hokkaido subduction earthquakes. Considering the characteristics of all the tested earthquakes and the probability of detection deduced from the synthetic tests, these results are reasonable.

The interpretation of satellite gravity data only is not easy, for the following reasons: (1) although the EOF decomposition is efficient to extract signal out of noisy data, both the time-series and the geographical pattern contain residual noise components and (2) the low resolution of the GRACE data smoothes away an important part of the signal. We should also keep in mind that satellite gravity provides a different information from seismic observations, since it includes both seismic and aseismic motions, and from GPS measurements, since it also includes effect of density changes at depth in addition to surface displacements effects. The combination of all these complementary information should lead to a better understanding of the geodynamic processes at stake.

Finally, space gravity is a tool that should not be forgotten in geodynamics study, as it allows to get information with a regular time and space distribution. We can expect an increase of the data precision and of the resolution in the next few years, which would allow better resolved geoid signal, allowing to complement the motion observed by the other geodetic techniques.

## ACKNOWLEDGMENTS

The authors are grateful to Richard Biancale for useful discussions and clarifications. The very constructive reviews helped to make this paper more understandable. This study was supported by CNES through the TOSCA ('Solid Earth, Ocean, continent surfaces, atmosphere' Terre solide, ocean, surfaces continentals, atmosphere) committee. This paper is the Institut de Physique du Globe de Paris (IPGP) contribution #2353.

## REFERENCES

* Banerjee P., Pollitz, F. & Bürgmann, R., 2005. Size and duration of the great 2004 Sumatra-Andaman earthquake from far-field static offsets, *Science*, **308**, 1769–1772.

* Banerjee, P., Pollitz, F., Nagarajan, B. & Bürgmann, R., 2007. Co-seismic slip distribution of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Nias earthquakes from GPS static offsets, *Bull. seism. Soc. Am.*, **97**(1A), S86–S102.

* Barnes, D.F., 1997. Gravity changes during the 26 years following the 1964 Alaskan earthquake, Geological studies in Alaska. *USGS professional paper*, **1614**, 115–122.

* Bettadpur, S., 2006. Release notes for CSR RL01 GRACE L2 products, May 1st 2006.

* Biancale, R., Lemoine, J.-M., Balmino, G., Bruinsma, S., Perosanz, F., Marty, J.-C., Loyer, S. & Gégout, P., 2005. 3 years of geoid variations from GRACE and LAGEOS data at 10-day intervals over the period from July 29th, 2002 to March 24th, Data CD, CNES/GRGS product.
Carrère L. & Lyard, F., 2003. Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing—comparisons with observations, *Geophys. Res. Lett.*, **30**(6), 1275, doi:10.1029/2002GL016473.

Chen, J.L., Wilson, C.R., Tapley, B.D. & Grand, S., 2007. GRACE detects co-seismic and postseismic deformation from the Sumatra-Andaman earthquake, *Geophys. Res. Lett.*, **34**(13), L13302, doi:10.1029/2007GL030356.

de Viron, O., Panet, I. & Diament, M., 2006. Extracting low frequency climate signal from GRACE data, *Earth Discuss.*, **1**, 21–36.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary Reference Earth Model (PREM), *Phys. Earth planet. Inter.*, **25**, 297–356.

Flechtner, F., 2006. Release notes for GFZ RL03 GRACE L2 products, May 9th 2006.

Gray, J.E. & Allan, D.W., 1974. A method for estimating the frequency stability of an individual oscillator, *Proc. 28th Freq. Contr. Symp.*, 243–246.

Gross, R. & Chao, B., 2001. The gravitational signature of earthquakes, in *Gravity, Geoid and Geodynamics 2000, IAG Symposia*, Vol. 123, pp. 205–210, Springer-Verlag, New-York.

Gutenberg, B. & Richter, C.F., 1959. *Seismicity of the Earth*, Princeton University Press, Princeton, New Jersey, 310 pp.

Han, S.C., Shum, C.K., Bevis, M., Ji, C., & Kuo, C.Y., 2006. Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake, *Science*, **313**(5787), 658–662.

Kanamori H. & McNally K.C., 1982. Variable rupture mode of the subduction zone along the Ecuador-Colombia coast. *Bull. seism. Soc. Am.*, **72**(4), 1241–1253.

Lefèvre, F., Lyard, F. & Le Provost, C., 2000. FES98: a new global finite element solution independent of altimetry, *Geophys. Res. Lett.*, **27**(17), 2717–2720.

Mikhailov, V., Tikhotsky, S., Diament, M., Panet, I. & Ballu, V., 2004. Can tectonic processes be recovered from new gravity satellite data? *Earth planet. Sci. Lett.*, **228**, 281–297.

Miura, S., Suwa, Y., Hasegawa, A. & Nishimura, T., 2004. The 2003 M8.0 Tokachi-Oki earthquake—how much has the great event paid back slip debts? *Geophys. Res. Lett.*, **31**, L05613, doi:10.1029/2003GL019021.

Ogawa R. & Heki K., 2007. Slow postseismic recovery of geoid depression formed by the 2004 Sumatra-Andaman Earthquake by mantle water diffusion, *Geophys. Res. Lett.*, **34**, L06313, doi:10.1029/2007GL029340.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **75**, 1135–1154.

Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space., *Bull. seism. Soc. Am.*, **82**, 1018–1040.

Panet, I. et al., 2007. Co-seismic and post-seismic signatures of the Sumatra December 2004 and March 2005 earthquakes in GRACE satellite gravity, *Geophys. J. Int.*, **171**, 177–190 doi:10.1111/j.1365-246X.2007.03525.x.

Plafker, G. & Savage, J.C., 1970. Mechanism of the Chilean Earthquakes of May 21 and 22, *Geol. soc. Am. Bull.*, **81** 1960, 1001–1030.

Preisendorfer, R.W., 1988. *Principal Component Analyses in Meteorology and Oceanography*, Elsevier, Amsterdam, the Netherlands.

Robinson, D.P. & Sandron D., 2006. A Mw 8.1 Earthquake on a fossil fracture zone: the December 23, 2004 Macquarie Ridge Earthquake, in *Proceedings of the AGU Fall Meeting*, San Francisco.

Savage, J.C. & Hastie, L.M., 1966. Surface deformation associated with dip-slip faulting, *J. geophys. Res.*, **71**, 4897–4904.

Schmidt, R. et al., 2006. GRACE observations of changes in continental water storage, *Global planet. Change*, **50**, 112–126.

Sun, W. & Okubo, S., 1993. Surface potential and gravity changes due to internal dislocations in spherical Earth, 1993, I: theory for a point dislocation, *Geophys. J. Int.*, **114**, 569–592.

Tanioka, Y. & Hasegawa, Y., 2007. Analysis of the tsunami generated by the 2006 great Kurile earthquake (Mw 8.2) EGU topic conference ‘The Role of Geophysics in Natural Disaster Prevention’, Lima, Peru, 05–09 March 2007, 2nd Alexander von Humboldt International Conference.

Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins, M.W., 2004. GRACE measurements of mass variability in the Earth system, *Science*, **305**, 503–505.

Wahr, J., Molenaar, M. & Bryan, F., 1998. The time variability of the Earths gravity field: hydrological and oceanic effects and their possible detection using GRACE, *J. geophys. Res.*, **103**, 30 205–30 229.

Wright, T., Lu, Z. & Wicks, C., 2004. Constraining the slip distribution and fault geometry of the Mw 7.9, 3 November 2002, denali fault earthquake with interferometric synthetic aperture radar and global positioning system data, *Bull. seism. Soc. Am.*, **94**(6B), S175–S189.