Improvement of Near-field Tsunami Forecasting Method Using Ocean Bottom Pressure Sensor Network (S-net)

Since the installation of a dense cabled observation network around the Japan Trench (S-net) by the Japanese government that includes 150 sensors, several tsunami forecasting methods that use the data collected from the ocean floor sensors were developed. One of such methods is the tsunami forecasting method which assimilates the data without any information of earthquakes. The tsunami forecasting method based on the assimilation of the ocean-bottom pressure data near the source area was developed by Tanioka in 2018. However, the method is too simple to be used for an actual station distribution of S-net. To overcome its limitation, we developed an interpolation method to generate the appropriate data at the equally spaced positions for the assimilation from the data observed at sensors in S-net. The method was numerically tested for two large underthrust fault models, a giant earthquake (Mw8.8) and the Nemuro-oki earthquake (Mw8.0) models. Those fault models off Hokkaido in Japan are expected to be ruptured in the future. The weighted interpolation method, in which weights of data are inversely proportional to the square of the distance, showed good results for the tsunami forecast method with the data assimilation. Furthermore, results indicated that the method is applicable to the actual observed data at the S-net stations. The only limitation of the weighted interpolation method is that the computed tsunami wavelengths tend to be longer than the actual tsunamis wavelength.


Introduction
The 2011 Tohoku-oki earthquake (Mw9.0) generated a large tsunami along the Pacific coast of northern Japan (Mori et al., 2012). Although the Japan Meteorological Agency (JMA) issued a major tsunami warning along the Pacific coast of Japan immediately after the earthquake (Ozaki, 2011), the tsunami caused a catastrophic disaster with approximately 19,000 casualties. After this event, the development of a more accurate and rapid tsunami warning system became a highpriority focus area in Japan.
For this purpose, the Japanese government installed a dense cabled observation network, called the seafloor observation network for earthquakes and tsunami around the Japan Trench (S-net), in 2017. The network is operated by the National Research Institute for Earth Science and Disaster Resilience (NEID). In this network, 150 observation stations consisting of ocean-bottom pressure sensors and seismometers are connected by cables at the 30-km intervals (Uehira et al.,2012;Kanazawa, 2013). The cables are distributed on the seafloor offshore off northern Japan ( Figure   1).
Recently, various early forecasting methods for tsunami using the dense seafloor observation networks have been proposed. Tsushima et al. (2012)

developed a method called tsunami
Forecasting based on Inversion for initial sea Surface Height (tFISH) which estimates the initial sea surface deformation from tsunami waveforms observed at ocean bottom pressure gauges. Yamamoto et al. (2016) developed a tsunami forecasting method that rapidly estimates tsunami inundation based on a multi-index method and compares the observed tsunami heights with the pre-computed tsunami heights at ocean bottom sensors. Moreover, new tsunami computation methods based on the assimilation of tsunami observations without the tsunami source information were developed by Maeda et al., (2015), Tanioka (2018), and Tanioka and Gusman (2018). Gusman et al. (2016) applied the method developed by Maeda et al. (2015) to the tsunami generated by the 2012 Haida Gwaii earthquake and computed the tsunami wave field successfully by assimilating data observed at ocean bottom pressure sensors in Cascadia. However, it is difficult to apply this method with ocean bottom pressure data within the source area because the ocean surface displacement above the source area just after the earthquake cannot be observed by the ocean bottom pressure sensors. The near-field tsunami inundation forecasting method developed by Tanioka (2018) and Tanioka and Gusman (2018) uses the time derivative of the pressure waveforms observed at the ocean bottom pressure sensors near the source area. Therefore, tsunami computation by assimilating data can be performed without any tsunami source information as soon as the earthquake or tsunami generation is completed. Tanioka and Gusman (2018) tested ocean-bottom pressure sensors equally distributed at 15 min intervals or approximately 30 km aparts. In reality, the S-net ocean-bottom sensors are not installed at uniform intervals (Figure 1), particularly in the north-south direction, and have a lesser number than those tested by Tanioka and Gusman (2018). They finally concluded that it is necessary to improve their method by using the exact locations of the S-net sensors for real time tsunami inundation forecasts.
In this paper, we improve the tsunami forecasting method proposed by Tanioka (2018) and Tanioka and Gusman (2018) to be able to use the exact locations of the S-net sensors. This improved method was tested for two expected large underthrust earthquakes off the Pacific coast of Hokkaido, Japan: an expected magnitude 9 class earthquake, and the Nemuro-oki earthquake with a magnitude of 8.

