GRACE gravitational measurements of tsunamis after the 2004, 2010, and 2011 great earthquakes

The 2004 Sumatra, 2010 Maule, and 2011 Tohoku great earthquakes triggered tsunamis as large as a few decimeters over a few 100 km in the open ocean. The transient ocean mass redistribution propagating as tsunamis changed the Earth’s gravity field enough to perturb the GRACE satellites’ orbits at ~ 500 km above the surface. The on-board microwave ranging system detected inter-satellite acceleration anomalies of up to 1.0–4.0 nm/s2. There is good agreement between GRACE measurements and tsunami models for the three events. Complementarily to buoys, ocean bottom pressure sounders, and satellite altimeters, GRACE is sensitive to the long-wavelength spatial scale of tsunamis and provides an independent source of information for assessing alternate early earthquake and tsunami models. Our study demonstrates an innovative way of applying GRACE and GRACE Follow-On data to detect transient geophysical mass changes which cannot be observed by the conventional monthly Level-2 and mascon solutions.


Introduction
The Gravity Recovery and Climate Experiment (GRACE) mission measured time-variable global gravity fields with a spatial resolution of 300-400 km nearly every month for 15 years from 2002 to 2017 (Tapley et al. 2004). The conventional data products of the GRACE mission include monthly averaged global snapshots of Earth's gravity field (Level-2 data product, L2) and of surface mass change (Level-3 data products, L3) inferred from the L2 data after removing the elastic Earth loading effect and non-surficial processes such as glacial isostatic adjustment and earthquakes. The monthly surface mass redistribution data have been used to quantify variabilities in terrestrial water storage, ice-sheets and glaciers, and ocean mass (e.g., Tapley et al. 2019). Geophysical processes at such spatial scales are generally characterized by seasonal and inter-annual changes. The monthly solutions seem sufficient to recover the majority of the signals.
However, there are substantial gravitational signals not captured in the conventional monthly L2 data, but that can be obtained from a tailored processing of the GRACE data. For example, the sub-monthly barotropic gyre in the Argentine basin is observed in 10-day gravity solutions (Bruinsma et al. 2010;Han et al. 2014;Yu et al. 2018). Other higherfrequency ocean variability is visible in the residual Level-1B (L1B) inter-satellite tracking data (Chambers and Bonin 2012; see also Fig. 1). Various oceanic mass variabilities at different time scales are manifested in inter-satellite ranging data measured by the GRACE K-band ranging (KBR) Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0019 0-020-01395 -3) contains supplementary material, which is available to authorized users. system. In this study, we examined unusually intense but short-lasting gravitational variability caused by tsunami waves, discernible at satellite altitudes for less than a day after the earthquakes. Such transient geophysical signals are rarely addressed with GRACE data. Recently, Ghobadi-Far et al. (2019) studied the Earth's free oscillations excited by the 2004 Sumatra earthquake, which lasted for a few days after the earthquake rupture, using the GRACE KBR residuals.
Transoceanic tsunamis are mainly triggered by sudden displacement of the ocean floor by earthquakes. During the GRACE era from 2002 to 2017, three great earthquakes [2004 Sumatra (Indonesia), 2010 Maule (Chile), and 2011 Tohoku (Japan)] generated large transoceanic tsunamis. In the open ocean, changes in the ocean mass load caused by the tsunami waves were detected by ocean bottom pressure gauges (e.g., Saito et al. 2010;Yoshimoto et al. 2016) and the tsunami sea surface height variations were observed by satellite altimetry (e.g., Smith et al. 2005;Song et al. 2005;Satake 2007, 2013;Song et al. 2012;Hamlington et al. 2012). These observations have been used for studying the mechanism of large tsunamis in the open ocean, improving tsunami model simulation and prediction, as well as constraining earthquake sources. However, these measurements are local samples of extensive tsunamis and may not necessarily provide information on the overall sea surface disturbance.
In this study, we introduce a new space-borne technique with a wide field of view for offshore tsunami observation that measures the gravitational effect of tsunami-induced ocean mass change. The gravitational signature of tsunamis is manifested in the inter-satellite ranging between the two GRACE satellites. We demonstrate how inter-satellite distances were perturbed by transient gravity changes from the 2004 Sumatra, 2010 Maule, and 2011 Tohoku tsunamis, and present unequivocal evidence of the tsunami waves detected in the actual GRACE KBR measurements over 1-14 h after the earthquakes.

