The source detection of 28 September 2018 Sulawesi tsunami by using ionospheric GNSS total electron content disturbance

The 28 September 2018 magnitude Mw7.8 Palu, Indonesia earthquake (0.178° S, 119.840° E, depth 13 km) occurred at 10:02 UTC. The major earthquake triggered catastrophic liquefaction, landslides, and a near-field tsunami. The ionospheric total electron content (TEC) derived from records of 5 ground-based global navigation satellite system (GNSS) receivers is employed to detect tsunami traveling ionospheric disturbances (TTIDs). In total, 15 TTIDs have been detected. The ray-tracing and beamforming techniques are then used to find the TTID source location. The bootstrap method is applied in order to further explore the possible location of the tsunami source based on results of the two techniques, which show the beamforming technique has a slightly better performance on finding possible locations of the tsunami source. Meanwhile, the circle method is employed to examine tsunami signatures of the sea-surface height and video records, and find possible tsunami origin locations. The coincidence of the TTID source location and the tsunami location shows that the ionospheric TEC recorded by local ground-based GNSS receivers can be used to confirm the tsunami occurrence, find the tsunami location, and support the tsunami early warning.


Introduction
Seismic waves and tsunami of a large earthquake can significantly induce traveling atmospheric disturbances (TADs) near the Earth's surface in the epicenter and the tsunami source areas ). These TADs propagate vertically at the local sound speed traveling from the atmosphere into the ionosphere, forming traveling ionospheric disturbances (TIDs) and simultaneously inducing internal gravity waves within. Liu et al. (2019a) examined ionospheric electron density profiles away from the epicenter and the tsunami source areas during the 2011 M9.0 Tohoku earthquake, and found that underneath Rayleigh and tsunami waves can induce more prominent TIDs. Data recorded from ground-based global navigation satellite system (GNSS) receivers have been used to observe ionospheric total electron content (TEC) disturbances triggered by seismic waves (Jung et al. 2006;Liu et al. 2010Liu et al. , 2012Liu et al. , 2016Sun et al. 2016;Tsai et al. 2011) and tsunami waves (Artru et al. 2005;Liu et al. 2006aLiu et al. , b, 2011Liu et al. , 2019aTsai et al. 2011;Galvan et al. 2011;Kamogawa et al. 2016;Komjathy et al. 2012). Since the horizontal speed of seismic waves of 2.5-3.5 km/s very differs from that of tsunami waves of 200-250 m/s, we can clearly identify and easily separate tsunami-induced and earthquake-induced TEC disturbances. Liu et al. (2006a) employ the ray-tracing technique (Aki and Richards 2002) to estimate the arrival times at locations of tsunami traveling ionospheric disturbances (TTIDs) of TEC for finding the tsunami source, as well as to see whether the observed disturbances of the ionospheric Open Access *Correspondence: tigerjyliu@gmail.com 1 Center for Astronautical Physics and Engineering, National Central University, Taoyuan, Taiwan Full list of author information is available at the end of the article GNSS TEC is triggered by the tsunami or not during the 26 December 2004 M9.3 Sumatra earthquake (Liu et al. 2006a). They reported that the source of the 26 December 2004 Indian Ocean tsunami was about 580 km southwest from the Sumatra epicenter. Recently, Savastano et al. (2017) showed the possibility of detecting TTIDs in a real-time scenario. Liu et al. (2019a) applied the ray-tracing technique (Lee and Stewart 1981), and the beamforming technique (Huang et al. 1999) to analyze ionospheric TEC disturbances of the 2004 Indian Ocean tsunami, showing that the TTIDs can be rapidly identified within 20 min after the earthquake allowing the tsunami source to be located. They further proposed to use existing ground-based GNSS (including GPS, GLONASS, Galileo, and BeiDou) receiving stations, which have been routinely operated by International GNSS Service (IGS), to monitor TTIDs in real time for tsunami early warning applications.
On 28 September 2018, a Mw 7.8 earthquake occurred at 10:02:44 UT at 0.18° S, 119.84° E in the mountainous Donggala Regency, Central Sulawesi, and generated a tsunami (http://itic.ioc-unesc o.org/index .php). International Tsunami Information Center reported that 5 min after the earthquake, the Indonesia Agency for Meteorology, Climatology and Geophysics (Badan Meteorologi Klimatologi dan Geofisika-BMKG) issued a tsunami warning for a local tsunami. Recently, a retrieved marigram from the Pantoloan-Sulteng tide gauge (0.695° S, 119.827° E) shows a 3.8-m trough-to-peak tsunami that arrived 6 min after the earthquake occurrence. There were no other nearby instrumental observations. Several analyses based on pictures, post-disaster information, and video clips suggest the first tsunami wave hit the Palu beach area (0.885° S, 119.857° E) 7-12 min after the earthquake (Heidarzadeh et al. 2019). Preliminary field surveys conducted by the BMKG and IRIDeS/Indonesia Ministry of Environment and Forestry/Chuo University report an eyewitness height of up to 11.3 m in Palu and 3-10 m on the west and east sides of the bay. Detailed field surveys can be found in Omira et al. (2019). However, based on the above information, it is still difficult to locate the tsunami's source.
In this paper, the time delay method, the circle method, the ray-tracing technique, and the beamforming technique are employed. Figure 1 illustrates locations of the Sulawesi epicenter and 15 TTID signatures identified from TEC measurements of 5 ground-based GNSS receiving stations. The time delay method is applied on distances and arrival times of the TTIDs to estimate the associated horizontal speed and onset time, which are further given, respectively, to the ray-tracing technique and the beamforming technique to locate the tsunami's source area. A statistical bootstrap analysis is also conducted to evaluate the reliability of the computed tsunami's source location.

Ionospheric tsunami disturbance and source location
Tsunami waves create mechanical disturbances in the near sea-surface atmosphere, which propagate upward into the ionosphere, interact with the ionized gas within, and disturb the ionospheric plasma, electron density, TEC, etc. In general, it takes 8-10 min for disturbances activated by seismic or tsunami waves traveling from the earth's surface upward into the ionosphere at 300 km altitude, indicating that the TTID upward propagation speed is of about 500-625 m/s, which is near the averaged sound speeds in the atmosphere of between the ground and 300 km altitude (Liu et al. 2006a(Liu et al. , b, 2011(Liu et al. , 2019b. In this study, for simplicity, we assume the TTID average upward speed is 600 m/s. When radio signals transmitted from GNSS satellites at about 20,200 km altitude pass through the ionosphere, and reach ground-based receivers, they will experience changes in propagation speed and refraction from the ionospheric plasma, electron density, TEC, etc. Since the largest electron density is at the F2-peak, which is usually located at about 250-350 km altitude, the ionosphere then can be treated as a thin shell at about 325-375 km altitude (Sardon et al. 1994;Liu et al. 1996). From recorded broadcast ephemeris and assuming the thin-shell ionosphere at 325 km (Davies 1990;Hofmann-Wellenhof et al. 1992;Liu et al. 1996;Sardon et al. 1994;Wu 1992), the intersection point, which is termed the ionospheric pierce point (IPP), between the signal path from a GNSS satellite to a ground-based receiver and the thin-shell ionosphere give the longitude and latitude location of TTIDs. Thus, the TEC at the IPP acts as a space tide gauge or a space buoy floating at 325 km altitude to detect TTIDs. Figure 1a displays locations of Sulawesi epicenter and the 15 TTIDs of the 5 GNSS ground-based receivers. To detect TTIDs, we first apply the time delay method (Liu et al. 2006a(Liu et al. , 2019aTsai et al. 2011) to estimate the possible arrival time range of TEC fluctuations observed at each GNSS receiving station, and then identify the precise arrival time. For example, the IPP point observed kat1 21 was at 2135 km from epicenter, and assuming the tsunami horizontal speed of about 190-300 m/s, it requires the travel time of about 2.1 (= 2135 km/300 m/s + 8 min) to 3.3 (= 2135 km/190 m/s + 8 min) hours, where 8 min of the vertical travel time should be included. We further add up the travel time and the earthquake time, and search the possible arrival time range between 12.1 and 13.3 UT. The precise arrival time (12.28 UT) of TTIDs is identified, when the TEC starts drastically changing. Figure 1b illustrates 15 TEC variations filtered using a high-pass filter with the 20-min cut-off period, the arrival time of each TTID, and the arrival time versus the epicenter distance (i.e., the distance between the epicenter and each observed TTID) (for detail see, Table 1). The average horizontal propagation speed of the TTIDs usually is close to that of the typical speed of tsunami waves of 200-250 m/s (Dean and Dalrymple 1984;Dias and Dutykh 2007;Murty 1977) in the open sea. A slope of the distance versus the arrival time of the 15 TTIDs shows the average horizontal speed being about 228 m/s. To further search the TTID's source location, we apply the ray-tracing technique and beamforming technique on the arrival time and the location of observed TTIDs (Table 1). The ray-tracing technique (Lee and Stewart 1981) and the beamforming technique (Huang et al. 1999) are used to globally search for the TTID's source location inside a rectangular area between ± 30° N and 70-150° E, by shifting along a 0.1° grid in the longitudinal and then latitudinal direction or vice versa. In the ray-tracing technique, for a given velocity model, we can guess a source location, and calculate the travel time from the guessed source location to each TTID, which is the distance between the guessed source location and the TTID location divided by the given velocity, including the upward speed and the horizontal speed. We then subtract the calculated travel time from the arrival time to obtain the origin time for each TTID. Since there are 15 origin times, their associated average and standard deviation at the guessed source location can be computed. We then repeat the above process by shifting the guessed source location 0.1° in the longitudinal direction and then the latitudinal one to cover the entire rectangular area. Finally, we contour the standard deviations of the origin time in the rectangular area to find the location of the minimum, where can be considered as the TTID source location. Since the most prominent TTID signature generally is triggered by its underlying tsunami wave (Liu et al. 2019a), the location of the TTID source should be very close to that of the tsunami origin but the appearance time with 8-10 min delay. Assuming the thin-shell ionosphere at 325 km altitude, the averaged vertical speed of 600 m/s, and the horizonal speed of 228 m/s, the ray-tracing yields that the minimum of the standard deviation is ± 14.3 min with the averaged origin time of 10:02 UTC located at 0.7° S, 119.6° E, where is 63 km southsouthwest of the epicenter and the most likely location of the TTID or tsunami source (Fig. 2a). In fact, various horizonal speeds of 200-250 m/s have been tested, and the speed of 228 m/s results in the TTID onset time of 10:02 UTC. This suggests that the ray-tracing technique is very sensitive to the velocity model.
Similarly, in the beamforming technique, for a given TTID origin time at 10:02 UTC, we can guess a source location, find the distance, and calculate the travel time, as well as compute the averaged velocity and standard deviation of the 15 TTIDs for each guessed source location. Again, we plot contours of standard deviations of the velocity in the rectangular area to find the minimum, which again can be considered as the TTID source location. Assuming the thin-shell ionosphere at the same 325 km altitude and an averaged vertical travel time of 542 s (i.e., 325 km divided by 600 m/s), the beamforming technique finds that the minimum of the standard deviation of ± 23.7 m/s with the averaged horizontal speed of 227.5 m/s is located at 0.6° S, 119.6° E, which is about 53 km south-southwest of the epicenter, to be the most probable the TTID and/or tsunami source location (Fig. 2b). Note that the averaged horizontal speed of 227.5 m/s computed by the beamforming technique is very close to that the 228 m/s assumed by the ray-tracing technique.

Statistical analyses of TTID source location
It has been demonstrated that the possible location of the tsunami source can be found by using both the ray-tracing and beamforming techniques based on the 15 TTIDs. To further explore the possible location of the tsunami source, the bootstrap method (Efron and Tibshirani 1986) is used to find the confidence region of the source location centered at the original detected TTID source location. In the bootstrap method, 1000 samples of size 15 are resampled with replacement from the 15 TTIDs to obtain the associated 1000 TTID's source locations based on either the ray-tracing or beamforming technique (yellow dots in Fig. 2c, d).
Note that the bootstrap source locations show a northeast-southwest-ward trend and have different variations in longitude and latitude. The covariance matrix of the longitude and latitude of bootstrap source locations is then computed, hence the Mahalanobis distance (Dillon and Goldstein 1984) between each bootstrap source location and the original detected source location is calculated. Moreover, the Mahalanobis distances are ordered from the smallest to the largest so that the 90-percentile of the distances is obtained. The 90% confidence regions for the possible tsunami source locations detected by the ray-tracing and beamforming techniques (white ellipses in Fig. 2c, d) based on the associated 90-percentile of the distances are finally constructed which can be, respectively, expressed as: and where x and y denote the longitude and latitude. The 90% confidence region of the beamforming techniques is slightly more concentrated than that of the ray-tracing technique.

Discussion and conclusion
Scientists (Liu et al. 2006a(Liu et al. , b, 2011(Liu et al. , 2019aTsai et al. 2011) have found that the distance between the epicenter and the tsunami origin location can be on the order of hundreds km. For the 2011 Tohoku tsunami,  using dense IPPs of 12,000 space buoys and/or tide gauges directly imaged the tsunami origin location, 0.5787(x − 119.5) 2 +2.7379(y + 0.8) 2 − 1.8358(x − 119.5)(y + 0.8) ≤ 3.9820, 1.3100(x − 119.6) 2 +3.0104(y + 0.6) 2 − 2.6814(x − 119.6)(y + 0.6) ≤ 4.6647, which is 200 km eastward of the epicenter, while Tsai et al. (2011), assuming a horizontal speed of 210 m/s and applying the ray-tracing technique on a limited number of 27 TTIDs, found that the tsunami origin was 195 km southeast of the epicenter. Thus, the discrepancy of the tsunami origin locations between the two studies might result from the limited number of 27 TTIDs being used by Tsai et al. (2011). Additionally, the ray-tracing and beamforming techniques have been employed to find the most prominent vertical surface motion caused by Rayleigh waves of earthquakes (Jung et al. 2006;Liu et al. 2010Liu et al. , 2012, and however, the beamforming technique has not yet been often applied on TTIDs in detail. Liu et al. (2010) found that results of the ray-tracing technique are a function of the velocity model, while the beamforming technique can successfully locate the most prominent vertical surface motion of the Chi-Chi earthquake. Note that the horizontal speed might significantly disturb the detection of the TTID source location. Several factors, such as temperatures, densities, wind speeds, etc., in the atmosphere, can easily affect TTID propagation speeds, and therefore it is difficult to issue a suitable and/ The white ellipses represent the 90% confidence interval for the possible tsunami source locations or correct velocity model. By contrast, the beamforming technique assumes the origin time of TTIDs near the earth's surface being close to that of the co-located tsunami waves, which is the same as the earthquake occurrence time. It might be the assumption being suitable and close to the reality, the beamforming technique can quickly locate the TTID and/or tsunami location. Figure 2c, d reveals that the 90% confidence region of the beamforming technique is more concentrated than that of the ray-tracing one, which suggests that the beamforming technique might yield a better performance in locating the source origin of TTIDs or tsunami waves. To further find if the TTID source locations derived from the ray-tracing and beamforming techniques are close to the tsunami origin location, we apply the circle method (Lay and Wallace 1995; Liu et al. 2006b) on data of the Pantoloan-Sulteng tide gauge (0.695° S, 119.827° E) that the tsunami arrived at 6 min after the earthquake occurrence, and of video records showing the time that the first tsunami wave hits the Palu beach area (0.885° S, 119.857° E), roughly 7-12 min after the earthquake. Assuming a tsunami wave with 80 m/s in Palu bay (let the bay depth of about 700 m), we apply the circle method on the tsunami arrival time of 6 min and 7-12 min after the earthquake and compute and draw two rings of possible tsunami source locations centering at Pantoloan-Sulteng tide gauge and at the Palu beach area, respectively. Figure 3 displays the ring centered at Pantoloan with radii of the inner edge of 24 km (= 80 m/s × 5 min) and the outer edge of 29 km (= 80 m/s × 6 min) km, and the ring centered at the Palu beach area with radii of the inner edge of 34 km (= 80 m/s × 7 min) and the outer edge of 58 km (= 80 m/s × 12 min) km. The intersection area of the two rings can be considered as the tsunami source location. Since the two TTID source locations fall inside of the intersection area, both the beamforming and raytracing techniques are suitable to find the tsunami origin location.
The source of 2018 Sulawesi tsunami is still under debate. Heidarzadeh et al. (2019) used tsunami raytracing method to find out the potential landslide locations inside the Palu Bay. Omira et al. (2019) identified locations of the coastal landslides by combing the data of video analysis, field survey, and satellite images. Song et al. (2019) proposed a two-segment fault model by using geodetic data. By contrast, our results of the ray-tracing technique, the beamforming technique, and the circle method show that the possible tsunami source locations are near but just outside the Palu Bay. The discrepancy might result from the fact that the limited number of the 15 TTIDs and that the used ground-based GNSS receivers are too far away from the Sulawesi area. It shall be able to drastically shorten the time in detecting TTIDs and confirming a tsunami occurrence, as well as significantly improve the precision in locating the tsunami source, if we make GNSS TEC data derived from existing nearby and/or local ground-based GNSS receivers available in real time.
In conclusion, the ionospheric TEC derived from 5 ground-based GNSS receivers in the South Asia region detected 15 TTIDs induced by the 28 September 2018 Sulawesi tsunami and shed some light on the tsunami source location.