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 non-overlapping 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 east-west 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 east-west direction, resulting in the recombination of ions and electrons, causing the TIH to spread in the east-west 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 P-waves (longitudinal tremors) of about 7 kilometers per second and S-waves (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 sea-level change, and the initial sea-level change cannot be known.

In addition to the above methods, Japan has developed a system called REal-time GEONET Analysis system for Rapid Deformation monitoring (REGARD), which analyzes GEONET data in real-time 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 real-time 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 back-calculate 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. Two-dimensional simulation of ionospheric variations in the vicinity of the epicenter of the Tohoku-Oki 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 Infrasonic-Acoustic 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 3-D plot with a 2-D 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 space-based 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 ocean-bottom 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 co-seismic source: Comparison to far-and near-field observations, Pure and Applied Geophysics, 170, 1333–1359, https://doi.org/https://doi.org/10.1007/s00024-012-0528-y, 2013.

• Ohta, Y., Kobayashi, T., Tsushima, H., Miura, S., Hino, R., Takasu, T., et al.: Quasi real-time fault model estimation for near-field tsunami forecasting based on RTK-GPS analysis: Application to the 2011 Tohoku-Oki 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 Tohoku-Oki 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.

Thank you very much for your comments.

Our answers are as below.

Also, a PDF file of this answer is attached to show clearer figures.

General comments:

[About the first general comment.]

To display more understandable and convincing pictures, we created figures in which the estimated uncertainty of TEC values in Figure 5 is projected onto 2-dimension as below.

Figure 1：Left-hand side is for the full data and right-hand side is for the sparse data using only 5% of the GEONET receivers. (a) and (b) are measured TEC data. (c) and (d) are measured TEC data (blue dots) and the fitting surface (red surface). (e) and (f) are 2D projection of the fitting surface. (g) and (h) are 2D projection of the 99% one-sided confidence interval of the fitting surface. The fitting surface is computed using the INLA-SPDE method. The black dashed line indicates the Japan Trench.

[About the second general comment.]

Tatsuhiko SAITO, the first author of Saito et al. (2011b), shared the initial tsunami height data used to depict the figure in their paper.

Using the data, we compared the TEC depression with the initial tsunami shown in the below figures.

Firstly, we created the contours of the initial tsunami. From the figures in the first row, the region where the TEC deepest reaches its minimum is almost the same with the region in which the initial tsunami shows its peak value.

Next, we show the initial tsunami and the TEC depression at three different times at the fixed latitude and longitude where the initial tsunami reaches its peak value,

From the figures in the second and third rows, it is clear that the TIH deepest location and the initial tsunami highest points are the same at the fixed latitude and longitude.

In the next figure, the time series (from 5:46:30 to 6:16:30 (UTC)) data of the TIH is shown with the initial tsunami shape at the fixed latitude and longitude.

In the figure, the time series of the observed data in the atmosphere shows a strong correlation with the initial tsunami in its shape for the first time.

Figure2: (a1) – (a3) are TEC values and the initial tsunami contours. (b1) – (b3) are the TIH shape and the initial tsunami shape at the fixed latitude where the initial tsunami reaches its peak value. (c1) – (c3) are the TIH shape and the initial tsunami shape at the fixed longitude where the initial tsunami reaches its peak value.

The initial tsunami source presented by other researchers, shown in Ohta (2012) and Takagawa (2012), also shows that the initial tsunami region overlaps the TIH region as in Figure 3. In Ohta (2011), they used the algorithm developed by Okada [1985] to compute the initial sea-surface displacement based on their fault-determination procedure. On the other hand, in Takagawa (2012), they investigated the effect of the rupture process on a tsunami source inversion. Estimated sea-surface elevation of tsunami source by their inversion method based on the assumption of finite rupture velocity of 2 km/sec is shown.

Figure3: TIH and the initial tsunami comparison. In the first row, the initial tsunami estimated by Ohta (2012) is shown. In the second row, the initial tsunami estimated by Takagawa (2012) is shown. In both cases, the TIHs with TEC values less than -4 are described.

[Reference]

Ohta (2012):

Ohta, Yusaku, et al. "Quasi real‐time fault model estimation for near‐field tsunami forecasting based on RTK‐GPS analysis: Application to the 2011 Tohoku‐Oki earthquake (Mw 9.0)." Journal of Geophysical Research: Solid Earth 117.B2 (2012).

Takagawa (2012):

Takagawa, Tomohiro, and Takashi Tomita. "Effects of rupture processes in an inverse analysis on the tsunami source of the 2011 Off the Pacific Coast of Tohoku earthquake." The Twenty-second International Offshore and Polar Engineering Conference. OnePetro, 2012.

[About the third general comment.]

We showed the averaged results of the sparse data experiments in our manuscript.

To show the fitting method can detect the surface in each sparse data experiment in our analysis, several sparse data experiment results are described as below.

If there is no data point in the area where TIH is formed, it is impossible to capture the decrease in TEC, but if there are several observation points in the TEC decrease region, the fitting method works properly.

Figure4: Three different sparse data and its fitting surface mapped onto 2 dimension.

TIH overlapping the initial tsunami:

[About the first comment related to TIH overlapping the initial tsunami]

In our manuscript, we used the initial tsunami presented in Saito 2011 to show the overlap with TIH.

The region where the wave height of the initial tsunami is higher than 2 m is shown in blue as the initial tsunami range.

In order to compare the tsunami with TIH, including the shape of the tsunami, the data was provided by Saito and the comparison is shown in Figure2 and Figure5.

Here is a comparison of them using contour lines.

In addition, the shape of the initial tsunami and the shape of the TIH were compared when fixed at the latitude and longitude where the initial tsunami height is the highest.

Moreover, the time series of TIH for fixed latitudes and longitudes are also presented.

From these results, it is inferred that the TIH has information about the initial tsunami.

The papers that the reviewer gave in the comment are not about the initial tsunami height but about the distribution of slips, so it cannot be used for comparison with TIH as Saito 2011 did.

On the other hand, Ohta (2012) and Takagawa (2012) are papers that mention the distribution of initial tsunami heights.

By referring to the figures in these two papers, we can compare the initial tsunami with TIH, and we can confirm that there is an overlap between TIH and the initial tsunami as shown in Figure3.

[About the second comment related to TIH overlapping the initial tsunami]

In order to make a detailed comparison between the TIH and the initial tsunami, the time series of the shape of the TIH and the shape of the initial tsunami are compared.

The comparison is made by fixing the latitude and longitude at the position where the initial tsunami is the highest.

Figure5: The red curves are the time series from 5:46:30 to 6:16:30 of TIH at fixed latitude and longitude where the initial tsunami estimated by Saito (2011) reaches its peak. The blue curve is the initial tsunami.

[About the third comment related to TIH overlapping the initial tsunami]

Our manuscript shows for the first time that our method can properly fit the TIH and obtain the information of the initial tsunami from the fitted data.

The Tohoku-Oki earthquake, which occurred in Japan with a dense GNSS network, is the best case study for the proper application of our method.

The Tohoku-Oki earthquake is the only case in which the tsunami-induced ionospheric changes were observed in such a dense GNSS network.

In addition, this earthquake is the only case where such a huge TIH was observed.

These facts are the reason why we applied our method to the Tohoku-Oki earthquake.

As commented by the reviewer, we consider the application of our method to other cases as a future task.

Since our method has shown its applicability to sparse GNSS networks, we would like to apply it to other cases when the scientific validity of this manuscript is recognized after this peer review.

Minor comments:

[About the first minor comment]

We will add the definition of TECu when the word appears the first time on page 3.

Also, we will delete the definition described on page 4.

[About the second minor comment]

We will add the definition of O(n³), that is, the computational cost is the cubic of the number of data points.

[About the third minor comment]

The accuracy is the closeness to the true value, the uncertainties are the dispersion of values in the spatial interpolation, and the robustness is against the absence of measurements.

[About the fourth minor comment]

We will put a space between elements and to.

[About the fifth minor comment]

We will delete the sentence.

[About the sixth minor comment]

We will delete the words.

[About the seventh minor comment]

We will delete the sentence.

[About the eighth minor comment]

According to Shinagawa (2013), a study that runs a physics-based simulation model, the initial tsunami reproduces TIH. Then, TIH is spread out over the initial tsunami.

Kamogawa (2016) uses the TEC values around the tsunami source area to derive the relationship between TIH and the initial tsunami.

The initial tsunami is shown as a yellow circle in the figure.

Figure 6: TIH expansion with the tsunami source which is expressed as the yellow circle.

[About the ninth minor comment]

For example, an article, Sato, F., Tanabe, T., Murase, H., Tominari, M., & Kawai, M. (2017). Application of a wearable GPS unit for examining interindividual distances in a herd of Thoroughbred dams and their foals. Journal of equine science, 28(1), 13-17., shows the Hubeny’s distance formula in the section Materials and Methods, the paragraph Calculation of the distance between GPS units.

DOI https://doi.org/10.1294/jes.28.13

[About the tenth minor comment]

Due to the influence of the geomagnetic field, the plasma tends to move more to the south.

In addition, the backward moving average frequency filter applied in this study does not completely exclude the high frequency component.

For these reasons, we detected a decrease in the electron density as a result of the recombination caused by the oscillation of the plasma due to the high-frequency component on the south side. This temporary decrease in electron density is separate from the TIH formation caused by the low frequency component.

[About the eleventh minor comment]

In Figure 8, TEC values less than -3 and -4 are depicted. In this context, the larger means TEC values less than -3. However, there is no explicit definition about large TEC reduction.

Based on the atmospheric physics, the background TEC value is around 20 TECu, then 5% of its value, for example, 1 TECu can be mentioned as a large TEC change.

In this case, the magnitude of earthquakes has a big impact on TEC change, so defining the large TEC change is difficult.

About this matter, we created Figure 5 shown above.

The TEC dip corresponds to the initial tsunami shape is described in the figure.

[About the twelfth minor comment]

In figure 8, panels (a1), (a2), and (a3), the TIH with TEC values less than -3 is located on the region which is almost the same with the initial tsunami region. In other words, the TIH is formed above the initial tsunami region.

[About the thirteenth minor comment]

TIH is formed by the extremely low frequency component of acoustic waves that travel from the sea surface to the ionosphere.

The recombination of plasma and electrons caused by this low-frequency component is the direct cause of TIH, and this recombination takes time, so the volume continues to increase over a long period of time.

Simulations of TIH formation using a physical model presented by Shinagawa (2013) also show that TIH forms over a similar time period.

[About the fourteenth minor comment]

Huge simulated tsunami in this context refers to tsunamis with a simulated initial tsunami height of around 5 m in Kamogawa (2016). On the other hand, smaller simulated tsunami in this context means tsunamis with simulated initial tsunami heights of 1 m or less.

[About the fifteenth minor comment]

We will add that The CI means the confidence interval, which shows the degree of uncertainty.

[About the sixteenth minor comment]

In this context, a huge tsunami means a tsunami with the same magnitude with Tohoku-Oki earthquake tsunami.

However, since there is no general definition of a huge tsunami, it is impossible to describe it in a strict and explicit manner.

By convention, tsunamis triggered by magnitude about 9 earthquakes are considered to be huge tsunamis.

[About the seventeenth minor comment]

We will delete the phrase.

[About the eighteenth minor comment]

For the visual representation, we created Figure 2 and Figure 5 above based on the initial tsunami estimated by Saito (2011).

In Figure 2, the initial tsunami contours and the fixed latitude and longitude tsunami shape are compared with the TIH. In Figure 5, the time series of TIH shape at the fixed latitude and longitude are compared with the initial tsunami.

In Figure 3 shown above, the comparisons between TIH and the initial tsunami estimated by other research groups (Ohta (2012) and Takagawa (2012)) are shown.

[About the nineteenth minor comment]

We will correct the brackets.

[About the twentieth minor comment]]

In Kamogawa (2016) Figure 3, we can see the positive correlation between the percentage of TEC depression and simulated initial tsunami height.

In our manuscript, we use the volume of TIH as a measure for TIH, but we applied our new method for the Tohoku-Oki tsunami to show our method can detect the TIH volume.

To relate our measure to a corresponding measure of the tsunami size, it is necessary to apply the method to other tsunami cases.

We would like to apply it to other cases when the scientific validity of this manuscript is recognized after this peer review.