GRACE inter-satellite ranging measurements
We first computed dynamic reference orbits of the GRACE satellites based on the data processing strategy and background models used in ITSG-Grace2018 monthly gravity solutions (Mayer-Gürr et al. 2018) for December 2004, February 2010, and March 2011, covering the time periods of the Sumatra, Maule, and Tohoku tsunamis, respectively. The reference orbits were computed on a daily basis using the non-gravitational accelerometer observations and star camera orientation data as well as the ITSG-Grace2018 monthly mean solution along with several other background models such as the high-frequency (non-tidal) atmosphere and ocean mass variability AOD1B RL06 (Dobslaw et al. 2017), and the updated ocean tides model FES2014b (Carrère et al. 2015) using several years of the GRACE data (Mayer-Gürr et al. 2012). Furthermore, to isolate the gravity change associated with the tsunami from the coseismic gravity change, the coseismic earthquake geopotential models (Han et al. 2013) were also included in the background force models used in the orbit integration (see Figure S1 in supporting information).
The reference range-rate data computed from the dynamic orbits were then subtracted from the GRACE range-rate measurements to obtain the range-rate residuals. As the last step, we applied numerical differentiation to compute range-acceleration residuals from range-rate residuals. The choice of using range-acceleration was made simply because accelerations are better collocated with the location of mass redistribution on the Earth's surface, and thus, they are easier to interpret for regional gravity changes. However, the numerical differentiation amplifies random noise at high frequencies, and hence, low-pass filtering is necessary to accentuate the signals from GRACE range-acceleration data. In our convention, a positive (or negative) mass anomaly on the surface is shown as a trough (or crest) in range-acceleration (Ghobadi-Far et al. 2018b).
The range-acceleration residuals over the oceans represent the high-frequency (sub-monthly) ocean variability not modeled in the AOD1B model as well as the KBR and accelerometer measurement noise and background model errors. Note that GRACE monthly mean solutions included in the The largest anomalies are found mostly near the coastal regions. The residuals represent the sub-monthly gravity variation (due to mass redistribution) and are shown to provide context for the larger residuals observed in the open ocean due to tsunamis background models removed the majority of slowly varying (e.g., seasonal and inter-annual) oceanic signals. The root mean square (RMS) of the range-acceleration residuals is typically ~ 0.4 nm/s 2 in the deep ocean and as large as 1.0 nm/s 2 in some coastal regions (see Fig. 1).

Tsunami modeling
We simulated the tsunami wave propagation using the approach of Allgeyer and Cummins (2014) which solves the shallow-water (2D) equations while accounting for the elastic response of the Earth to tsunami mass load as well as seawater density stratification. Implementation of the elastic loading and density stratification considerably improves the predictions of the travel time and waveform of tsunami waves (see also Inazu and Saito 2013;Tsai et al. 2013;Watada 2013). Moreover, the initial sea level condition for tsunami wave propagation was computed from the water column change by both vertical and horizontal seafloor deformation. The importance of the horizontal forcing in tsunami generation is extensively discussed by Song and Han (2011). We computed the seafloor deformation based on the earthquake source models of Ammon et al. (2005) and Hébert et al. (2007) for the 2004 Sumatra, Fujii and Satake (2013) for the 2010 Maule, and Satake et al. (2013) for the 2011 Tohoku events. These source models for the three earthquakes were, respectively, obtained by geodetic and seismological inversions, tsunami waveforms plus coastal geodetic data, and tsunami waveforms. The tsunami models were computed on a 2′ × 2′ grid with a temporal resolution of 1 s for 24 h of propagation. The output has been downsampled to a 0.5° × 0.5° grid (average filter) with a temporal resolution of 1 min for the purpose of this study.

Method
We developed a method based on forward modeling of the gravitational effect of tsunami waves along the satellite orbit and direct comparison with the GRACE L1B observation of inter-satellite ranging perturbation. We carried out the model-data comparison in terms of line-of-sight gravity difference (LGD). LGD can be directly estimated from range-acceleration residual, however with a limited accuracy (Weigelt 2017;Ghobadi-Far et al. 2018a). Recently, Ghobadi-Far et al. (2018b) developed a transfer function based on correlation-admittance spectral analysis of rangeacceleration and LGD that can be applied to residual rangeacceleration, enabling accurate LGD determination with an error of 0.15 nm/s 2 over the frequency band higher than 5 cycles per revolution (CPR; 1 CPR = 0.18 mHz) or 1 mHz (wavelengths of 8000 km and smaller). The method is most appropriate for regional applications because of limited accuracy at the low frequency band (less than 5 CPR).
To compute the gravitational effect of the tsunami wave field along the GRACE orbits in terms of LGD, the grids of tsunami sea surface height change (with a temporal resolution of 1 min) were converted into spherical harmonic coefficients of gravitational potential to degree 200. Then, the gravity vector was computed for each GRACE satellite along its reference orbit using the tsunami geopotential coefficients. The tsunami-induced LGD changes were simply obtained by projecting the difference of the gravity vectors between the two satellites in the line-of-sight direction (Ghobadi-Far et al. 2018b).

