A possible space-based tsunami early warning system using observations of the tsunami ionospheric hole

Ionospheric plasma disturbances after a large tsunami can be detected by measurement of the total electron content (TEC) between a Global Positioning System (GPS) satellite and its ground-based receivers. TEC depression lasting for a few minutes to tens of minutes termed as tsunami ionospheric hole (TIH) is formed above the tsunami source area. Here we describe the quantitative relationship between initial tsunami height and the TEC depression rate caused by a TIH from seven tsunamigenic earthquakes in Japan and Chile. We found that the percentage of TEC depression and initial tsunami height are correlated and the largest TEC depressions appear 10 to 20 minutes after the main shocks. Our findings imply that Ionospheric TEC measurement using the existing ground receiver networks could be used in an early warning system for near-field tsunamis that take more than 20 minutes to arrive in coastal areas.

The intensity of southward (northward) propagation of the TEC disturbance was clearly visible in the northern (southern) hemisphere, while that of northward (southward) propagation was not clearly observed, because the asymmetrical intensity of northward and southward propagation is caused by the different plasma motion forced to the direction of magnetic field lines 11 .
For the Tohoku EQ, the TEC disturbance began above the tsunami source area and the TEC disturbance originating from acoustic waves subsequently propagated radially. Subsequently, the acoustic and gravity waves generated by Rayleigh waves and tsunami propagations caused the TEC disturbances 9,10,12,13 . In addition to the origin of these disturbances, a TEC depression with a radius of a few hundred kilometres, which remained stationary for approximately 1 hour, was observed 14 . This was termed as the tsunami ionospheric hole (TIH). A TIH was also visible after the Sumatra EQ and the 2010 M8.8 Chile EQ. Further investigations revealed that TIHs appeared after other megathrust and shallow global EQs with magnitudes greater than 7.2 15 . A numerical simulation of the TIH 16 elucidated the plasma dynamics. When the acoustic waves reached the ionosphere, the amplitude of the acoustic waves was intensified due to the very low density of neutral atmosphere. The acoustic waves forced the ionospheric plasma to considerably migrate only along the magnetic field line because of neutral-ion drag. Subsequently, in the decompression phase of the acoustic waves, the downward plasma caused dissociative recombination and suppressed ion production, creating a TEC depression, i.e. TIH. The immediate recombination and slow production caused the TIH to exist for a long time.
Monitoring of tsunami-generated TEC disturbances using GPS receiving network has been expected to be a practical tsunami early warning system, because the disturbance can be detected about 10 minutes after the main shock and the GPS receiving network is improving 17,18 Moreover, the quantitative relationship between initial tsunami height and TEC disturbance amplitude expected from the physical point of view might provide a numerical warning. In fact, Astafyeva et al. 15 reported a correlation not only between initial tsunami height and initial TEC enhancement amplitude, but also between initial tsunami height and TIH duration for 11 large global EQs.
Thus, the quantitative relationship between the initial tsunami height and amplitude of the TEC disturbance immediately above the tsunami source area is required to assess the possibility of using GPS-TEC ionospheric monitoring for a practical early warning system. In this paper, we focus on the value of a TEC depression in a TIH to estimate initial tsunami height, because the TIH remains stationary. We discuss the possibility of constructing an early warning system using only the existing GPS network.

Results
In the analysis herein, we investigate large EQs with accompanying tsunamis in Japan and Chile because the dense and intermediate-dense networks of GPS receiving stations enables us to obtain TEC data above the tsunami source area in magnetic latitudes ranging from 20°N to 35°N in Japan and from 20°S to 35°S in Chile, and the inclination of the magnetic field can be regarded as approximately constant. After commencement of the GPS-TEC observation, seven EQs with sub-ionospheric points (SIPs; the footprint of intersection between the ionosphere modelled as a thin layer and the path of a GPS satellite and a receiver, see information on GPS data in the Method section) around the tsunami source area occurred in Japan and Chile  Table 1). The initial tsunami heights of these EQs were seismologically estimated [19][20][21][22][23][24][25][26][27][28][29][30][31][32] .
To discriminate the TEC trend (mainly originating from the variable length of the slant path between the satellite and the receiver) and the tsunami-generated TEC disturbance including the TIH and the acoustic waves, the difference between the TEC trend and the raw TEC data was analysed. Considering the elevation angle between the ground and the slant path, the residual component (termed Δ vTEC) at the SIP was obtained. After removing the acoustic wave component using a low-pass filter (LPF), the TIH component of Δ vTEC, termed LPF Δ vTEC, was obtained. As an example, spatial images of LPF Δ vTEC before and after the main shock of the Tohoku EQ are illustrated in Fig. 2.
We focused on the percentage TEC depression caused by the TIH and the interval δ T min between the time of the main shock and the minimum value of LPF Δ vTEC (namely the largest TIH). The relationship between the percentage of the TEC depression, δ T min , initial tsunami height and EQ magnitude are illustrated in Fig. 3. As shown in Figs 3a and b, the percentage of TEC depression versus the initial tsunami height are correlated. In contrast, the percentage of TEC depression and the magnitude are not correlated, which is convincing because, in general, the magnitude of the tsunami-generating EQ is not always proportional to the initial tsunami height 33 .
In Fig. 3a, the Maule EQ (D) is slightly large, probably because the night-time small background TEC yielded relatively large error of the percentage of TEC depression, while δ T min of the Maule EQ shown in Fig. 3c seem appropriate. On other hand, in Fig. 3c, the four small EQs with small simulated initial tsunami height (A, B, C, and E) seems inappropriate, probably because of difficult identification for the small depression, while the percentage of TEC depression for them are relatively appropriate. Thus, the EQ with more than 1 m initial tsunami height during daytime might show precise percentage of TEC depression and δ T min .

