12 Jul 2021
12 Jul 2021
Robust uncertainty quantification of the volume of tsunami ionospheric holes for the 2011 TohokuOki Earthquake: towards lowcost satellitebased tsunami warning systems
 ^{1}Department of Statistical Science, University College London, London, UK
 ^{2}Department of Space and Climate Physics, University College London, London, UK
 ^{3}Global Center for Asian and Regional Research, University of Shizuoka, Shizuoka, Japan
 ^{4}Institute of Oceanic Research and Development, Tokai University, Shizuoka, Japan
 ^{1}Department of Statistical Science, University College London, London, UK
 ^{2}Department of Space and Climate Physics, University College London, London, UK
 ^{3}Global Center for Asian and Regional Research, University of Shizuoka, Shizuoka, Japan
 ^{4}Institute of Oceanic Research and Development, Tokai University, Shizuoka, Japan
Abstract. We develop a new method to analyze the total electron content (TEC) depression in the ionosphere after a tsunami occurrence. We employ Gaussian process regression to accurately estimate the TEC disturbance every 30 s using satellite observations from the GNSS network, even over regions without measurements. We face multiple challenges. First, the impact of the acoustic wave generated by a tsunami onto TEC levels is nonlinear and anisotropic. Second, observation points are moving. Third, the measured data is not uniformly distributed in the targeting range. Nevertheless, our method always computes the electron density depression volumes, along with estimated uncertainties, when applied to the 2011 TohokuOki Earthquake, even with random selections of only 5 % of the 1,000 GPS Earth Observation Network System receivers considered here over Japan. Also, the statistically estimated TEC depression area mostly overlaps the range of the initial tsunami, which indicates that our method can potentially be used to estimate the initial tsunami. The method can warn of a tsunami event within 15 minutes of the earthquake, at high levels of confidence, even with a sparse receiver network. Hence, it is potentially applicable worldwide using the existing GNSS network.
 Preprint
(3842 KB) 
Supplement
(1603 KB)  BibTeX
 EndNote
Ryuichi Kanai et al.
Status: open (extended)

RC1: 'Comment on nhess2021119', Anonymous Referee #1, 06 Aug 2021
reply
The paper developed the method to identify ionospheric hole generated by tsunami more precisely than the previous method. I am interested in their method and believe that this method is scientifically important. There are some problems need to be answers clearly before the paper is published.
The major comments
 They used sparse data where 95 % of the GNSS receivers are randomly removed from the observed data. However, they only show one example of the sparse data set. If you randomly removed 95 % of data. You can generate a large number of sparse data sets. Therefore, you can analyze how the variation of data sets affects to the results. If you remove 95 % of data randomly, some of your data sets may have only a few receivers in the ionospheric hole. We want to know how those data sets affect to the results.
 In all of the maps in Figures, the position of the Japan trench should be shown because we all knew that the tsunami initial surface uplift of the 2011 Tohokuoki earthquake was located at landward from the trench because it was underthrust earthquake, Therefore, the ionospheric hole was also better to be located at landward from the trench. By looking at Figure 8, a part of ionospheric hole at 6:12:00 and 6:16:00, (a2), (a3), (b2), and (b3), are located at oceanward for the Japan trench. Please explain reasons for those. At 6:08:00 is 21 minutes after the earthquake, the ionospheric hole is still the same as the initial tsunami surface uplift zone (a1and b1 in Figure 8). What are reasons that the hole increased the areas to seaward (eastward) at 25 and 29 minutes after the earthquake?
 I am sure that it is important to identify the initial uplift area for tsunami early warning purpose. However, it takes 2029 minutes to estimate the area. Therefore, I believe that this method should be more effective and powerful by combining the existing method or some other method recently developed. The authors should discuss those in the paper.
Minor comments.
 In page 3, “Is the assumption of altitude of 300 km sufficient for day and night times?”
 In page 4, “6:46:30 (UTC) and 6:46:18(UTC)” should be “5:46:30 (UTC) and 5:46:18(UTC)”
 In Figure 5, three dimensional plots (c, d, e, and f) are difficult to find exact the confidence levels. They may be better to plot in the map views (2D) with some crosssections.
 Please explain clearly how do you chose the locations of eight triangles in Figure 6 and locations of each time at eight directions in Figure 7.