2004 Sumatra tsunami
We first extracted the GRACE ground tracks in the region from longitudes of 35° E to 130° E and latitudes of 40° S to 26° N over the 24 h after the Sumatra earthquake, which occurred on 26 December 2004 00:58:53 UTC ( Figure S2). We then selected those arcs that matched the tsunami wave propagation in time and space ( Figure S3). The GRACE satellites first flew over the tsunami 1 h 15 min after the earthquake in an ascending track, when the tsunami wave height was as large as 40 cm and the induced surface gravity change was ~ 20 µGal ( Figure S4).
We compared the GRACE observed LGD with the tsunami model LGD computed using the earthquake source model from Ammon et al. (2005), for the arc shown in Fig. 2a. Figure 2b shows the model LGD (in blue) and GRACE observation (in red) time series as a function of latitude. We applied low-pass filtering at 50 CPR to suppress the high-frequency noise in the GRACE LGD data. The two negative anomalies of − 2.5 and − 2 nm/s 2 found in the GRACE LGD data at latitudes 5° S and 15° N, respectively, correspond to the positive tsunami wave around the same latitudes where the southern and northern leading edges of the tsunami are located. The large crest of 4 nm/s 2 found at latitude 5° N is caused by the notable negative sea level anomalies found in the model in that region. The remarkable agreement between the GRACE data and tsunami model from latitude 10° S to 20° N expresses the exceptionally large tsunami (40 cm over wavelengths of few hundred km). The orbits of the GRACE satellites were significantly perturbed by the transient gravity change due to the tsunami so that their inter-satellite distance changes were detectable by the KBR system.
If the GRACE LGD observations along the ascending arc in Fig. 2b were to be associated with transient tsunami waves, they should have not appeared along the descending arc crossing the same region with ~ 12 h of time difference because tsunamis would not last over the region for such a long period. To demonstrate this, we compared the GRACE LGD data along the ascending track ( Fig. 2a) with the descending track (Fig. 2c) covering the same region most closely but roughly 12 h later. The LGD data from the descending track and the corresponding model synthetics are shown in Fig. 2d. With the diffusive tsunami propagation, we do not expect to observe the same size of gravitational perturbation 12 h later. As expected, both the GRACE data and the model show no significant anomaly in the descending track (Fig. 2d). For other cases of GRACE overpassing tsunamis (such as arc numbers 5, 6, and 7 in Figure S3), the tsunami signal did not emerge from the GRACE data due to the measurement noise as well as natural oceanic variability.
To quantify the agreement between GRACE and the tsunami model, we applied consistent low-pass filters to both the GRACE and tsunami LGD time series at 40, 50, and 60 CPR frequencies (see Figure S5). The correlation coefficients between GRACE and the tsunami model from latitude 10° S to 20° N, where GRACE satellites were over the tsunami waves, are 0.99, 0.98, and 0.98. The RMS signal reduction, which is computed as 1 − RMS(LGD GRACE −LGD tsunami ) RMS(LGD GRACE ) , is 0.79, 0.76, and 0.74, respectively. We note that the same results in terms of correlation and RMS reduction were obtained when we applied a band-pass filter between 5 and 50 CPR to GRACE and tsunami data along the arc number 1. To demonstrate the statistical significance of the results, we also compared the GRACE LGD along the arc number 5 sampling the same region with ~ 12 h of difference ( Fig. 2c) with the tsunami LGD along arc number 1. The c The GRACE descending track covering the same region ~ 12 h after the ascending track in a along with the tsunami wave field at the same time. d Comparison between the GRACE observed and tsunami simulated LGD along the descending track shown in c. The comparison of ascending and descending arcs demonstrates that the ocean mass change detected by GRACE along the ascending arc #1, which matches the tsunami model, is indeed associated with the transient tsunami waves, because the same GRACE signal does not appear along the descending arc #5 passing over the same region. The star indicates the epicenter of the earthquake, and LGD denotes line-of-sight gravity difference. For the reference to the GRACE arc numbers, see Figure S2 in supporting information. GRACE satellites completed one orbital revolution around the Earth in ~ 90 min, so it took them ~ 15 min to cross the Indian Ocean from latitude 40° S to 20° N correlation and RMS reduction in this case is reduced to − 0.50 and − 1.84, respectively (see Table S1 in supporting information). This is due to the dissipation of tsunami waves over time.
We examined the GRACE tsunami data for contrasting different earthquake source models, as done in Han et al. (2010) but using the solid Earth deformation signals (not tsunamis). We computed the model LGD perturbation with the seismic source of Hébert et al. (2007), shown in black in Fig. 2b. It is seen that the tsunami model with the seismic source from Ammon et al. (2005), which is closer to the moment release estimate by Okal and Stein (2009) using ultra-long period normal mode data, compares better with the GRACE data in terms of both amplitude and phase. This example shows a potential benefit of GRACE-like data for assessing alternate early source models, improving tsunami model simulations, and providing enhanced constraints on seismic source parameters.