Previous Method to Determine the Tsunami Height Field
Tanioka (2018) presented a tsunami simulation method that assimilates dense ocean bottom pressure data near the tsunami source regions. Tanioka and Gusman (2018) suggested a near-field tsunami inundation forecasting method based on the method presented by Tanioka (2018). Here, we briefly explain the method used to determine the tsunami height field. A governing equation for the tsunami propagation, the wave equation of the linear shallow-water approximation, is as where d is the depth of the ocean, ( , , ) is the ocean surface displacement or tsunami height where i and j are indices of points for the x-and y-directions in the computed domain, k is the index for time steps, ∆ and ∆ are the spacing intervals in the x and y directions, respectively, and ∆ is the time interval. Afterwards, we assumed that there were hypothetical pressure sensors at each point (marked in red dots in Figure 2). The spacing interval for both the x and y directions is 10 arc-min or approximately 18 km. As such, the total number of hypothetical sensors was assumed to be 275. The left-hand side of Eq.
(2) was obtained from the pressure data at each hypothetical pressure sensor. The unknown parameters were ocean surface displacement field or tsunami height filed, , , above each hypothetical pressure gauge at a particular time, k. The tsunami heights outside the hypothetical pressure gauge network were assumed to be 0.
Subsequently, we were able to compute the tsunami height field, , , above the hypothetical pressure sensor network by solving Eq.
(2) for all hypothetical pressure sensors by using the inversion technique as described by Tanioka (2018).

