# Rapid identification of tsunamigenic earthquakes using GNSS ionospheric sounding

Fabio Manta, Giovanni Occhipinti, Lujia Feng, Emma M. Hill

► To cite this version:

Fabio Manta, Giovanni Occhipinti, Lujia Feng, Emma M. Hill. Rapid identification of tsunamigenic earthquakes using GNSS ionospheric sounding. Scientific Reports, 2020, 10, p. 469-489. 10.1038/s41598-020-68097-w. insu-03584751

HAL Id: insu-03584751

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

Submitted on 24 Feb 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
<figure><img src="image_2.png" /></figure>

OPEN

# Rapid identification of tsunamigenic earthquakes using GNSS ionospheric sounding

Fabio Manta¹,²,⁴✉, Giovanni Occhipinti¹,³,⁴, Lujia Feng¹,⁰ & Emma M. Hill¹,⁰,¹

The largest tsunamis are generated by seafloor uplift resulting from rupture of offshore subduction-zone megathrusts. The rupture of the shallowest part of a megathrust often produces unexpected outsize tsunami relative to their seismic magnitude. These are so called 'tsunami earthquakes', which are difficult to identify rapidly using the current tsunami warning systems, even though, they produce some of the deadliest tsunami. We here introduce a new method to evaluate the tsunami risk by measuring ionospheric total electron content (TEC). We examine two $M_w$ 7.8 earthquakes (one is a tsunami earthquake and the other is not) generated in 2010 by the Sunda megathrust, offshore Sumatra, to demonstrate for the first time that observations of ionospheric sounding from Global Navigation Satellite System (GNSS) can be used to evaluate the tsunamigenic potential of earthquakes as early as 8 min after the mainshock.

'Tsunami earthquakes', as originally defined by Kanamori¹, are events generating tsunami with larger amplitude than expected from their seismic magnitude. Most tsunami earthquakes are generated by high levels of slip on the shallow megathrust, which results in large seafloor uplifts and hence very dangerous tsunami. The shallow location of the slip—close to the subduction trench—means that the ruptures generating tsunami earthquakes are at significant distance from land-based monitoring networks, limiting our ability to quickly and accurately assess their magnitude and source parameters. Conventional approaches using various seismological methods²⁻⁴ or rapid inversion of GNSS (Global Navigation Satellite System) estimates of ground motion⁵ regularly encounter difficulties in accurately estimating the uplift of the seafloor and consequently fail in predicting the tsunamigenic nature of tsunami earthquakes. These events are generated in the relatively soft accretionary wedge and are characterized by slow rupture velocities, and therefore lower levels of seismic energy radiation⁶,⁷; they are thus only weakly felt by nearby populations, prohibiting natural warnings.

In addition to seismic waves propagating within the solid Earth, tsunamigenic earthquakes can induce waves in the atmosphere, including the ionosphere. Firstly, the uplift of the seafloor at the source generates an acoustic-gravity pulse ($AGW_{epi}$) that propagates through the atmosphere above the epicentral area (Fig. 1). This perturbation travels at the speed of high-atmospheric acoustic waves (~1,000 m/s) and reaches the ionosphere approximately 8 min after the event (Fig. 1)⁸. Due to the short response time, this ionospheric perturbation can provide critical information for rapidly assessing the tsunami potential. Secondly, if an ocean tsunami is generated, another perturbation later reaches the ionosphere. A tsunami is essentially an oceanic gravity wave; consequently, by linear dynamic coupling, it produces an atmospheric internal gravity wave ($IGW_{tsuna}$), which reaches the ionosphere approximately 40 min after the event (Fig. 1)⁹. During the upward propagation, the generated $AGW_{epi}$ and $IGW_{tsuna}$ are both strongly amplified by the double effects of the exponential decrease of the atmospheric density $\rho$ with altitude and the conservation of the kinetic energy $\rho v^2$, where v is the local speed oscillation induced by wave crossing. At ~300 km altitude in the ionosphere, the $AGW_{epi}$ and $IGW_{tsuna}$ are amplified enough to induce strong perturbations in the plasma density. These perturbations can be observed by ionospheric monitoring techniques, such as measurements of Total Electron Content (TEC) by GNSS (see¹⁰ for a review) or airglow cameras¹¹,¹². The TEC represents the number of electrons integrated along the ray-path between the GPS satellite and the station (see “Material and methods” for more details about TEC calculation). The signature of tsunamigenic earthquakes in the GNSS-derived TEC is consequently detectable in both the

¹Earth Observatory of Singapore, Nanyang Technological University, Singapore, Singapore. ²Asian School of the Environment, Nanyang Technological University, Singapore, Singapore. ³Institut Universitaire de France, Paris, France. ⁴Present address: Université de Paris, Institut de Physique du Globe de Paris, CNRS, 75005 Paris, France. ✉email: manta@ipgp.fr
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Schematic of the two ionospheric signals associated with a tsunamigenic earthquake. At the time of the earthquake ($t_0$), the sudden rupture of the fault causes an uplift of the sea floor and subsequently triggers a tsunami. Acoustic Gravity Waves (AGW$_{epi}$) associated with the uplift propagate through the atmosphere and reach the ionosphere ~8 min after the rupture. The tsunami is an oceanic internal gravity wave; consequently, by a linear coupling, it produces an atmospheric Internal Gravity Wave (IGW$_{tsuna}$). The IGW$_{tsuna}$ reaches the ionosphere at time $t_1$, ~40 min after the tsunami is initiated. (b) Main geological features of the Sumatran subduction zone. The contoured areas represent estimated coseismic rupture patches for the 6 April 2010 $M_w$ 7.8 Banyaks earthquake in the north¹⁹, and the 25 October 2010 $M_w$ 7.8 Mentawai earthquake in the south²⁰. The yellow dots represent the SuGAr stations existing and available at the time of the events, while the green dots represent the current total extension of the SuGAr network.</figcaption></figure>