2010 Maule tsunami
We found a total of 14 GRACE ground tracks in the region from longitudes of 130° E to 300° E and latitudes of 70° S to 70° N in the 24 h after the 2010 Maule earthquake, which occurred on 27 February 2010 06:34:11 UTC ( Figure S6). Unfortunately, the first two tracks that match the tsunami propagation in time and space do not occur until 11, and 12 h 30 min after the earthquake (arc numbers 7 and 8 in Figure  S7). By that time, the tsunami wave energy had dissipated, with the largest remaining amplitude being < 10 cm at the wave front.
We examined in detail arc number 9, where the GRACE satellites passed over the center of the leading edge of the tsunami waves in an ascending track ~ 14 h after the earthquake. Figure 3a shows that the leading edge at latitude of ~ 30° E propagates northwest with an amplitude of ~ 10 cm and an equivalent surface gravity change of ~ 4 µGal ( Figure S8). The observed and simulated LGD agrees over latitudes 20° N to 40° N (Fig. 3b). GRACE detected  Figure S10 in supporting information for another similar case the positive sea level change at the leading edge of the tsunami as well as other oceanic variability in other regions with similar sea level perturbation. Such a signal is not found in the GRACE LGD from latitude 20° N to 40° N along the descending track crossing the same region ~ 12 h before (Fig. 3c, d). The correlation and RMS reduction between the low-pass filtered GRACE and tsunami model LGD at 40 CPR along the ascending arc number 9 from latitude 20° N to 40° N are ~ 0.9 and 0.60, respectively, and slightly less agreement is found when higher CPR is used (see Figure S9). To show the statistical significance of these values, we note that the correlation and RMS reduction between GRACE LGD along the descending arc number 2 (Fig. 3c), which passes over the same region as the arc number 9, and the transient tsunami LGD along the arc number 9 are − 0.30 and − 0.53, respectively (see Table S1).
There is another arc of the GRACE satellites crossing over the tsunami (arc number 10 in Figure S7), an ascending arc 90 min after arc number 9 shown in Fig. 3. The results of GRACE and model comparison for this pass are shown in Figure S10. The largest signals are observed near the equator, where the GRACE satellites detected the wave front.

2011 Tohoku tsunami
The 2011 Tohoku earthquake triggered a catastrophic tsunami on 11 March 2011 at 05:46:24 UTC. The GRACE satellites first flew over the tsunami on a descending track about 3 h 45 min after the earthquake ( Figure S11 and S12). Figure 4a shows that the GRACE satellites flew over the tsunami wave field into the northward propagating wave front at latitude of 55° N, and exited from the southward propagating wave front at latitude 10° N. The southward leading edge is particularly prominent with an amplitude of 20 cm and surface gravity change of 8 µGal ( Figure S13). The GRACE LGD observations agree well with the tsunami model over and S16 in supporting information for other possible cases of tsunami detection by GRACE latitudes of 5° N to 60° N where the model predicts a large LGD perturbation (Fig. 4b). The sharp and intense peak of the LGD anomaly of up to − 2 nm/s 2 at 10° N is a clear indication of the southern wave front.
The GRACE arc number 9 sampled the same area in an ascending track 12 h after arc number 3 (Fig. 4c). By that time, the tsunami wave field had diminished in that area to less than ± 5 cm, resulting in a less distinguishable tsunami LGD anomaly relative to inherent ocean variability (Fig. 4d). Low-pass filtering of the observed and simulated LGD at the same cutoff frequencies along the descending arc number 3 from latitude 0° N to 60° N results in correlation ≥ 0.8 and RMS reduction ≥ 0.40 (see Figure S14). These high correlation and RMS reduction values are reduced to − 0.16 and − 0.67, respectively, when GRACE LGD along the ascending arc number 9 (Fig. 4c), which covers the same region as arc number 3, is compared to transient tsunami LGD along the arc number 3, illustrating that the gravitational anomalies from latitude 0° N to 60° N were indeed transient (less than ~ 12 h) as expected from tsunami-induced ocean mass changes (see Table S1). The additional comparison for two other possible arcs (number 7 and 8) is also shown in Figures S15 and S16, where the signals from the GRACE data and the model are smaller but still positively correlated.