Interpolation Method for Data Assimilation
A problem with the method mentioned above is that in practice, the S-net sensors are not distributed at the10 arc-min intervals in the x and y direction, and the number of the S-net sensors is much less than 275 ( Figure 2). To solve this problem, we need to develop an appropriate interpolation method to obtain pressure data at the hypothetical points situated at the 10 arc-min intervals (red dots in Figure 2) from the observed pressure data at the S-net sensors (black dots in , at each hypothetical point (xi,j, yi,j) at the 10 arc-min intervals (red dots in Figure   2) was obtained by the weighted interpolation equation as follows: ( 3) where wp,i,j, is the weight factor for each observed point (p) to calculate , at each hypothetical point (i, j). Then, the weight factor was defined as follows: where is the control factor. If control factor equals one, the weights are inversely proportional to the distance between the observed points in S-net and hypothetical points for assimilation. In this study, we varied control factor from 1 to 3 (integers only) to find the best factor that provides an appropriate tsunami height distribution near a tsunami source area using the assimilation technique developed by Tanioka (2018).

４．Test Region and Fault Models
Great underthrust earthquakes repeatedly occurred in the past at the plate interface in the subduction zone where the Pacific plate subducts beneath the Kurile Trench off Hokkaido ( Figure   1). Historically, the Pacific coast of Hokkaido has suffered from large disasters due to those earthquakes and their associated tsunamis (Tanioka et al., 2007). Additionally, the paleoearthquake studies using tsunami deposits along the Pacific coast of Hokkaido revealed that much larger tsunamis repeatedly occurred in the past 6000 years compared with historical tsunamis (Nanayama et al., 2003;Hirakawa et al, 2005). The most recent large paleo-tsunami that occurred in the early 17 th century ( Figure 1) was well studied by Satake et al., (2008) and Ioki and Tanioka (2016). The fault model of the earthquake was constructed by reproducing the observed tsunami deposit distribution along the Hokkaido coast with the computed one. Ioki and Tanioka (2016) concluded that the great underthrust earthquake with a moment magnitude of 8.8 ruptured a large area of the plate interface off Hokkaido. The Headquarters for Earthquake Research Promotion, Japan, published a long-term evaluation report on the occurrence of large subduction earthquakes along the Kurile Trench (Headquarters of Earthquake Research Promotion, 2018). They reported that the probability of a giant earthquake with a magnitude exceeding Mw 8.8 along the Kurile trench in 30 years is between 7% and 40%. Therefore, mitigation of tsunami disasters due to the giant earthquake (M8.8) has become an urgent challenge in Japan, particularly in Hokkaido.
In this study, we chose the hypothetical giant earthquake (Mw8.8) along the Kurile trench off Hokkaido to test our tsunami simulation by assimilating pressure data observed by the S-net sensors near the source region. The fault length and fault width for this hypothetical earthquake ( Figure 2) were set to be 300 km and 100 km, respectively, as per the estimations of the 17 thcentury giant earthquake fault model (Satake et al., 2008;Ioki and Tanioka, 2016). A strike of 228°, a dip of 15°, and a rake of 90° were set to be the same as those in the 17 th -century giant earthquake fault model (Table 1). A slip amount was calculated to be 15 m by assuming a moment magnitude of 8.8 and a rigidity of 4 x 10 10 N/m 2 .
The Headquarters for Earthquake Research Promotion also reported that the probability of a great Nemuro-oki earthquake (Mw7.8-8.5) along the Kurile trench in 30 years is approximately 70%, being higher than that of the giant earthquake (Mw8.8). In this region, the 1973 Nemuro-oki earthquake (Mw 7.8) (Figure 1) and the 1894 Numuro-oki earthquake (Mw8.2) were previously occurred. Tanioka et al. (2007) showed that the fault size of the 1894 event (200km x 100 km) was larger than that of the 1973 event (80 km x 80 km). In this study, we chose the Nemuro-oki earthquake (Mw 8.0) to test our method. The fault length and fault width of the earthquake ( Figure   2) were set to be 100 km and 80 km, respectively. A strike of 240°, a dip of 15°, and a rake of 90° were set the same as those in the fault model of the Nemuro-oki earthquakes (Table 1). A slip amount was calculated to be 4 m by assuming a moment magnitude of 8.0 and a rigidity of 4 x 10 10 N/m 2 .

Method for Numerical Computation Tests
Two reference tsunamis, used as the original tsunamis for synthetic tests, were computed in the area off Hokkaido and Tohoku (Figure 2). The coseismic vertical deformation was calculated using the Okada's (1985) equations from the fault models of the two reference earthquakes (Table 1) and used as the initial conditions of tsunamis. The tsunamis were numerically computed by solving the linear long-wave equations by using the finite difference method with a staggered grid system (Satake, 2015). The grid size of the tsunami computation was set at 1 arc-min (approximately 1.8 km). The time step of the computation was set to be 2 s to satisfy the stability condition. The duration of the earthquake was chosen to be 90 seconds.
To test the assimilation method, the tsunami waveforms at the S-net observation points (black dots in Figure 2) were computed numerically from the reference fault models and transferred to the water-depth fluctuation, hf(x,y,t) in Eq. (1), as described in Tanioka (2018). Afterwards, the left-hand side of Eq.
(2), = ℎ +1 −2ℎ +ℎ −1 ∆ 2 , at each sensor position (xp, yp) in S-net at a time k was calculated from the water-depth fluctuation data. Time intervals, ∆ , for the giant earthquake (Mw8.8) and the Nemuro earthquake (Mw8.0) cases were set at 6 s and 30 s, respectively. Time interval, ∆ , needs to be larger for a smaller earthquake because term, ℎ +1 − 2ℎ + ℎ −1 , needs to be large enough to be observed at the S-net sensors. Then, those at equally distributed hypothetical points (xi,j, yi,j) at the 10 arc-min intervals (red dots in Figure 2), , = ℎ , +1 −2ℎ , +ℎ , −1 ∆ 2 , were obtained by the above interpolation methods. Subsequently, the tsunami height fields at equally distributed hypothetical points at a time k were computed using the assimilation method as described by Tanioka (2018). For tsunami simulation, the tsunami height field obtained at a time k from the observed depth fluctuation data at the S-net sensors replaced the tsunami height field at a time k during the tsunami numerical simulation. The replacements of the tsunami height field as data assimilations were completed at 90 s after the initiation of the earthquake as the duration of the earthquake was set at 90 s. Finally, the tsunami numerical simulation was continued until 90 min after the earthquake by solving the linear long-wave equations without data assimilation.

Results
In the first test case, a tsunami caused by the giant earthquake ( Figure 4. The maximum tsunami heights at these seven locations are also presented in Table 2. The tsunami waveforms ( Figure 4) and maximum tsunami heights (Table   2) from the reference tsunami computation were comparable to those from the assimilation of simulated observation data using the interpolation, where control factor α equaled two. The maximum tsunami heights computed from the assimilation of the data using the interpolation, where control factor α equaled one, were underestimated compared to the maximum tsunami heights computed directly from the reference fault model. Some of the tsunami waveforms from the assimilation of those data were different from those of the reference model, especially the waveforms at Hanasaki in Hokkaido. On the contrary, the maximum tsunami heights computed from the assimilation of the data using interpolation, where control factor α equaled 3, were overestimated in comparison with those computed directly from the reference fault model.
Therefore, we conclude that the weighted interpolation method, where control factor α equaled two in Eq. (4), is appropriate for our assimilation method using the ocean-bottom pressure data observed at the S-net sensors.
As for the second test case, the tsunami caused by the fault model of the Nemuro-oki earthquake (Mw 8.0) was numerically computed as a reference tsunami. The computed tsunami height fields at 90 and 1200 s after the initiation of the earthquake are shown in Figure 5a. The tsunami height fields were computed by the assimilation of simulated observation data at the S-net sensors ( Figure   2) using the interpolation method as described by Eq. (3), where control factor α equaled two ( Figure 5b). Seven computed tsunami waveforms from the reference tsunami and from the data assimilations using the interpolation are compared in Figures 6a-b. The maximum tsunami heights at these seven locations are also presented in Table 3. As the assimilation area was much larger than the source area in this case, small tsunami height fields in the southwest of the assimilation area, where no tsunami was computed as a reference tsunami, were realized for the data assimilation ( Figure 5b). Therefore, the computed tsunami waveforms from the data assimilation using the interpolation (Figure 6b) have longer period characteristics and higher maximum tsunami heights than those waveforms of the reference tsunami. These results indicate that the data assimilation area needs to be limited to an area with reasonable water depth fluctuations observed at the S-net sensors.
Furthermore, we used a small assimilation area shown in Figure 7. For sensors in this small assimilation area, the calculated assimilation data = ℎ +1 −2ℎ +ℎ −1 ∆ 2 were more than 1/20 of the maximum assimilation data of all the sensors. The computed tsunami height fields at 90 and 1200 s after the initiation of the earthquake are shown in Figure 5c. Figure 6 shows that seven computed tsunami waveforms from the reference tsunami are better fitted by the computed tsunami waveforms from the data assimilation in the small area using interpolation than in the large area.
Table 3 also shows that the maximum tsunami heights of the reference tsunami better agree with those computed from the data assimilation in the small area than in the large area.

Discussion
The case study for the Nemuro-oki earthquake (Mw 8.0) showed that the reference tsunami better agreed with the computed tsunami using the small assimilation area limited to the tsunami source area. It is better to check this effect for the giant earthquake (Mw 8.8) case. We recomputed the tsunami height field using the largest assimilation area shown in Figure 8a by the assimilation of simulated observation data at the S-net sensors using the interpolation method as described by Eq.
(3), where control factor α equaled two. The computed result was compared with the previous result using the assimilation area shown in Figure 2 for the giant earthquake (Mw 8.8) case ( Figure  8b) is almost the same as the result shown in Figure 3 for the control factor, α equaled two. Figure   8c shows that seven computed tsunami waveforms using the largest assimilation area are almost same as those computed waveforms using the assimilation area shown in Figure 2. Those results indicated that the effect of limiting a size of the assimilation area is negligible to be consider for the giant earthquake (Mw 8.8) case and only significant for the Nemuro-oki earthquake (Mw 8.0) case which has a smaller tsunami source area than the giant earthquake case.
To use our improved assimilation method in practice, non-tsunami components such as seabottom acceleration changes due to a fault motion of an earthquake should be removed from the original ocean bottom pressure data. Saito and Tsushima (2016) showed that the first ~90 s of ocean bottom pressure data is highly affected by acoustic waves and seismic waves generated by the earthquake. Mizutani et al. (2020) recently developed an early tsunami detection method for near-fault ocean bottom pressure data by comparing the seismic data and pressure data at the same site to eliminate non-tsunami components. This tsunami detection method should be applied to original pressure data before our improved assimilation method is applied to the data. Another difficulty for practical use of our method is that the end of the earthquake rupture is unknown.
Although replacements of the tsunami height field as data assimilations were completed at 90 s in this study, those replacements should continue more than the duration of the earthquake in practice.
We may need to decide the ending time of those replacements using the seismic records at the same sites as the pressure observations. Finally, after pressure sensors are suffered by the strong motion of a large earthquake, mechanical malfunction of the sensors may occur as suggested by Kubota et al. (2018). Those malfunctioned sensors need to be identified as soon as the strong motion ends and the data at those sensors should not be used in our data assimilation method.

Conclusion
We developed a tsunami forecasting method using the assimilation of data obtained by interpolation of observed data at the S-net sensors near the tsunami source area. The interpolation method, as described in Eq. (3), where control factor α equaled two, performed acceptable for the tsunami forecasting method for the giant earthquake (Mw 8.8) case. However, for the Nemuro-oki earthquake (Mw 8.0) case, the wavelengths of tsunami waveforms computed using the assimilation of data were slightly longer and larger than those of the reference tsunami. This indicates the limitation of data assimilation using the suggested interpolation method because the assimilation data at the 176 hypothetical points (red dots in Figure 7) were calculated from the interpolation of data at only 12 S-net sensors for the Nemuro-oki earthquake case. Then, the assimilation area should be limited to an area where the reliable water depth fluctuations are observed at the S-net sensors. Moreover, tsunami height fields from a smaller earthquake are more difficult to be computed using our data assimilation method.       for two assimilation tsunamis using a large assimilation area (black lines, the same as Figure   4) and the largest assimilation area ( red lines) for the giant earthquake using the control factor α=2.  Table 2. Comparison of the maximum tsunami heights at seven locations computed from the reference tsunami of the giant earthquake (Mw8.8) and the assimilations using the interpolation where control factors α ranges from 1 to 3 (see Figure 4). The ratios of the maximum tsunami heights are indicated in parentheses.  Table 3. Comparison of the maximum tsunami heights at seven locations computed from the reference tsunami of the Nemuro-oki earthquake (Mw8.0) and the assimilations using the large and small assimilation areas (for their location refer to Figure 6). The ratios of the maximum tsunami heights are indicated in parentheses.