near-field by observing the AGW<sub>epi</sub> directly related to the uplift of the seafloor<sup>13,14</sup> and the far-field by observing the IGW<sub>tsuna</sub> related to tsunami waves<sup>9,11,15-17</sup>.

Based on the relation between the AGW<sub>epi</sub>, which reaches the ionosphere approximately 8 min after an event, and the uplift of the seafloor, we investigate the potential of using GNSS observations of ionospheric TEC to discriminate tsunami earthquakes. We propose a technique to validate the risk estimated by conventional seismological methods, based on quantification of the intensity of the AGW<sub>epi</sub>, to discriminate the tsunamigenic nature of an earthquake. As a case study, we examine two $M_w$ 7.8 megathrust earthquakes that occurred in 2010 along the Sumatran subduction zone, which were both recorded by the Sumatran GPS Array (SuGAr) (Fig. 2). Although both events had the same magnitude, and were issued similar tsunami early warnings (later canceled), the 6 April 2010 Banyaks earthquake<sup>18,19</sup>, did not generate a significant tsunami and caused only minor earthquake-related damage to local infrastructure. On the contrary, the 25 October 2010 Mentawai event produced a large tsunami—recorded within 10 min of the earthquake by the Mentawai GPS buoy—which claimed more than 400 lives from communities along the coastlines of the Pagai Islands<sup>7,20,21</sup> (see “Material and methods” for a comprehensive description of the two events and the corresponding responses of tsunami warning systems).

Our analysis of the ionospheric response to these two events provides a unique natural experiment to test if the ionosphere is sensitive to the rupture properties of earthquakes and tsunami genesis, and whether GNSS-TEC observations are capable of identifying the events that have a high level of tsunami potential.

At the time of the Mentawai earthquake, 29 stations of the SuGAr network were in operation (Figs. 1b and 2a), while 25 stations were operating at the time of the Banyaks earthquake (Figs. 1b and 2b). In both cases, more than eight GPS satellites were visible in the sky, and we selected the four that better detected the events with comparable observation geometries. Cahyadi and Heki<sup>22</sup> observed that changes of the penetration angles of the line-of-sight (LOS) with respect to the AGW<sub>epi</sub> wave front propagating in the ionosphere can influence the amplitude of the detected signals. The authors referred to “good LOS” as the most favorable scenario for which the ionospheric piercing point (IPP) is located between the epicenter and the GNSS station (Fig. 2), in this configuration the LOS crosses mainly the positive part of the wave front returning a strong TEC signal. On the contrary they referred to “poor LOS” as the scenario for which the GNSS station is located between the epicenter and the IPP. With this second geometric configuration the LOS simultaneously penetrates the positive and negative parts of the anomaly, partially cancelling each other and consequently returning a weak signature in the GNSS-derived TEC. Moreover, as a consequence of the integrated nature of TEC along the station-satellite LOS, satellites at low elevation angles have a better detection capability<sup>8</sup> (see Supplemental Figs. 1 and 2).

Accordingly, to be able to compare the TEC signatures of the two events analyzed in this work, we selected (from Fig. 2) four satellite-station pairs (Fig. 3) with similar angles and distances from the epicenter area that better represent the different LOS scenarios: Good LOS, with receiver-satellite pairs BTHL-PRN29 and UMLH-PRN17 for the Mentawai and Banyak Islands respectively; and Poor LOS, with BSAT-PRN14 and BSIM-PRN2 for the Mentawai and Banyak Islands respectively. During the Mentawai event the TEC shows a perturbation related to the AGW<sub>epi</sub> (Fig. 2) with a frequency signature between 1 and 7 mHz (Fig. 3a,b) appearing 8 min after the main shock (21:42:23 local time, UTC +7) and propagating away from the epicenter at a horizontal speed
<figure><img src="image_4.png" /></figure>

**Figure 2.** (a) Map of the distribution of ionospheric piercing points (IPPs) for the pairs of 29 GNSS stations and the four closest satellites at the time of the Mentawai earthquake (local time: 21:42:23 UTC +7) (b) Map of the distribution of ionospheric piercing points (IPPs) for the pairs of 29 GNSS stations and the four closest satellites at the time of the Banyaks earthquake (local time: 5:15:02 UTC +7). The IPPs are represented by four different symbols corresponding to different satellites; the symbols are located at positions along the trace of the IPPs corresponding to the event time. Note the bold traces highlighting the trajectory of the satellite-receivers pairs discussed in Fig. 3. The red stars represent the epicenter of the earthquakes. The cartoon of “Good LOS” (green) and “Poor LOS” (red) configurations is shown between a and b. (c) Hodochrones of the TEC perturbation observed by GPS satellites PRN 21, 9, 14, and 29. (d) Inset with color-scale 5 times smaller to highlight the gravity waves related to the IGWtsuna recorded by PRN 21. (e) Hodochrones of the TEC perturbation observed by GPS satellites PRN 2, 28, 5, and 17. Color-scale in (e) is 5 time smaller than (c). Vertical dashed lines represent the time of the event. The gray dashed lines show the speed of the AGWepi (800 m/s) and IGWtsuna (250 m/s).