AC1: 'Reply on RC1', Ryuichi Kanai, 03 Sep 2021
reply
Thank you very much for your comments.
Our answers are as below.
[About the first major comment.]
In this study, the surface fitting for the sparse data is performed using the data received by the remaining 5% of receivers after randomly excluding 95% of the GNSS receivers.
Therefore, in principle, it is possible that randomly chosen 5% of the receivers may not be able to observe the data points that exist in the TIH region.
In our analysis, the minimum number of GNSS receivers within the target area detecting data from the satellite between 5:46:30 and 6:16:30, which is the period of analysis, is 832 (at 5:46:30), and 40 receivers were randomly selected as a conservative estimate of 5% of that number.
Consider the aforementioned case where only receivers that do not measure TIH are randomly selected to make up these 40 units.
For example, at 6:12:00, there are 66 receivers receiving data with a TEC value of 4 or less, and 25 receivers receiving data with a TEC value of 5 or less.
At 6:12:00, there are 917 receivers detecting data from the satellite in the target area, and the probability that at least one of the 40 receivers randomly selected from these receivers contains one of 66 receivers that receive 4 or less TEC value is around 95%.
Also, the probability that at least one of those 40 receivers contains one of 25 receivers that receive 5 or less is approximately 70%.
Therefore, a random selection of receivers that cannot measure TIH can happen though the probability of that is very low.
In this experiment, we analyze the usefulness of surface fitting in the sparse case, where there are receivers that measure TIH.
The specific details of the situation of this experiment are described below.
Ten experiments were conducted to randomly select 5% of receivers, and the minimum observed values measured in each case are shown in the table below.
With a rare exception, the receiver number by which the minimum value is measured in each case is different.
For example, at 06:12:00, in the case of random9, the minimum value is about 0.8 TECu larger than the case with all data.
Also, at 06:16:00, the minimum value of the full data is 5.19, but there is only one case, random5, where the minimum value is less than 5.
How TIH volume, i.e. the volume of the area where the value of TEC is less than 0 in the target area, the calculation is affected in this situation is analyzed in Figure below.
The volume of TIH is used as a measure to determine the effect.
It can be seen from the Figure below that in each case, there are very few measurement points with a value of 4 or lower, and even the number of measurement points with a value of 3 or less is fewer than 10.
However, when compared to the volumes calculated using all the data represented by the horizontal lines, it can be seen that each case has a very close value at each time. Furthermore, although the number of data points below 2 varies, the computed volumes are almost similar.
Obviously, if there is no data point in the TIH area, volume calculation itself is impossible. However, if data points shown in the figure are measured in the TIH, volume calculation is possible. Since a single receiver can receive data from multiple satellites, the number of receivers needed to receive data in TIH is extremely small.
[About the second major comment.]
Japan Trench data has been added to all 2D figures.
The figure below shows the uncertainty in estimating the values of TEC using the full data.
In this figure, the uncertainty is defined as three times the standard deviation.
In general, interpolation between data points can be performed with a small uncertainty, while extrapolation has a larger uncertainty.
Even in the case of interpolation, if the data points are sparse, the uncertainty will be large.
Therefore, it can be seen from the Figure that the uncertainty is larger in areas where the measurement points are sparse or where extrapolation is performed.
In this regard, for example, looking at the uncertainty values for each region listed in Figure below, at 6:12:00, the uncertainty for the region east of the initial tsunami region and west of the dotted line at longitude 145.37 is only slightly larger.
However, it is unlikely that this was the reason for capturing the nonoverlapping phenomenon of the initial tsunami and TIH regions.
Heki and Ping (2005) shows that acoustic waves propagate upward, which are gradually refracted, and their effects propagate horizontally in the ionosphere.
According to this principle, TIH is expected to spread evenly in the eastwest direction of the tsunami generation area.
Kakinami et al. (2012), who analyzed the measured data, showed that Slant TEC decreased in the area east of the tsunami generation area after the Great East Japan Earthquake.
The reason why the TIH calculated by our method appears to be extended to the east of the Japan Trench when compared to the initial tsunami area is that the initial tsunami height is highest on the Japan Trench.
The acoustic waves from the highest region on the trench propagate into the atmosphere and affect the neutral atmosphere evenly in the eastwest direction, resulting in the recombination of ions and electrons, causing the TIH to spread in the eastwest direction around the trench.
Therefore, we believe that the estimation by our method correctly captures the variation of the electron density.
In future reverse calculations of the initial tsunami area and height based on TIH information, it is necessary to take into account that the area with the highest initial tsunami has a large impact on TIH formation.
 Heki, K. and Ping, J.: Directivity and apparent velocity of the coseismic ionospheric disturbances observed with a dense GPS array, Earth and Planetary Science Letters, 236, 845–855, https://doi.org/https://doi.org/10.1016/j.epsl.2005.06.010, 2005.
 Kakinami, Y., Kamogawa, M., Tanioka, Y., Watanabe, S., Gusman, A. R., Liu, J.Y., et al.: Tsunamigenic ionospheric hole, Geophysical Research Letters, 39, https://doi.org/https://doi.org/10.1029/2011GL050159, 2012.