Summary and prospect
We presented a novel approach of analyzing GRACE intersatellite ranging measurements to model transient gravitational changes associated with the 2004 Sumatra, 2010 Maule, and 2011 Tohoku tsunamis. The method is based on the residual KBR data with respect to the monthly mean gravity field. The transient signals are not modeled in the monthly mean L2 gravity fields but remain in the residual L1B data. We exploited the linear relationship in the spectral domain between inter-satellite range-acceleration (second time-derivative of range) and line-of-sight gravity difference between the two GRACE satellites to evaluate the gravitational perturbation by tsunamis along the satellite orbits.
All three tsunamis were as large as a few decimeters over wavelengths of hundreds of km in the open ocean. Expectedly, the largest anomalies of sea surface height changes were observed at the wave front. Despite the polar orbit of GRACE satellites, there were still a few instances when the satellites were flying over the tsunami waves before they dissipated after 10-15 h. The KBR residuals in the open ocean are generally on the order of 0.4 nm/s 2 in terms of RMS variability, roughly translated into several cm of ocean mass variability. The tsunami gravitational changes measured by GRACE were up to 4.0, 1.0, and 1.5 nm/s 2 , respectively, for the 2004, 2010, and 2011 tsunamis. There are high correlations between GRACE measurements and tsunami models (up to 0.98) for all events and all cases of GRACE overpasses over the tsunami wave field, when greater than ~ 10 cm in amplitude.
It should be mentioned that the discrepancy between GRACE observations and tsunami models is mainly due to modeling uncertainty (tsunami and earthquake) as well as mis-and un-modeled high-frequency atmosphere and ocean signals in the GRACE KBR residuals. Comparison of GRACE observations along the ascending and descending arcs covering the same region (but ~ 12 h apart) showed the high-frequency nature of such observations and provided further evidence for tsunami detection. We provided further analysis to demonstrate the statistical significance of our results by comparing tsunami signals with all GRACE data along the arcs covering the same region in December 2004, February 2010, and March 2011 (see Table S2).
Whereas GPS can measure coseismic land deformation, and DART sensors or satellite altimeters observe tsunami sea surface height changes, GRACE data are sensitive to gravitational changes caused by both tsunami and land deformation. In contrast to point-wise ocean bottom pressure measurements and cross-sectional satellite altimeter profiling of tsunami heights, the GRACE gravimetric observations provide regional averages of tsunami-induced gravity changes in the open ocean. The GRACE data essentially represent the long-wavelength components of tsunami (typically ≥ 2000 km; see Fig. 5); in the shallowwater approximation, this corresponds to frequencies of 10 −4 Hz and under, which transcend the domain accessible by seismological studies, intrinsically limited to the gravest mode of the Earth (0.31 mHz). This is a result of the layer of space between the Earth's surface and the satellites' orbit acting as a powerful low-pass spatial filter for perturbations of the gravitational potential, which must satisfy Laplace's equation in a medium with zero density, and thus decay rapidly with height for smaller lateral wavelengths. In this context, the GRACE methodology introduced here has the potential to explore the properties of tsunamis (and hence of their parent earthquakes) in a hitherto unexplored domain of frequencies. Being unusually sensitive to the regional scale of tsunamis, the GRACE data can be used for assessing early models of earthquakes and tsunamis after the events to provide independent estimates of the overall tsunami size which are possibly less biased than those derived from sparsely located in situ buoy records.
Our study demonstrated a new application of GRACE microwave inter-satellite tracking data to study gravitational changes from transient (high-frequency) geophysical processes with tsunamis as an example. Equipped with an improved KBR system and laser ranging interferometer (Landerer et al. 2020), the GRACE Follow-on satellites will reduce the detection threshold for measuring gravity changes and may lead to more accurate detection of tsunamis as they propagate across major ocean basins, as well as other transient geophysical mass changes. The L1B processing strategy can be optimized so as to deliver low-latency GRACE Follow-on measurements. LGD in nm/s 2 and tsunami height in cm, respectively. Note that the sign of LGD is reversed for better visual comparison. This figure illustrates that the GRACE gravimetric observable, as opposed to satellite altimetry measurements, is especially sensitive to the long-wavelength components of the tsunami wave field