of 600–800 m/s (Fig. 2c). Later, a second weaker TEC perturbation at a lower frequency (~1.5 mHz), related to the IGWtsuna appears to travel at a horizontal speed of ~250 m/s (Fig. 2d), consistent with the tsunami speed of ~220 m/s observed by DART buoy 56,001, located ~1,600 km from the epicenter.⁶ The arrival time of IGWtsuna observed at the distance of 480 km from the epicenter (Fig. 3a) is coherent with the tsunami propagation. A much smaller TEC perturbation followed the Banyaks earthquake (Figs. 2e, 3c,d), which generated a much weaker tsunami that did not cause any damage (see “Material and methods” section).

Figure 3 shows that regardless the LOS scenario, the Mentawai earthquake induces a greater signature in the GNSS-TEC highlighting the important amount of energy transferred in the ionosphere by this event. On the contrary for the Banyaks earthquake the TEC signature remains very weak also with a Good LOS geometry. Indeed, the GNSS-TEC traces in the hodochrones in Fig. 2c,f show that the satellite coverage is dense enough to assure the favorable observation of the AGWepi and avoid attenuation effects related to the GNSS LOS geometry²².

Spectral analysis (Fig. 4) of observed TEC during the days before and after the Mentawai event reveals the unique characteristics of the AGWepi compared to the mean background level (MBL). The Ionospheric Tsunami Power Index (ITPI) is here introduced to define the ratio between TEC perturbation and MBL in order to remove the ionospheric dynamics, and consequently be able to compare the ITPI computed in different ionospheric
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Filtered ionospheric TEC time series and related spectrograms extracted by (a) observations of satellite PRN29 with respect to station BTHL (b) observations of satellite PRN14 with respect to station BSAT at the time of the 2010 *M*<sub>w</sub> 7.8 Mentawai earthquake (c) observations of satellite PRN17 with respect to station UMLH and (d) observation of satellite PRN2 with respect to station BSIM at the time of the 2010 *M*<sub>w</sub> 7.8 Banyaks earthquake. Grey curves represent the elevation angle of the satellite-receiver ray path. Solid red and white vertical lines indicate the time of the events. Dashed vertical lines indicate the time of the first potential arrival of the AGW<sub>epi</sub> (8 min) and the minimum time of IGW<sub>tsuna</sub> observation (40 min) in the ionosphere. Horizontal dashed lines are the Brünt-Vaïsalla frequency that represents the limit between gravity and acoustic domains.</figcaption></figure>