Discussion
Estimation of the initial tsunami height after large EQs such as the Tohoku EQ takes approximately 20 minutes only from the correlations shown in Fig. 3a and c. This means that an early warning system using TEC measurements may be useful for some coastal areas if 20 minutes plus the necessary time for evacuation is longer than the time taken for the near-field tsunami to arrive at the coast.
The correlations shown in Fig. 3a and c will be useful to verify the simulation of Shinagawa et al. 16 . After parameterization of the simulation, it is possible to predict the initial tsunami height in other regions with different magnetic inclination. As a spatial Δ vTEC map gradually shifts with satellite movement and the number of visible GPS satellites is limited, the Δ vTEC map does not provide continuous coverage of an entire area even in Japan. As the number of other satellites and receivers of Global Navigation Satellite Systems (GNSS) such as GLONASS, BeiDou, and Galileo and Regional Navigational Satellite System (RNSS) such as QZSS are rapidly increasing, TEC measurements of GPS and other GNSS/RNSS are expected to continuously cover the entire area in the near future. In addition, real-time GPS data is currently provided in many countries. In the present  methodology, Δ vTEC can be immediately calculated even using personal computer after obtaining real-time slant TEC. This yields a feasible early warning system for the near-field tsunami.

Method
GPS data. GPS is a satellite navigation system that provides precise location and time information derived from the phase differences of carrier radio signals (1575.42 and 1222.60 MHz) between several satellites and their receivers. Ionospheric plasma, the dispersive medium for radio signals, causes a frequency-dependent phase delay during signal propagation. The phase difference between two carrier signals measured at the receivers is proportional to the TEC between the satellite and the receiver, termed the slant TEC. The unit of TEC value is TECu, corresponding to 1.0 × 10 16 electron m −2 . The altitude of the F region, which is the major contributor to TEC, is approximately 200 to 400 kilometres. For the sake of simplicity, we approximated the region as a thin layer at 300 kilometres. The point at which the slant path from the satellite to the receiver intercepts this ionospheric layer is termed the ionospheric point. In Japan, GPS data are provided by the Geographical Survey Institute of Japan, which has installed more than 1200 receivers in a nationwide GPS array, named the GPS Earth Observation Network (GEONET; ftp://terras.gsi.go.jp/). Its sampling period of GPS data is 30 seconds or 1 second. GEONET provides spatial-temporal ionospheric images. In Chile and Argentina, GPS data are provided by UNAVCO (ftp://data-out.unavco.org/) and the Instituto Geográfico Nacional (http://www.ign.gob.ar/NuestrasActividades/ Geodesia/Ramsac/DescargaRinex), which has installed several tens of receivers in a GPS array, having sampling periods of GPS data between 15 and 60 seconds.

Percentage of TEC depression in TIH.
To derive the TEC disturbance caused by the TIH, the following procedure was applied. First, the least-squares quadratic curve of the slant TEC time-series in the period from 30 minutes prior to 40 minutes after the main shock was obtained, except for the Tohoku EQ and the Maule EQ (Fig. S1a). For the Tohoku and Maule EQs, only the period from 30 minutes prior to 7 minutes after the main shock was used because the TIH was too large to obtain appropriate approximation after the TIH. Peculiar bias components for each satellite and the receivers were reduced by the difference between this fitting curve and the raw data. Second, the difference in the vertical TEC (i.e. Δ vTEC) was derived by multiplying the cosine between the slant direction and the satellite earth direction by the difference between the fitting curve and the raw data (Fig. S1b).
A spectrum analysis of Δ vTEC in the Tohoku EQ for all stations is illustrated in Fig. S2. The period of the spectrum analysis was 4096 seconds, consisting of 1024-second data of Δ vTEC located in a constant part of a boxcar window function. In Fig. S2, the vertical axis indicates the distance between the SIP at the time of the main shock and the centre of the tsunami source area, and the horizontal axis indicates the period of the spectrum. Several resonant modes of acoustic waves occurred; in particular, they were intense in the period from 150 to 210 seconds above the tsunami source area. According to Watada and Kanamori 34 , the observed distinct patterns from 150 to 210 seconds are inferred to be resonant modes of the acoustic waves excited by the tsunami. To exclude the variations of resonance modes including the initial TEC pulse, i.e. the acoustic wave components, and extract the component of the TIH, the ±150 sec running mean of Δ vTEC except the Niigataken Chuetsu-oki EQ was taken, i.e. LPF Δ vTEC (Fig. S2b). For the Niigataken Chuetsu-oki EQ, ± 90 sec running mean is used because of a small TEC disturbance.
The following procedure was employed to estimate the value of the TEC depression of the TIH above the tsunami source area. First, the tsunami source area was defined to be the average location of the Tohoku EQ (N38.0, E143.4), the Maule EQ (S34.8, W72.9), and the Illapel EQ (S30.9, W72.3) reported by several articles [28][29][30][31]26,27,32 , and equivalent to the epicentre of other small EQs. Time-series of LPF Δ TEC to estimate the TEC depression were selected such that the SIP at the time 7 min after the main shock was located within 100 kilometres from the tsunami source area of Mw ≥ 8 EQs and 50 kilometres from the epicentre of the other EQs from 7 minutes to 24 (the Tohoku and Maule EQs) and 16 (the other EQs) minutes after the main shock (Fig. S1). We defined δ vTEC TIH as the difference between the value of the LPF Δ vTEC 7 minutes after the main shock and the minimum value of LPF Δ vTEC, and δ T min as the interval between the main shock and the minimum time of LPF Δ vTEC. Considering the simulation of Shinagawa et al. 16 , the TEC depression of TIH was assumed to be proportional to the electron density in the F2 region immediately before the TEC disturbance, because the depressed plasma migrates downward from the F2 region. In this study, the percentage TEC depression with respect to absolute vTEC, i.e. δ vTEC TIH /vTEC × 100, was discussed for comparison with different solar conditions such as during different seasons and at different local times. For derivation of absolute vTEC, data of Global Ionosphere Maps Total Electron Content (GIMTEC: ftp://ftp.unibe.ch/aiub/CODE/) was used. Finally, the average values and standard deviations of the percentage TEC depression (δ vTEC TIH /vTEC × 100) and δ T min for each EQ from the selected time-series of LPF Δ vTEC (Figs S1 and S2) were obtained.
Geomagnetic storm. Ionospheric storm caused by geomagnetic storm causes TEC variation. If the dominant quarter period of TEC variation during ionospheric storm is much longer than 15 min corresponding to the interval between the acoustic wave arrival at the ionosphere and the minimum of TEC caused by the TIH, the percentage TEC depression is not significantly affected by the ionospheric storm because the background TEC values during the interval are regarded as the same. In fact, the dominant period of TEC variation during ionospheric storm is on the order of hours 35 . In contrary, the TEC variations with the period of less than 300 s are reduced by the running mean process. However, the periods for more than 300 s should be taken into account. Only in the Tohoku EQ for the present analysis, Dst index was large (Table 1). Spectrum analysis shown in Fig. S2 indicated no significant variation for period larger than 300 s. Earthquakes with tsunami. The earthquakes accompanying tsunami listed in Table 1 are selected as follows: 1) Geomagnetically mid-latitude are selected because the numerical simulation 16 was applied for these area. 2) Accessible GPS data has been provided. 3) At least one SIP 7 min after the mainshock existed within 50 km and 100 km epicentral distances for M < 8 and M ≥ 8. 4) The initial tsunami height (> 0.1 m) was reported in the refereed journal except the largest foreshocks of the Tohoku EQ, which is reported below. Initial tsunami height. Average values and standard deviations of the initial tsunami height were obtained from sources listed in references [19][20][21][22][23][24][25][26][27][28][29][30][31][32] . Most studies were consistent with our results, comparing tsunami heights observed at several coastal tide gauges and offshore GPS buoys with the simulated values and defining the initial tsunami height as the vertical deformation from fault motion.
Tsunami waveforms recorded at tide gauges and offshore tsunami gauges can be used in tsunami waveform inversion to estimate the initial tsunami height in the source area. For the case of the largest foreshock of the 2011 Tohoku EQ that occurred on 9 March, five tsunami waveforms at two offshore pressure gauges and three GPS buoys were used to estimate the slip distribution of the largest foreshock 36 . This slip distribution has a major slip region (measuring 45 km × 45 km) with slip amounts larger than 0.5 m. The seismic moment calculated from the slip distribution is 1.2 × 10 20 Nm (Mw 7.3), which is consistent with the United States Geological Survey Wphase solution. The maximum initial tsunami height calculated from the slip distribution is approximately 0.3 m (Fig. S5).