[About the third major comment.]
Based on the results of previous studies, which show that a small tsunami cannot form a large TIH while a large tsunami can form a large TIH, we believe that we can issue an alarm that a large tsunami is occurring when the volume of the TIH reaches a certain threshold value.
Further analysis is needed, but we can conclude that the timing of crossing this threshold is quite early, at least based on the current analysis of this study.
We also believe that our method can estimate the initial tsunami independently from the previous methods.
Previous methods include, for example, the method of estimating tsunami from the source and magnitude of an earthquake.
Seismic waves are transmitted as Pwaves (longitudinal tremors) of about 7 kilometers per second and Swaves (lateral tremors) of about 4 kilometers per second. By detecting seismic waves at multiple observation points, it is possible to calculate the approximate location of the earthquake and to some extent the exact size and magnitude of the earthquake within two minutes.
Based on this information, the initial tsunami can be estimated.
In addition, a method has been developed to calculate the initial shape of a tsunami in reverse by calculating the shift of vessel speed due to the passage of a tsunami from the data of the Automatic Identification System (AIS), which is required to be installed on ships sailing in the area. This method is expected to be used to predict the initial tsunami around the world.
However, some problems exist, such as the fact that the exact magnitude of a huge earthquake with a magnitude greater than 8 is not immediately known, the fault displacement (estimated from seismic waves) does not always match the initial sealevel change, and the initial sealevel change cannot be known.
In addition to the above methods, Japan has developed a system called REaltime GEONET Analysis system for Rapid Deformation monitoring (REGARD), which analyzes GEONET data in realtime and extracts crustal deformation during earthquakes to automatically estimate fault models and earthquake scale within 3 minutes after the earthquake.
In addition, the Seafloor Observation Network for Earthquakes and Tsunamis along the Japan Trench (SNET) has been established on the seafloor from off Boso to off Tokachi.
The data is collected in realtime 24 hours a day.
This system is expected to be used to accurately predict the initial tsunami.
As a matter of course, it is difficult to accurately estimate the initial tsunami information because no estimation method is perfect and the area covered may be limited.
Therefore, in addition to the various existing initial tsunami estimation methods mentioned above, if the initial tsunami shape and height estimation based on our developed TIH estimation can be realized in the future, it is expected that the combination of these methods will enable us to realize even more accurate initial tsunami height and range estimation.
From this point of view, even if it takes about 20 minutes to obtain the initial tsunami information, we believe it is beneficial.
Furthermore, the existence of the second and third waves can be estimated by looking at the time variation of TIH.
Since the presence of large second and third waves affects the shape of TIH, it can be used to determine whether the tsunami warning should be maintained or canceled after it is issued.
In fact, after the 2011 earthquake in Japan, it took a day and a half for the tsunami warning to be lifted.
This is an advantage of our method, even if it takes time to obtain useful information.
As for the method to obtain the initial tsunami information by calculating inverse from the TIH information, we expect that the combination of the TIH estimated by our method and the acoustic wave propagation model may be able to estimate the initial tsunami area.
At present, although H. Shinagawa et al. (2013) and M. D. Zettergren et al. (2019) have been able to reproduce the TIH quite accurately, they have not yet been able to reproduce it completely.
Therefore, we expect to be able to backcalculate the region of the initial tsunami at an early stage by supplementing the simulation model with model discrepancies (J. Brynjarsdo ́ttir and A. O′Hagan (2014)), which is a statistical method that takes into account the differences between actual measurements and model outputs for estimation, and history matching (I. Vernon et al. (2014)), which can limit the likelihood region of model parameters.
 Shinagawa, T. Tsugawa, M. Matsumura, T. Iyemori, A. Saito, T. Maruyama, et al. Twodimensional simulation of ionospheric variations in the vicinity of the epicenter of the TohokuOki earthquake on 11 March 2011. Geophysical Research Letters, 40(19):5009–5013, 2013.
 D. Zettergren and J. B. Snively. Latitude and Longitude Dependence of Ionospheric TEC and Magnetic Perturbations From InfrasonicAcoustic Waves Generated by Strong Seismic Events. Geophysical Research Letters, 46(3):1132–1140, 2019.
 Brynjarsdo ́ttir and A. O′Hagan (2014). Learning about physical parameters: The importance of model discrepancy. Inverse problems, 30(11):114007, 2014.
 Vernon, M. Goldstein, and R. Bower. Galaxy formation: Bayesian history matching for the observable universe. Statistical science, pages 81–90, 2014.