conditions (see "Material and methods" section). We highlight that the ITPI is related to the entire spectrum of energy of the AGW<sub>epi</sub> and consequently provides an estimate of the strength of a perturbation detected in the ionosphere. Compared to the dTEC amplitude estimation used by previous authors, our approach is not affected by frequency pass-band filtering or polynomial fit to obtain the perturbation of TEC. The ITPI for the Mentawai event was 14, as TEC observations revealed a perturbation 14% larger than the MBL (Fig. 4). A comparable energetic signature related to external localized phenomena (e.g., plasma bubbles, traveling ionospheric disturbances) appeared on 22 October 2010 but this was not related to a seismic event and consequent tsunami early warning, showing the importance to couple ionospheric observations with conventional techniques. Weaker TEC perturbations followed the Banyaks earthquake, at approximately 5% of the MBL (Supplemental Fig. 3), confirming the presence of a weaker tsunami (ITPI=5). The ITPI measured for the two events shows a correlation with the values of seafloor uplift reported by Cahyadi and Heki<sup>22</sup>, by the United States Geological Survey (USGS, https://earthquake.usgs.gov/earthquakes/), and by Hill et al.<sup>20</sup>. We complete our observations for the Mentawai and Banyak Islands with TEC observations for additional 17 tsunamigenic events extracted by Cahyadi and Heki<sup>22</sup> and listed in Table 1. We correlated the TEC amplitude (blue symbols in Fig. 5) with the seafloor uplift and, for comparison, with our new ITPI. To better quantify the contribution of the seafloor uplift, we introduce a weighting factor here referred as the seafloor maximum volume displaced (V<sub>max</sub>). This parameter is obtained from the product of the maximum seafloor uplift and the rupture area (A) (see "Material and methods" section). Seafloor uplift is in reality quite spatially variable, often poorly constrained by slip inversions, and likely to be supplemented by inelastic uplift that is not accounted for in most models. However, this parameter can provide a first-order proxy for comparison with the ITPI. To emphasize the originality of our work we highlight that Cahyadi and Heki<sup>22</sup> correlate the TEC (%) perturbations with the moment magnitude of the events and not with the consequent seafloor maximum volume displaced (V<sub>max</sub>), which takes into account also the extent of the seafloor uplift directly related to the tsunami genesis. Our empirical model, showed in Fig. 5, supports the possibility of obtaining information about the tsunami genesis from ionospheric perturbations.

Coupled with early seismic observations of an earthquake (estimating magnitude and area of slip), the detection of the AGW<sub>epi</sub> by GNSS-TEC could provide valuable additional information—at the optimistic scenario of 8 min following the initiation of earthquake rupture—on the level of seafloor uplift and thus the potential of a significant tsunami. We highlight that the amplitude of the ionospheric TEC perturbation induced by the AGW<sub>epi</sub> is directly related to the seafloor maximum volume displaced (V<sub>max</sub>) (Fig. 5) and we do not need to intermediatey evaluate the magnitude of the event or the source parameters to estimate the tsunami risk. Additionally, estimation of GNSS high-rate position has been demonstrated as a valid technique for source estimation and tsunami
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spectrograms of the filtered TEC signal recorded by the pair of station PTLO and satellite PRN29 at the epicenter of the 2010 Mw 7.8 Mentawai earthquake over seven consecutive days. Horizontal dashed lines are the Brünt-Vaïsalla frequency that represents the limit between the gravity and the acoustic domain. Vertical lines in the central panel (d) indicate, respectively from left to right: the time of the event (solid line); the first potential arrival of the AGW<sub>epi</sub> (8 min, dashed line); the minimum time of IGW<sub>tsuna</sub> observation (40 min, dashed line) in the ionosphere. Green contour lines mark where the intensity of the signal is above a threshold value that represents the mean background level (MBL) calculated on quiet days.</figcaption></figure>

risk evaluation to complement seismic networks (e.g.,⁵). Thus, we can use the same GNSS receivers for the source estimation and the consequent early tsunami potential evaluation using positioning techniques and validate it using the ITPI analysis. We can confirm the uplift and related tsunami risk by looking for the TEC perturbation induced by the AGW<sub>epi</sub> within 8 min after the mainshock. After 40 min the TEC perturbation induced by the IGW<sub>tsuna</sub> gives a direct estimate of the oceanic displacement induced by the tsunami¹⁷ propagating in the open ocean. The IGW<sub>tsuna</sub> signature in the TEC has been already observed in real-time²³, proving the potential role of ionospheric monitoring in real-time early warning.

With this work, we prove the maturity of ionospheric observations for evaluating tsunami potential, further showing the capability to discriminate the tsunamigenic nature of slow rupture events characterized by a low level of seismic energy radiation. In particular, estimations of ITPI could be fed into existing tsunami warning systems as an additional piece of information to alert for the presence of anomalous tsunami earthquakes and/or guard against false warnings. We hope to open a new paradigm in tsunami warning systems, based on the redundancy of information from the synergy of classic and new techniques to reduce damages and loss of life.
<table><caption>**Table 1.** List of the events discussed in Fig. 5. Values of Max water height, Damage, and Loss are from NOAA (https://ntwc.ncep.noaa.gov/). Mw moment magnitude, A rupture area, V<sub>max</sub> seafloor maximum volume displaced. <sup>a</sup>Obtained from (USGS). <sup>b</sup>Values from (Cahyadi and Heki<sup>22</sup>). <sup>c</sup>Value from (Hill et al.<sup>20</sup>).</caption><thead><tr><th>Date</th><th>Location</th><th>Mw</th><th>Max uplift (m)</th><th>A (km²)</th><th>V<sub>max</sub> (km³)</th><th>Max water height (m)</th><th>Damage</th><th>Loss</th></tr></thead><tbody><tr><td>04/10/94</td><td>Hokkaido-toho-oki</td><td>8.3</td><td>2<sup>b</sup></td><td>5,400</td><td>1.08E+4</td><td>10.40</td><td>Yes ~ $1M</td><td>-</td></tr><tr><td>25/09/03</td><td>Tokachi-oki</td><td>8</td><td>0.88<sup>a</sup></td><td>78,720</td><td>6.93E+4</td><td>4.40</td><td>Yes ~ $1M</td><td>-</td></tr><tr><td>05/09/04</td><td>Kii-hanto-oki(fore)</td><td>7.1</td><td>0.3<sup>b</sup></td><td>2,400</td><td>7.20E+2</td><td>0.63</td><td>No</td><td>-</td></tr><tr><td>05/09/04</td><td>Kii-hanto-oki</td><td>7.4</td><td>0.7<sup>b</sup></td><td>3,200</td><td>2.24E+3</td><td>0.93</td><td>No</td><td>-</td></tr><tr><td>26/12/04</td><td>Banda Aceh</td><td>9.2</td><td>3.4<sup>a</sup></td><td>240,000</td><td>1.2E+6</td><td>50.90</td><td>Yes ~ $10B</td><td>227,899</td></tr><tr><td>03/05/06</td><td>Tonga</td><td>7.9</td><td>0.53<sup>a</sup></td><td>10,075</td><td>2.41E+3</td><td>0.27</td><td>No</td><td>-</td></tr><tr><td>15/11/06</td><td>Central Kuril</td><td>8.2</td><td>0.46<sup>a</sup></td><td>32,000</td><td>7.38E+3</td><td>21.90</td><td>Yes ~ 1-$5 M</td><td>-</td></tr><tr><td>16/07/07</td><td>Chuetsu-oki</td><td>6.6</td><td>0.3<sup>b</sup></td><td>1,213</td><td>3.64E+2</td><td>Few cm</td><td>No</td><td>-</td></tr><tr><td>12/09/07</td><td>Sumatra Bengkulu</td><td>8.5</td><td>1.08<sup>a</sup></td><td>70,000</td><td>7.56E+4</td><td>5</td><td>Yes ~ $1M</td><td>-</td></tr><tr><td>12/09/07</td><td>Sumatra Bengkulu(after)</td><td>7.9</td><td>0.75<sup>a</sup></td><td>21,120</td><td>1.58E+4</td><td>-</td><td>No</td><td>-</td></tr><tr><td>15/07/09</td><td>New Zeland</td><td>7.8</td><td>1.15<sup>a</sup></td><td>3,200</td><td>3.68E+3</td><td>0.47</td><td>Yes</td><td>-</td></tr><tr><td>27/02/10</td><td>Maule</td><td>8.8</td><td>2.9<sup>a</sup></td><td>60,000</td><td>1.74 E+5</td><td>29</td><td>Yes ~ $30B</td><td>153</td></tr><tr><td>06/04/10</td><td>Banyaks</td><td>7.8</td><td>0.42<sup>a</sup></td><td>20,736</td><td>8.71E+3</td><td>0.44</td><td>No</td><td>-</td></tr><tr><td>25/10/10</td><td>Mentawai</td><td>7.8</td><td>1.9<sup>c</sup></td><td>31,500</td><td>5.39E+4</td><td>16.9</td><td>Yes ~ $39M</td><td>431</td></tr><tr><td>09/03/11</td><td>Tohoku-oki(fore)</td><td>7.3</td><td>0.17<sup>a</sup></td><td>2,457</td><td>4.18+2</td><td>0.60</td><td>No</td><td>-</td></tr><tr><td>11/03/11</td><td>Tohoku-oki</td><td>9</td><td>5.27<sup>a</sup></td><td>107,100</td><td>5.64E+5</td><td>38.90</td><td>Yes ~ $220B</td><td>18,434</td></tr><tr><td>28/10/12</td><td>Queen Charlotte Islands</td><td>7.8</td><td>0.32<sup>a</sup></td><td>4,200</td><td>1.3E+3</td><td>12.98</td><td>No</td><td>-</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relation between seafloor maximum volume displaced (V<sub>max</sub>), coseismic TEC amplitude and ITPI for 17 events. Blue symbols represent earthquakes from the literature<sup>22,24,25</sup>. Red symbols represent the Mentawai and Banyak earthquakes discussed in this paper for which both the coseismic TEC amplitude and the ITPI have been computed. The left axis represents the TEC coseismic amplitudes normalized by background vertical TEC estimated from Global Ionospheric Maps (GIM) as defined by Cahyadi and Heki<sup>22</sup>. The right axis represents the ITPI for the Mentawai and Banyak events scaled proportionally to the computed coseismic TEC amplitude. X symbols represent events that generated a very small tsunami that did not cause damage. O symbols represent events that did generate a tsunami and damage.</figcaption></figure>
## Material and methods

To investigate the TEC variations caused by the 2010 Mw 7.8 Mentawai and Banyaks earthquakes, we analyzed 15-s data from 29 and 25, respectively, continuous GNSS stations of the SuGAR network that were in operation on the day of each earthquake (Fig. 2). We analyzed the GNSS data for the week centered on the day of each earthquake to obtain the ionospheric perturbations before, during, and after the events.

We calculated the TEC by applying a method similar to the one described by Calais and Minster²⁶. To obtain a more accurate measure of the apparent distance between a satellite and receiver, we used the carrier phase (L1 & L2), neglecting the less precise pseudo-range measurements (P1 & P2):

$$TEC = \frac{(L1 - L2)}{40.3} \frac{f_1^2 f_2^2}{f_1^2 - f_2^2} \qquad (1)$$

where, $f_1$ and $f_2$ are, respectively, the corresponding high and low GPS frequency and TEC is the total electron content measured in TEC units (1 TECU = 10¹⁶ el/m²). We detected and corrected the cycle slips using the program *Ninja* within the GPS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) version 6.2²⁷.

As a consequence of the integrated nature of TEC, the observed ionospheric perturbations result from a large range of altitudes from the satellite to the receiver. However, the main contribution to the observed TEC variations is located around the height of the maximum of ionosphere ionization (the F2 layer). For this study, we fixed the altitude of the F2 layer at 300 km, which can be considered as the average altitude for day and night for the Equator. The intersection of the line of sight between each GNSS satellite-receiver pair at this altitude is termed the ionospheric pierce point (IPP), which we use to spatially visualize the observed TEC perturbations. In order to identify perturbations caused by the earthquake and the subsequent tsunami, we filtered all the initial TEC time series on a defined range of frequencies following the same methodology as Rolland et al.²⁸. We applied a 1 to 10 mHz bandpass finite impulse response (FIR) Butterworth filter in order to remove the contributions from daily ionospheric variabilities, satellite motions, and instrumental biases.

We introduce for the first time the Ionospheric Tsunami Power Index (ITPI) to rapidly identify tsunamigenic earthquakes:

$$ITPI = \left( \frac{\text{PSD}_{\text{MAX}}}{\text{MBL}} - 1 \right) \cdot 100. \qquad (2)$$

The ITPI corresponds to the ratio between the maximum level of power spectral density of the TEC (PSDMAX) recorded during the event, and the mean background level (MBL; defined as the maximum level of PSD recorded during the 6 days preceding the earthquake, averaged by the number of days). By averaging the background noise for several days prior to the event, we are able to remove outliers and smooth the noisy parts related to other perturbations traveling in the ionosphere. This new parameter could allow us to automatically identify the seafloor uplift and the consequent tsunamigenic potential of an event (Fig. 5). The relation for the Mentawai and Banyak events can be observed in Fig. 5, which shows, on the right side, the ITPI scale for the two events. The same relation with the seafloor maximum volume displaced $V_{max}$ can be observed for the Coseismic Ionospheric Disturbance (CID), which is shown on the left side for both the Mentawai and Banyak events and 17 earthquakes described in the literature²² and listed in Table 1. To calculate the seafloor maximum volume displaced ($V_{max}$), we multiply the maximum uplift value for the rupture area (A). Values of uplift are provided by the US Geological Survey (USGS, https://earthquake.usgs.gov/earthquakes/) resulting from finite fault, calculated using Okada-style deformation codes²⁹. We highlight that the ITPI is a more accurate way to estimate the ionospheric perturbation as it is calculated in the spectral domain and it is consequently free of eventual pass-band filter applied to CID.

**Tsunami warning response to the two Mw 7.8 Mentawai and Banyaks earthquakes.** We analyzed two Mw 7.8 earthquakes (Fig. 2) that occurred in 2010 along the Sunda megathrust (Indonesia): The 6 April Banyaks earthquake¹⁸,¹⁹ and the 25 October Mentawai earthquake²⁰,²¹. Although they had the same magnitude, only the Mentawai earthquake produced an anomalously large tsunami.

The Banyaks earthquake ruptured a deeper portion (20–30 km) of the megathrust, producing relatively small uplift of the seafloor and thus only a small tsunami (max water height: 44 cm, max run up: 6 m, source NOAA database). On the other hand, the Mentawai earthquake was a shallow event with estimated slip concentrated at depths of <6 km, no more than ~50 km away from the trench. Given the lower rigidity of the shallow sediments at this location, the Mentawai event generated considerably higher levels of slip²⁰,³⁰ than the Banyaks event¹⁸ while maintaining the same moment magnitude. The reported shake intensity at 150 km distance from the epicenter was MMI 5 and MMI 6 for the Mentawai, and for the Banyaks events, respectively. The high slip at shallow depths resulted in seafloor uplifts of up to several meters, producing the outsize tsunami with a reported maximum run up > 16 m²⁰.

For the Mentawai event, the Badan Meteorologi Klimatologi Dan Geofisika (BMKG) of Indonesia initially estimated the magnitude as 7.2 and issued a tsunami warning within 5 min of the earthquake. The warning was afterward cleared without receiving information of tsunami damage²¹. The German Indonesian Tsunami Early Warning System (GITEWS) recorded the tsunami on a surface GPS buoy located in the Mentawai islands (GITEWS SUMATRA-03) within 10 min of the earthquake, with amplitudes of ~15 cm and a period of a few minutes²¹,³⁰. Later the tsunami arrival was recorded by the Padang tide gauge about 1 h after the earthquake, and by the DART station 56,001 located around 1,600 km southeast of the epicenter about 2 h after the earthquake²¹,³⁰. However, no tsunami warning was issued because some buoys had been vandalized³¹ and were not in operation. For the Banyaks event, BMKG also recognized a tsunami potential and immediately issued a tsunami warning
that was canceled later, as only a minor tsunami was generated. For both events, the Pacific Tsunami Warning Center (PTWC) did not issue a warning, but a less severe tsunami watch, which was also canceled later (~2 h after). The GEOFON global seismological broadband network operated by the German GeoForschungsZentrum (GFZ) produced a real-time model for the two earthquakes, estimating a magnitude of $M_w$ 7.6 for the Banyaks event and a magnitude of $M_w$ 7.8 for the Mentawai event.

The Mentawai event did not provide a natural warning: witnesses reported feeling only weak, long-period shaking²⁰. While many of the people on the islands knew to evacuate to higher ground on feeling strong earthquake shaking, they did not recognize this event was so dangerous. The characteristics of tsunami earthquakes make them particularly dangerous and exemplify the need for additional warning systems.

**The current state of the art for tsunami early warning.** Due to difficulties in accurately estimating the source extent, conventional approaches – such as seismological methods—encounter difficulties in estimating the tsunamigenic potential of an earthquake. This problem has been addressed by various studies. Convers and Newman³² developed a new method to rapidly discriminate between normal and slow ruptures events, such as those of tsunami earthquakes³,⁷. Kanamori and Rivera² developed a method based on W-phases to determine seismic source parameters for tsunami warning purposes. This method can be used to reliably determine magnitude within ~20 min of the event³³. Recent studies have shown that real-time GPS and seismic data can also be used in a combined approach⁵. This technique allows for rapid determination of the size of the source within few minutes from the rupture³⁴, resulting in a more effective analysis compared to the seismic analysis alone, which suffers from saturation at large magnitudes. Melgar et al.³⁴ tested their algorithm for rapid magnitude estimation with high-rate GPS data for a large number of events, with the Mentawai event being the least successful one. The uncertainties in the magnitude estimation for this particular event highlight again that additional techniques for rapid identification of tsunami earthquakes could improve current warning systems.

Several techniques that imply direct approaches, such as the use of offshore instruments, have been extensively applied for tsunami detection. NOAA developed and deployed the Deep-ocean Assessment and Reporting of Tsunami (DART). DART consists of a sea floor bottom pressure recorder and a moored surface buoy (https://www.ndbc.noaa.gov/dart/dart.shtml). Nevertheless, a DART buoy is primarily suitable to detect a far-field tsunami, as they are deployed and anchored at least 250 miles away from the shore³³ thus being less useful for near-field events. GPS buoys and tide gauges, such as these deployed by the Nationwide Ocean Wave information network for Ports and Harbours (NOWPHAS) near the Japanese coastline, and by the GITEWS, which was operating along the Sumatran coast between 2005 and 2011³⁵, have the potential to provide short time tsunami warning within 10 min of the earthquake. However, in addition to damage caused by weather conditions, floating buoys may be also subject to vandalism, particularly in marine water that is favorable to fishing³⁶. Servicing damaged buoys can be challenging and costly.

Other ocean-bottom seismographic and tsunami observation systems, based on fiber-optics submarine cables, have also been developed for tsunami detection³⁷. These techniques are based on monitoring seafloor earthquakes to detect tsunami. However, these systems come with deployment costs that are usually prohibitively high. More recently, Xerandy et. al.³⁸ proposed a new method based on an underwater communication system composed of fiber optic cables and an undersea network of sensors that could give a 20-min warning time.

It is clear that any additional information that gives rapid identification of the tsunami potential resulting from an earthquake would be useful for the success of tsunami early warning systems. In this study, therefore, we analyze the ionospheric response to two test events of the same seismic magnitude, which provide a unique natural experiment to test if the ionosphere is sensitive to the rupture properties and tsunami genesis, and whether GNSS-TEC observations are able to identify the events that have a high level of tsunami potential.

**Data availability**

The datasets generated during and/or analysed during the current study are available in the DR-NTU repository, https://doi.org/10.21979/N9/3KUEM5.

Received: 22 October 2018; Accepted: 27 May 2020

Published online: 06 July 2020

**References**

1. Kanamori, H. Mechanism of tsunami earthquakes. *Phys. Earth Planet. Int.* **6**, 346–359 (1972).

2. Kanamori, H. & Rivera, L. Source inversion of W phase: Speeding up seismic tsunami warning. *Geophys. J. Int.* **175**, 222–238 (2008).

3. Convers, J. A. & Newman, A. V. Global evaluation of large earthquake energy from 1997 through mid-2010. *J. Geophys. Res.* **116**, B08304 (2011).

4. Duputel, Z., Rivera, L., Kanamori, H. & Hayes, G. W phase source inversion for moderate to large earthquakes (1990–2010). *Geophys. J. Int.* **189**, 1125–1147 (2012).

5. Melgar, D., Bock, Y. & Crowell, B. W. Real-time centroid moment tensor determination for large earthquakes from local and regional displacement records. *Geophys. J. Int.* **188**, 703–718 (2012).

6. Lay, T. et al. The 25 October 2010 Mentawai tsunami earthquake (Mw 7.8) and the tsunami hazard presented by shallow megathrust ruptures. *Geophys. Res. Lett.* **38**, L06302 (2011).

7. Newman, A. V., Hayes, G., Wei, Y. & Convers, J. The 25 October 2010 Mentawai tsunami earthquake, from real-time discriminants, finite-fault rupture, and tsunami excitation. *Geophys. Res. Lett.* **38**, L05302 (2011).

8. Occhipinti, G., Rolland, L., Lognonné, P. & Watada, S. From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes. *J. Geophys. Res. Space Phys.* **118**, 3626–3636 (2013).

9. Occhipinti, G., Kherani, E. A. & Lognonné, P. Geomagnetic dependence of ionospheric disturbances induced by tsunamigenic internal gravity waves. *Geophys. J. Int.* **173**, 753–765 (2008).
10. Occhipinti, G. The seismology of the planet Mongo: The 2015 ionospheric seismology review. Subduction dynamics: From mantle flow to mega disasters. *Geophys. Monogr.* **211**, 169–182 (2016).

11. Occhipinti, G. et al. Three-dimensional numerical modeling of tsunami-related internal gravity waves in the Hawaiian atmosphere. *Earth Planets Space* **63**, 847–851 (2011).

12. Makela, J. J. et al. Imaging and modeling the ionospheric airglow response over Hawaii to the tsunami generated by the Tohoku earthquake of 11 March 2011. *Geophys. Res. Lett.* **38**, L00G02 (2011).

13. Astafyeva, E., Lognonné, P., & Rolland, L. First ionospheric images of the seismic fault slip on the example of the Tohoku-oki earthquake. *Geophys. Res. Lett.* **38**, L22104 (2011).

14. Astafyeva, E., Rolland, L., Lognonné, P., Khelfi, K. & Yahagi, T. Parameters of seismic source as deduced from 1 Hz ionospheric GPS data: Case study of the 2011 Tohoku-oki event. *J. Geophys. Res. Space Phys.* **118**, 5942–5950 (2013).

15. Occhipinti, G., Lognonné, P., Kherani, E. A. & Hébert, H. Three-dimensional waveform modeling of ionospheric signature induced by the 2004 Sumatra tsunami. *Geophys. Res. Lett.* **33**, L20104 (2006).

16. Rakoto, V., Lognonné, P. & Rolland, L. Tsunami modeling with solid Earth–ocean–atmosphere coupled normal modes. *Geophys. J. Int.* **211**, 1119–1138 (2017).

17. Rakoto, V., Lognonné, P., Roll, L., & Coïsson, P. Tsunami wave height estimation from GPS-derived ionospheric data. *J. Geophys. Res. Space Phys.* **123**(5), 4329–4348 (2018).

18. Feng, L. et al. A unified GPS-based earthquake catalog for the Sumatran plate boundary between 2002 and 2013. *J. Geophys. Res. Solid Earth.* **120**, 3566–3598 (2015).

19. Morgan, P. M. et al. The diverse slip behavior of the Banyak Islands section of the Sunda Megathrust offshore Sumatra. *Am. Geophys. Union, Fall Meet. 2015, Abstr. id. T21D-2858* (2015).

20. Hill, E. M. et al. The 2010 Mw7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data. *J. Geophys. Res. Solid Earth.* **117**, B06402 (2012).

21. Satake, K. et al. Tsunami source of the 2010 Mentawai, Indonesia earthquake inferred from tsunami field survey and waveform modeling. *Pure Appl. Geophys.* **170**, 1567–1582 (2013).

22. Cahyadi, M. N. & Heki, K. Coseismic ionospheric disturbance of the large strike-slip earthquakes in North Sumatra in 2012: Mw dependence of the disturbance amplitudes. *Geophys. J. Int.* **200**, 116–129 (2015).

23. Savastano, G. et al. Real-time detection of tsunami ionospheric disturbances with a stand-alone GNSS receiver: A preliminary feasibility demonstration. *Sci. Rep.* **7**, 46607 (2017).

24. Astafyeva, E., Shalimov, S., Olshanskaya, E. & Lognonné, P. Ionospheric response to earthquakes of different magnitudes: Larger quakes perturb the ionosphere stronger and longer. *Geophys. Res. Lett.* **40**, 1675–1681 (2013).

25. Astafyeva, E., Rolland, L. M. & Sladen, A. Strike-slip earthquakes can also be detected in the ionosphere. *Earth Planet. Sci. Lett.* **405**, 180–193 (2014).

26. Calais, E. & Minster, J. B. GPS detection of ionospheric perturbations following the January 17, 1994, Northridge Earthquake. *Geophys. Res. Lett.* **22**, 1045–1048 (1995).

27. Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M. & Webb, F. H. Precise point positioning for the efficient and robust analysis of GPS data from large networks. *J. Geophys. Res. Solid Earth* **102**, 5005–5017 (1997).

28. Rolland, L. M. et al. The resonant response of the ionosphere imaged after the 2011 off the Pacific coast of Tohoku Earthquake. *Earth Planet Space* **63**, 853–857 (2011).

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

30. Yue, H. et al. Rupture process of the 2010 Mw 7.8 Mentawai tsunami earthquake from joint inversion of near-field hr-GPS and teleseismic body wave recordings constrained by tsunami observations. *J. Geophys. Res. Solid Earth.* **119**, 5574–5593 (2014).

31. Suppasri, A. et al. A decade after the 2004 Indian Ocean tsunami: The progress in disaster preparedness and future challenges in Indonesia, Sri Lanka, Thailand and the Maldives. *Pure Appl. Geophys.* **172**, 3313–3341 (2015).

32. Convers, J. A. & Newman, A. V. Rapid earthquake rupture duration estimates from teleseismic energy rates, with application to real-time warning. *Geophys. Res. Lett.* **40**, 5844–5848 (2013).

33. Comfort, L. K., Znati, T., Voortman, M. & Freitag, L. E. Early detection of near-field tsunamis using underwater sensor networks. *Sci. Tsunami Hazards* **31**, 231–243 (2012).

34. Melgar, D. et al. Earthquake magnitude calculation without saturation from the scaling of peak ground displacement. *Geophys. Res. Lett.* **42**, 5197–5205 (2015).

35. Münch, U., Rudloff, A. & Lauterjung, J. Postface “The GITEWS Project - results, summary and outlook”. *Nat. Hazards Earth Syst. Sci.* **11**, 765–769 (2011).

36. C.-C. Teng et al., Buoy vandalism experienced by NOAA National Data Buoy Center. *Ocean. 2009 MTS/IEEE Biloxi Mar. Technol. Futur. Glob. Local Challenges.* **1**, 26–29 (2009).

37. T. Schmitz, W. Rutzen, W. Jokat, In *International Symposium on Underwater Technology, UT 2007—International Workshop on Scientific Use of Submarine Cables and Related Technologies 2007*, 301–304 (IEEE, New York, 2007).

38. Xerandy, X., Znati, T. & Comfort, L. K. Cost-effective, cognitive undersea network for timely and reliable near-field tsunami warning. *Int. J. Adv. Comput. Sci. Appl.* **6**, 224–233 (2015).

## Acknowledgements

Some of the results presented here are reported in the PhD thesis of Dr. Fabio Manta (https://hdl.handle.net/10356/107577), archived by Nanyang Technological University, Singapore. We are grateful to the many scientists and field technicians who have kept the SuGAR network in operation. These include Paramesh Banerjee, Iwan Hermawan, Danny Natawidjaja, Jeffrey Encillo, Imam Suprihanto, Dudi Prayudi, and Bambang Suwargadi. SuGAR GPS data are available for download from ftp://ftp.earthobservatory.sg/SugarData. This research was supported by the National Research Foundation Singapore under its Singapore NRF Investigatorship scheme (Award No. NRF-NRFI05-2019-0009) and by the EOS and the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative. G.O. is supported by Programme National de Télédétection Spatiale (PNTS), grant n PNTS-2014-07; by the CNES, grant GISnet&Back; and by the Institut Universitaire de France (IUF). We thank Benoit Taisne for useful discussions, and the reviewers for constructive comments. This is EOS contribution 300 and IPGP contribution 4054.

## Author contributions

F.M. analyzed and interpreted data, made all the figures, and took the lead in drafting the manuscript. G.O. and E.H. designed the research. L.F. helped to analyze the GNSS data. All authors contributed to drafting of the manuscript.
**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41598-020-68097-w.

**Correspondence** and requests for materials should be addressed to F.M.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020