[About the first minor comment.]
According to the results of Maruyama's analysis (Maruyama et al., 2011) of ionogram information on the day of 2011 off the coast of the Tohoku earthquake, the electron density peak of the ionosphere was at an altitude of 306 km in the data from Kokubunji, which is the closest to the epicenter at 440 km.
This result suggests that the assumption of a hypothetical thin ionosphere at an altitude of 300 km in our analysis is reasonable.
Note that Maruyama's analysis was conducted for the March 11, 2011 earthquake, and a similar analysis is needed to determine whether setting the ionosphere at 300 km is the most appropriate assumption for other earthquake cases.
 Maruyama, T., Tsugawa, T., Kato, H., Saito, A., Otsuka, Y., and Nishioka, M.: Ionospheric multiple stratifications and irregularities induced by the 2011 off the Pacific coast of Tohoku Earthquake, Earth, Planets and Space, 63, 65, https://doi.org/10.5047/eps.2011.06.008, 2011.
[About the second minor comment.]
Yes, thank you very much for your comment. I modified these expressions.
[About the third minor comment.]
To make it easier to understand the uncertainty, we replaced the 3D plot with a 2D plot of the 99% confidence interval as below.
[About the fourth minor comment.]
The location of the tsunami source is the same value used in Kamogawa et al. (2016), and its coordinates are calculated by referring to Maeda et al. (2011), Grilli et al. (2013), Ohta et al. (2012), and Saito et al. (2011b). The eight directions are evenly divided into North, Northwest, West, Southwest, South, Southeast, East, and Northeast from its tsunami center position.
 Kamogawa, M., Orihara, Y., Tsurudome, C., Tomida, Y., Kanaya, T., Ikeda, D., et al.: A possible spacebased tsunami early warning system using observations of the tsunami ionospheric hole, Scientific reports, 6, 37 989, https://doi.org/https://doi.org/10.1038/srep37989, 2016.
 Maeda, T., Furumura, T., Sakai, S., and Shinohara, M.: Significant tsunami observed at oceanbottom pressure gauges during the 2011 off the Pacific coast of Tohoku Earthquake, Earth, Planets and Space, 63, 53, https://doi.org/https://doi.org/10.5047/eps.2011.06.005, 2011.
 Grilli, S. T., Harris, J. C., Bakhsh, T. S. T., Masterlark, T. L., Kyriakopoulos, C., Kirby, J. T., and Shi, F.: Numerical simulation of the 2011 Tohoku tsunami based on a new transient FEM coseismic source: Comparison to farand nearfield observations, Pure and Applied Geophysics, 170, 1333–1359, https://doi.org/https://doi.org/10.1007/s000240120528y, 2013.
 Ohta, Y., Kobayashi, T., Tsushima, H., Miura, S., Hino, R., Takasu, T., et al.: Quasi realtime fault model estimation for nearfield tsunami forecasting based on RTKGPS analysis: Application to the 2011 TohokuOki earthquake (Mw 9.0), Journal of Geophysical Research:Solid Earth, 117, https://doi.org/https://doi.org/10.1029/2011JB008750, 2012.
 Saito, T., Ito, Y., Inazu, D., and Hino, R.: Tsunami source of the 2011 TohokuOki earthquake, Japan: Inversion analysis based on dispersivetsunami simulations, Geophysical Research Letters, 38, https://doi.org/https://doi.org/10.1029/2011GL049089, 2011b.
In the case that the embedded images are not shown clearly, I attached a PDF file including all the images in this reply.
Ryuichi Kanai et al.
Ryuichi Kanai et al.
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

243  62  8  313  17  2  3 
 HTML: 243
 PDF: 62
 XML: 8
 Total: 313
 Supplement: 17
 BibTeX: 2
 EndNote: 3
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1