Hazard and Risk-Based Tsunami Early Warning Algorithms for Ocean Bottom Sensor S-Net System in Tohoku, Japan, Using Sequential Multiple Linear Regression

: This study presents robust algorithms for tsunami early warning using synthetic tsunami wave data at ocean bottom sensor (OBS) arrays with sequential multiple linear regression. The study focuses on the Tohoku region of Japan, where an S-net OBS system (150 pressure sensors) has been deployed. To calibrate the tsunami early warning system using realistic tsunami wave proﬁles at the S-net stations, 4000 stochastic tsunami simulations are employed. Forecasting models are built using multiple linear regression together with sequential feature selection based on Akaike Information Criterion and knee-point method to identify sensors that improve the accuracy most signiﬁcantly. The study considers tsunami wave amplitude at a nearshore location and regional tsunami loss for buildings to develop hazard-based and risk-based tsunami warning algorithms. The models identify an optimal conﬁguration of OBS stations and waiting time for issuing tsunami warnings. The model performance is compared against a base model, which only uses the earthquake magnitude and epicenter location. The result indicates that estimating the tsunami amplitude and loss via S-net improves accuracy. For the hazard-based forecasting, adding six sensors from the S-net improves the accuracy of the estimation most signiﬁcantly with an optimal waiting time of 3 min. For the risk-based forecasting, a longer waiting time between 5 and 10 min is suitable. at multiple


Introduction
A tsunami can cause great destruction to a coastal region by damaging buildings and infrastructure and by depriving people's lives. Taking the 2011 Tohoku tsunami in Japan as an example, it was generated by a moment magnitude (M) 9.0 earthquake and recorded a high run-up of 40+ m in Miyako [1]. The impact of the 2011 Tohoku earthquake and tsunami on coastal communities was tremendous and even caused a serious accident at Fukushima Daiichi Nuclear Power Plant [2]. The total death tolls were over 19,700 people, and the damage costs exceeded USD 170 billion [3,4]. Although tsunamis cannot be prevented, their severe impact can be mitigated through enhanced tsunami preparedness, timely warning, and effective response [5,6]. Therefore, the tsunami early warning system is an effective tool to reduce victims [7][8][9].
The Japan Meteorological Agency (JMA) is responsible for forecasting and issuing warnings that relate to earthquake and tsunami threats in Japan. JMA originally used an automatic monitoring method to estimate the earthquake magnitude and epicenter location and to forecast the impact of earthquakes and tsunamis. During the 2011 Tohoku tsunami, JMA was unable to forecast the tsunami triggered by this large earthquake of M9.0 using the original JMA method [10]. The inaccuracy was due to the underestimation of the earthquake magnitude caused by saturation of seismic signals, and as a result, the tsunami waves were initially underpredicted, and the tsunami warning update was too late [11].
To improve the accuracy of tsunami forecasting and early warning, JMA proposed to use more ocean bottom sensors (OBS) and increase the number of off-shore sensors in the ocean. The S-net project started in 2011 [12,13], involving a large network of submarine seismic and tsunami sensors around Japan Trench. The S-net has 150 real-time monitoring pressure sensors (tsunami meters) in total, which cover an area of 1000 × 300 km 2 from off-Hokkaido to off-Kanto [14]. The accuracy of tsunami forecasting depends on how the S-net sensors are geographically located [5,6].
In calibrating a robust tsunami early warning system, it is important that the system is tested to work well for a variety of possible earthquake ruptures. However, there is no abundant historical tsunami data that can be directly used for developing such a warning system, because no major tsunami events occurred since the 2011 Tohoku tsunami. As it is uncertain how the future tsunami events may unfold, it is sensible to use synthetic tsunami wave data to calibrate a tsunami early warning algorithm based on conventional earthquake source information (e.g., JMA's magnitude estimate and epicenter location) as well as S-net wave data. For this purpose, adopting stochastic tsunami simulation models that are developed for the Tohoku region [15,16] is a suitable approach. The stochastic tsunami simulations facilitate the consideration of realistic heterogeneous earthquake slips, which are seismologically and statistically compatible with inferred earthquake slips via source inversion analyses (e.g., [17]). Once a large set of synthetic earthquake source models is simulated, wave profiles at off-shore and some coastal locations can be generated by evaluating physics-based governing equations, which capture both spatial and temporal resolutions. Subsequently, coastal tsunami inundation can be simulated, and the corresponding regional tsunami loss for buildings can be calculated by using suitable exposure and fragility models [18]. This numerical set-up is particularly useful because it facilitates the development of tsunami early warning algorithms for not only conventional hazard metrics (the probability of occurrence of the tsunami event, which is quantified by wave amplitude at a shallow coastal location near a town/city) but also risk metrics (the product of hazard, exposure, and vulnerability, which is quantified by the total tsunami loss of damaged buildings in coastal communities).
This study presents new effective and robust algorithms for tsunami early warning using synthetic tsunami wave data at S-net OBS arrays with sequential multiple linear regression (MLR). The calibration data are generated through stochastic tsunami simulations that have been carried out for the Tohoku region of Japan by considering the 2011 Tohoku tsunami hazard and risk data [15,16,18]. Two types of forecasting models are developed: the hazard-based warning system uses tsunami wave amplitude at a nearshore location, whereas the risk-based warning system adopts regional tsunami loss of buildings along the coast. For both hazard-and risk-based warning algorithms, an optimal configuration of OBS stations and waiting time for issuing tsunami warnings are identified by implementing sequential MLR together with Akaike Information Criterion (AIC) [19] and knee-point method [20]. The model performance is compared against a base model, which only uses the earthquake magnitude and epicenter location. The identification of a more compact OBS network system for tsunami early warning is motivated because the installation and maintenance of the off-shore sensors are expensive [21], and a smaller network of stations can be more easily deployed [22]. In such a system, the waiting time should be relatively short, since a trade-off exists between data acquisition time and evacuation time [23,24]. One of the important conclusions that this paper aims to draw is to demonstrate that the forecasting performance is improved by using tsunami data from the S-net, in comparison with the conventional forecasting method based on earthquake information alone. The novelty of this work includes: (i) a comprehensive set of stochastic tsunami simulation results is adopted for developing tsunami early warning systems for the S-net; (ii) risk-based tsunami warning approaches are investigated; and (iii) optimal networks for different forecasting metrics are obtained through extensive statistical analyses.

Tsunami Simulation Data
Hypothetical data, obtained from numerical tsunami simulations, are commonly used in tsunami scenarios [25][26][27][28][29]. Because historical events of earthquakes and tsunamis are lacking to perform statistical analyses, 4000 cases of simulated stochastic earthquakes are used in this study, although there is a tsunami event during the 2016 Fukushima earthquake [30]. The processes of how the data used in this study are generated are shown in Figure 1. The earthquake data, including magnitude, epicenter, and earthquake slip, are first simulated from the seismicity, fault, and scaling models [16]. Then, by solving nonlinear shallow water equations for given initial dislocation profiles of water, wave amplitude data are generated at selected recording locations (e.g., S-net stations and nearshore coastal locations). After simulating the inundation heights at building locations, the tsunami loss is evaluated by integrating hazard, exposure, and vulnerability models [18]. In the context of this study, the maximum on-shore tsunami amplitude is used as a response variable for the hazard-based early warning system, whereas the portfolio tsunami loss from buildings is adopted as a response variable for the risk-based system. Explanatory variables of the early warning system are the earthquake source information (i.e., magnitude and epicenter location) and the wave amplitude at the off-shore S-net sensor locations. networks for different forecasting metrics are obtained through extensive statistical analyses.

Tsunami Simulation Data
Hypothetical data, obtained from numerical tsunami simulations, are commonly used in tsunami scenarios [25][26][27][28][29]. Because historical events of earthquakes and tsunamis are lacking to perform statistical analyses, 4000 cases of simulated stochastic earthquakes are used in this study, although there is a tsunami event during the 2016 Fukushima earthquake [30]. The processes of how the data used in this study are generated are shown in Figure 1. The earthquake data, including magnitude, epicenter, and earthquake slip, are first simulated from the seismicity, fault, and scaling models [16]. Then, by solving nonlinear shallow water equations for given initial dislocation profiles of water, wave amplitude data are generated at selected recording locations (e.g., S-net stations and nearshore coastal locations). After simulating the inundation heights at building locations, the tsunami loss is evaluated by integrating hazard, exposure, and vulnerability models [18]. In the context of this study, the maximum on-shore tsunami amplitude is used as a response variable for the hazard-based early warning system, whereas the portfolio tsunami loss from buildings is adopted as a response variable for the risk-based system. Explanatory variables of the early warning system are the earthquake source information (i.e., magnitude and epicenter location) and the wave amplitude at the off-shore S-net sensor locations.

Stochastic Earthquakes
The stochastic tsunami simulations can be used to generate earthquake rupture models for future events [16]. The seismicity model (e.g., a statistical earthquake occurrencemagnitude model based on earthquake catalog data in the Tohoku region) characterizes the distribution of earthquakes. Figure 2 shows the two interested locations, Iwanuma and Onagawa, with their building asset distributions (see Figure 2a) and the 99 off-shore sensors that measure the wave amplitude ( Figure 2b). The regional fault source model represents the 2011 Tohoku earthquake fault plane, covering an area spanning 650 km long and 250 km wide [31]. The earthquake source zone is discretized into 10 × 10 km 2 sub-fault areas (see Figure 3a) to achieve heterogeneous slips. Eight source parameters, including fault length, fault width, mean slip, maximum slip, Box-Cox power parameter, correlation length along dip, correlation length along strike, and Hurst number, are used to determine the size, position, and slip distribution of the earthquakes [16]. The whole range of the

Stochastic Earthquakes
The stochastic tsunami simulations can be used to generate earthquake rupture models for future events [16]. The seismicity model (e.g., a statistical earthquake occurrencemagnitude model based on earthquake catalog data in the Tohoku region) characterizes the distribution of earthquakes. Figure 2 shows the two interested locations, Iwanuma and Onagawa, with their building asset distributions (see Figure 2a) and the 99 off-shore sensors that measure the wave amplitude ( Figure 2b). The regional fault source model represents the 2011 Tohoku earthquake fault plane, covering an area spanning 650 km long and 250 km wide [31]. The earthquake source zone is discretized into 10 × 10 km 2 sub-fault areas (see Figure 3a) to achieve heterogeneous slips. Eight source parameters, including fault length, fault width, mean slip, maximum slip, Box-Cox power parameter, correlation length along dip, correlation length along strike, and Hurst number, are used to determine the size, position, and slip distribution of the earthquakes [16]. The whole range of the simulated earthquake magnitudes is from M7.5 to M9.1, with 0.2 bin width, the 4000 cases of earthquake ruptures are evenly split into eight classes (the first class contains 500 cases of earthquake within M7.5-M7.7, while the eighth class contains 500 cases of earthquake within M8.9-M9.1). The earthquake magnitude is used as a fundamental explanatory variable to forecast major tsunami wave arrivals since both maximum on-shore tsunami amplitude and tsunami loss have an increasing trend with it (see Figures 4 and 5). However, the variability of tsunami amplitude and loss is still large for a specific magnitude range. Additional information is necessary to achieve a robust forecasting. simulated earthquake magnitudes is from M7.5 to M9.1, with 0.2 bin width, the 4000 cases of earthquake ruptures are evenly split into eight classes (the first class contains 500 cases of earthquake within M7.5-M7.7, while the eighth class contains 500 cases of earthquake within M8.9-M9.1). The earthquake magnitude is used as a fundamental explanatory variable to forecast major tsunami wave arrivals since both maximum on-shore tsunami amplitude and tsunami loss have an increasing trend with it (see Figures 4 and 5). However, the variability of tsunami amplitude and loss is still large for a specific magnitude range. Additional information is necessary to achieve a robust forecasting.    There are three off-shore sensors, which are A, B and C representing near, middle, and far locations from the coastline. A greater earthquake magnitude tends to generate a larger earthquake rupture in size and higher tsunami waves.  . Earthquake rupture and tsunami wave plots for three individual cases of (a) M7.6, (b) M8.2, and (c) M9.1. The yellow triangle represents the midpoint of the earthquake rupture. There are three off-shore sensors, which are A, B and C representing near, middle, and far locations from the coastline. A greater earthquake magnitude tends to generate a larger earthquake rupture in size and higher tsunami waves.

Wave Amplitude
The wave amplitude data can be used as both the response variable, i.e., maximum amplitude at the interested locations, and the explanatory variable, i.e., wave amplitude at off-shore S-net locations ( Figure 2). A well-tested TSUNAMI computer code [32] is used to generate the tsunami wave amplitude by solving nonlinear shallow water equations, which provide both spatial and temporal resolutions. The initial water surface elevation was estimated using formulae by Okada [33] and Tanioka and Satake [34]; the latter is considered to account for the effects of horizontal deformation of the seafloor in the accretionary prism of the subduction areas. A complete dataset of bathymetry/elevation, coastal/riverside structures, and surface roughness, provided by the Miyagi Prefectural Government, is organized into nesting grids of four resolutions (i.e., 1350-m, 450-m, 150m and 50-m domains). In the tsunami simulation, the coastal/riverside structures are represented by a vertical wall at one or two sides of the computational cells. To evaluate the volume of water that overpasses these walls, Honma's weir formulae are employed [35]. The bottom friction is evaluated using Manning's formula following the Japan Society of Civil Engineers standard [35]. The simulation of the wave amplitude starts when an earthquake rupture is initiated, and the duration is set to 120 min [18]. Figure 3 shows three different earthquake scenarios and waveforms at three selected locations. The farther the off-shore sensors are from the coastline, the earlier they can detect the wave arrivals. The off-shore sensors contain useful information for calibrating the warning algorithm. Figure 4 shows the plot of the maximum on-shore tsunami amplitude against earthquake magnitude in Iwanuma and Onagawa. The average tsunami amplitude increases as the earthquake magnitude increases. The high magnitude earthquakes determine the tsunami, and the earthquake information can be used as a fundamental

Wave Amplitude
The wave amplitude data can be used as both the response variable, i.e., maximum amplitude at the interested locations, and the explanatory variable, i.e., wave amplitude at off-shore S-net locations ( Figure 2). A well-tested TSUNAMI computer code [32] is used to generate the tsunami wave amplitude by solving nonlinear shallow water equations, which provide both spatial and temporal resolutions. The initial water surface elevation was estimated using formulae by Okada [33] and Tanioka and Satake [34]; the latter is considered to account for the effects of horizontal deformation of the seafloor in the accretionary prism of the subduction areas. A complete dataset of bathymetry/elevation, coastal/riverside structures, and surface roughness, provided by the Miyagi Prefectural Government, is organized into nesting grids of four resolutions (i.e., 1350-m, 450-m, 150-m and 50-m domains). In the tsunami simulation, the coastal/riverside structures are represented by a vertical wall at one or two sides of the computational cells. To evaluate the volume of water that overpasses these walls, Honma's weir formulae are employed [35]. The bottom friction is evaluated using Manning's formula following the Japan Society of Civil Engineers standard [35]. The simulation of the wave amplitude starts when an earthquake rupture is initiated, and the duration is set to 120 min [18]. Figure 3 shows three different earthquake scenarios and waveforms at three selected locations. The farther the off-shore sensors are from the coastline, the earlier they can detect the wave arrivals. The off-shore sensors contain useful information for calibrating the warning algorithm. Figure 4 shows the plot of the maximum on-shore tsunami amplitude against earthquake magnitude in Iwanuma and Onagawa. The average tsunami amplitude increases as the earthquake magnitude increases. The high magnitude earthquakes determine the tsunami, and the earthquake information can be used as a fundamental factor for the tsunami early warning. However, the variation of tsunami amplitude is significant among different earthquake magnitudes, and the earthquake magnitude may not be the most suitable feature to forecast the tsunami amplitude. Wave amplitude collected by the off-shore sensors may be more suitable for forecasting the tsunami amplitude for tsunami early warning purposes.

Tsunami Loss
After the wave data are generated, inundation depths at building locations ( Figure 2) are estimated, and tsunami damage of these buildings is evaluated by applying the tsunami fragility functions. More specifically, the tsunami fragility model [36] is implemented to evaluate the probability of different severities of tsunami damage as a function of inundation depth. Using Monte Carlo simulations, tsunami damage states are determined, and then, tsunami loss is evaluated by considering the following damage ratios: minor (0.03-0.1), moderate (0.1-0.3), extensive (0.3-0.5), complete (0.5-1.0), and collapse (1.0). For a damage state with lower and upper limits (e.g., the moderate damage state is specified with damage ratios between 0.1 and 0.3), a damage ratio is sampled from a uniform distribution. The exposure model includes the portfolio of 6152 and 1760 building assets near the coastline in Iwanuma and Onagawa, respectively. The building dataset is based on the post-2011-Tohoku tsunami damage data compiled by the Ministry of Land Infrastructure and Transportation. The building data include information of the building location, structural material, and number of stories, which can be used to estimate the cost of buildings. Finally, tsunami loss is calculated by sampling a value of the total replacement cost for a building from the lognormal distribution and then by multiplying it by the final damage ratio [18]. Figure 5 shows the tsunami loss against earthquake magnitude in Iwanuma and Onagawa. Since Iwanuma has about three times more building assets than Onagawa, the building loss in Iwanuma tends to be higher than the loss in Onagawa. The tsunami loss starts to increase when the earthquake magnitude is over M7.9 and increases significantly when the magnitude exceeds M8.5. It is noteworthy that the variability of the tsunami loss is large, even larger than the tsunami amplitude. This means that a tsunami triggered by a large magnitude earthquake has high potential to cause a large loss, but it is difficult to predict the exact loss value due to its large uncertainty by using earthquake magnitude alone. Thus, it is important to investigate whether the use of off-shore sensors can help forecast the tsunami loss more accurately and issue a more accurate tsunami early warning based on this risk metric.
The tsunami loss is a rare event even when the earthquake magnitude is relatively large, and the positions of large earthquake slips with respect to the target city/town location have significant influence on the tsunami loss. In other words, the directivity of tsunami waves is an important factor in forecasting potential tsunami consequences. To visualize this effect clearly, Figure 6 plots the earthquake magnitude and the earthquake latitude by distinguishing tsunami loss cases and no tsunami loss cases. The plots show a separation of with-loss and without-loss cases in a fan-shaped manner with its origin corresponding to the target location of the town. The fan-shaped spread of the non-zero tsunami loss data for Iwanuma is wider than that for Onagawa. Since the logarithmic transformation is applied to the loss data to reduce the skewness in the loss data, nonzero tsunami loss data are employed in the MLR analyses (see Section 2.2). It is noted that these non-zero tsunami loss data can be distinguished effectively by implementing logistic regression using earthquake magnitude and latitude differences with respect to the target location. nami loss data for Iwanuma is wider than that for Onagawa. Since the logarithmic formation is applied to the loss data to reduce the skewness in the loss data, non tsunami loss data are employed in the MLR analyses (see Section 2.2). It is noted that non-zero tsunami loss data can be distinguished effectively by implementing logis gression using earthquake magnitude and latitude differences with respect to the location.

Multiple Linear Regression with AIC Forward Selection and Knee-Point
MLR is a linear regression model that uses two or more explanatory variables to predict the response variable. The response variables in this study are the maximum on-shore tsunami amplitude and the tsunami loss (see Figures 4 and 5). The explanatory variables are earthquake magnitude (x mag ), earthquake latitude (x lat ), earthquake longitude (x long ), and the maximum wave amplitude collected at the i off-shore sensors during a specific waiting time, where i = 1, 2, . . . , 99. The base (M1) and full (M2) models can be written as: Since the distribution of the data is assumed to be normal in MLR, logarithmic transformation is applied to the response variables to reduce the skewness of the data. Non-zero tsunami loss data are used for calibrating the risk-based tsunami early warning model ( Figure 6).
Since there are 102 explanatory variables (3 earthquake parameters plus 99 off-shore sensor parameters) in the full model, the variable reduction is needed to remove redundant variables to avoid collinearity in MLR and to identify the most informative off-shore sensors. After the MLR model is fitted, variable selection using AIC and knee-point are applied to select the most informative variables. In this study, the earthquake information is kept in the model, since using the earthquake information is a standard approach to issue early warning [37,38], and this information is usually available a few minutes after an earthquake. The variable selection is applied to the 99 off-shore sensors, so that a smaller number of informative off-shore sensors can be defined. AIC is an estimator of prediction error and accounts for the trade-off between the goodness-of-fit and the simplicity of the model. The lower the AIC value is, the better the model is. AIC uses the likelihood ratio to estimate the goodness-of-fit of each model and gives a penalty when more explanatory variables are used [19]. Denoting the number of parameters in the model by k, the total number of observations by n, and the maximum likelihood estimator (MLE) byθ, the formula of AIC can be written as: θ is calculated as by maximizing the loglikelihood L: The forward selection is a method of stepwise regression that starts with the base model and adds a variable that improves the most for the model one by one. The AIC value of the base model usually is large when the forward selection procedure is initiated, and as more variables are added to the model, the AIC value will decrease rapidly, reach a minimum value, and then increase when enough parameters are added. The knee-point method has been widely used in the optimal selection problem [20,39,40]. This method is useful for identifying an effective early warning system with a minimum number of sensors. The model equation for the minimum AIC can be written as: Whereas the model equation for the knee-point based model can be written as: where 0 < k2 < k1 < 99. Figure 7 shows an example of knee-point and minimum points. The trend of the AIC score decreases significantly from the base model (M1) to the model based on the knee-point method (M4), slowly decreases to the minimum point (M3), and increases slightly to the end (M2). The corresponding numbers of sensors that are picked by the minimum point and the knee-point are used to fit the model. The knee-point method with a smaller number of off-shore sensors provides a relatively good prediction, whereas the minimum method provides the optimal selection based on AIC scores. The model performances of M1, M2, M3 and M4 are compared for both hazard-based and risk-based analyses for Iwanuma and Onagawa. The most informative number of off-shore sensors and optimal waiting time are identified based on the model performances.

Comparison of Model Performances
The process of how to fit and compare the models is shown in Figure 8

Comparison of Model Performances
The process of how to fit and compare the models is shown in Figure 8. In Section 2.2.1, for a given waiting time (e.g., for 5 min waiting time, x i represents the maximum wave amplitude up to 5 min at the ith station), the process of developing M1, M2, M3 and M4 has been presented. The AIC score plot determines the number of selected sensors (k1 and k2) for M3 and M4 (Figure 7). To determine the suitable number of recording stations for all cases of waiting time, k1 and k2 are calculated by averaging the number of minimum points and knee-points of sensors, respectively. Subsequently, the waiting time is identified by comparing the mean squared error (MSE) of models with the same number of sensors but for different waiting times (note: the model performances for the optimized S-net warning system become stable after the waiting time of 5 min; see the Section 3). the mean of the sum of the square of the difference between the actual in the data and the estimated value ̂ for each observation, and is given by: Final models (fit based on CV-10) firstly compare the MSE with themselves under different waiting times. Once the waiting time is determined, the sensors can be found by tracing back the AIC selected variables under that specific waiting time.

Results
The performances of four models M1, M2, M3 and M4 are evaluated and compared in this section. The results are presented for the hazard-based and risk-based early warning models and for the two locations Iwanuma and Onagawa ( Figure 2). The output variable in the hazard-based early warning model is the maximum on-shore tsunami amplitude (Figure 4), while that in the risk-based early warning model is the tsunami loss (Figures 5 and 6). The objective of the developed early warning systems is to identify an optimal configuration of S-net recording stations that achieve a relatively short waiting time and small MSE (Figure 8). Figure 9 presents the results of the hazard-based early warning model for Iwanuma. The AIC score plot (Figure 9a) identifies an optimal configuration of recording sensors for different criteria (i.e., M3 and M4), noting that each point on the plot corresponds to a set of selected off-shore sensors. The AIC score (Figure 9a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and then decreases slowly to the model with the minimum point (around 25 to 42 sensors). The average minimum points and knee-points are 35 and 6 stations, and these correspond to 1 and 2 in Equations (5) and (6).

Maximum On-Shore Amplitude for Iwanuma
The MSE-waiting time plot shows a trade-off between the model performance against waiting time, and Figure 9b shows the MSE of M3 (base plus 35 sensors information) and M4 (base plus 6 sensors information). A time that achieves a stable forecasting performance is chosen as an optimal waiting time. By observing that MSE decreases significantly after 3 min and becomes relatively flat after that, the suggested waiting time is 3 min for both models. For M3 and M4, different combinations of selected off-shore sensors are obtained by AIC forward selection under different waiting times. For robust estimation of the optimal network configuration, a 10-fold cross-validation (CV-10) technique is applied to each model. CV-10 evenly splits the whole data into 10 folds that separate into training data and test data with a 9:1 proportion. Then, each model is trained ten times with ten different splits of training and test data. The feature of the final model is calculated by averaging the ten split models, and the evaluation of model performance is based on the final model. The mean square error (MSE), also regarded as the misfit obtained by the CV regression model, is used as the metric to evaluate how well the final model fits the data. The MSE is the mean of the sum of the square of the difference between the actual y in the data and the estimated valueŷ for each observation, and is given by: Final models (fit based on CV-10) firstly compare the MSE with themselves under different waiting times. Once the waiting time is determined, the sensors can be found by tracing back the AIC selected variables under that specific waiting time.

Results
The performances of four models M1, M2, M3 and M4 are evaluated and compared in this section. The results are presented for the hazard-based and risk-based early warning models and for the two locations Iwanuma and Onagawa (Figure 2). The output variable in the hazard-based early warning model is the maximum on-shore tsunami amplitude (Figure 4), while that in the risk-based early warning model is the tsunami loss ( Figures 5 and 6). The objective of the developed early warning systems is to identify an optimal configuration of S-net recording stations that achieve a relatively short waiting time and small MSE (Figure 8).

Hazard-Based Results
3.1.1. Maximum On-Shore Amplitude for Iwanuma Figure 9 presents the results of the hazard-based early warning model for Iwanuma. The AIC score plot (Figure 9a) identifies an optimal configuration of recording sensors for different criteria (i.e., M3 and M4), noting that each point on the plot corresponds to a set of selected off-shore sensors. The AIC score (Figure 9a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and then decreases slowly to the model with the minimum point (around 25 to 42 sensors). The average minimum points and knee-points are 35 and 6 stations, and these correspond to k1 and k2 in Equations (5) and (6). The sensor maps for the optimal system configuration based on the AIC score plot and the MSE-waiting time plot are shown in Figure 9c,d for M3 and M4, respectively. The selected sensors for M3 cover the whole region of the S-net sensors comprehensively, whereas the selected sensors for M4 tend to include one sensor close to Iwanuma and five sensors far from the coastline and spanning from north to south. The identified configuration of the sensors for M3 attempts to utilize all available real-time wave information, while the sensor configuration for M4 attempts to cover areas near the trench where large earthquake slip asperities are located and to ensure the monitoring and detection of significant local or near-shore events by having a close-by sensor. Figure 9. Results of the hazard-based early warning model for Iwanuma: (a) AIC score plot, (b) MSE-waiting time plot, (c) optimal sensor map for M3, and (d) optimal sensor map for M4. The pink star represents the target location (Iwanuma). The red and empty circles represent the selected and unselected sensors, respectively. Figure 10 shows how the spatial locations of the selected sensors for M4 change when the waiting time changes in forecasting the maximum on-shore tsunami amplitude in Iwanuma. Figure 10a shows the selected six sensors when the waiting time is 5 min. The sensors tend to be far from the coastline. As the waiting time increases to 10 min ( Figure  10b) and 20 min (Figure 10c), the selected sensors tend to be closer to the coastline, and more sensors are selected toward the direction of Iwanuma with respect to the center of the rupture area. From the MSE plot (Figure 10b), the MSE values of M3 and M4 have similar trends and are close to each other. The performance of M3 and M4 in terms of MSE is not significantly different. Table 1 summarizes the mean values of the MSE with regard to forecasting the maximum on-shore tsunami amplitude in Iwanuma. M2, M3 and M4 have significantly lower MSE values compared to M1, clearly demonstrating that a small number of pre-selected S-net sensors, which measure the wave amplitude during an early period after an earthquake, provide useful information for forecasting the tsunami hazard in Iwanuma. The MSE-waiting time plot shows a trade-off between the model performance against waiting time, and Figure 9b shows the MSE of M3 (base plus 35 sensors information) and M4 (base plus 6 sensors information). A time that achieves a stable forecasting performance is chosen as an optimal waiting time. By observing that MSE decreases significantly after 3 min and becomes relatively flat after that, the suggested waiting time is 3 min for both models.
The sensor maps for the optimal system configuration based on the AIC score plot and the MSE-waiting time plot are shown in Figure 9c,d for M3 and M4, respectively. The selected sensors for M3 cover the whole region of the S-net sensors comprehensively, whereas the selected sensors for M4 tend to include one sensor close to Iwanuma and five sensors far from the coastline and spanning from north to south. The identified configuration of the sensors for M3 attempts to utilize all available real-time wave information, while the sensor configuration for M4 attempts to cover areas near the trench where large earthquake slip asperities are located and to ensure the monitoring and detection of significant local or near-shore events by having a close-by sensor. Figure 10 shows how the spatial locations of the selected sensors for M4 change when the waiting time changes in forecasting the maximum on-shore tsunami amplitude in Iwanuma. Figure 10a shows the selected six sensors when the waiting time is 5 min. The sensors tend to be far from the coastline. As the waiting time increases to 10 min (Figure 10b) and 20 min (Figure 10c), the selected sensors tend to be closer to the coastline, and more sensors are selected toward the direction of Iwanuma with respect to the center of the rupture area. From the MSE plot (Figure 10b), the MSE values of M3 and M4 have similar trends and are close to each other. The performance of M3 and M4 in terms of MSE is not significantly different. Table 1 summarizes the mean values of the MSE with regard to forecasting the maximum on-shore tsunami amplitude in Iwanuma. M2, M3 and M4 have significantly lower MSE values compared to M1, clearly demonstrating that a small number of pre-selected S-net sensors, which measure the wave amplitude during an early period after an earthquake, provide useful information for forecasting the tsunami hazard in Iwanuma.  The general observations of the results for Onagawa are similar to those for Iwanuma. The AIC score (Figure 11a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors). The AIC score decreases slowly to the minimum point (around 30 to 45 sensors). Figure 11b shows the MSE for M3 (base plus 35 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 3 min for both models (same as Section 3.1.1). Given the number of sensors added to the model and the waiting time, the selected sensors for M3 and M4 are shown in Figure  11c,d. The selected sensors for M3 cover the entire region of the S-net, whereas the selected sensors for M4 are located far from the coastline only with one sensor close to Onagawa. Figure 12 shows how the spatial locations of the selected sensors for M4 change when the waiting time changes in forecasting the maximum on-shore tsunami amplitude in Onagawa. With the increase in waiting time, the locations of off-shore stations tend to be  The general observations of the results for Onagawa are similar to those for Iwanuma. The AIC score (Figure 11a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors). The AIC score decreases slowly to the minimum point (around 30 to 45 sensors). Figure 11b shows the MSE for M3 (base plus 35 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 3 min for both models (same as Section 3.1.1). Given the number of sensors added to the model and the waiting time, the selected sensors for M3 and M4 are shown in Figure 11c,d.
The selected sensors for M3 cover the entire region of the S-net, whereas the selected sensors for M4 are located far from the coastline only with one sensor close to Onagawa.

Tsunami Loss for Iwanuma
The AIC score (Figure 13a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and decreases slowly to the minimum point (around 20 to 42 sensors). Figure 13b shows the MSE of M3 (base plus 33 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 5 min for both models, similar to the optimal waiting time of 3 min using the hazard-based algorithm. Given the number of sensors added to the model and the waiting time, the Figure 11. Results of the hazard-based early warning model for Onagawa: (a) AIC score plot, (b) MSE-waiting time plot, (c) optimal sensor map for M3, and (d) optimal sensor map for M4. The pink star represents the target location (Onagawa). The red and empty circles represent the selected and unselected sensors, respectively. Figure 12 shows how the spatial locations of the selected sensors for M4 change when the waiting time changes in forecasting the maximum on-shore tsunami amplitude in Onagawa. With the increase in waiting time, the locations of off-shore stations tend to be close to Onagawa. The summary MSE values for the Onagawa are also listed in Table 1. The consistency of the results for Iwanuma and Onagawa give greater confidence in the stability of the identified optimal configurations of the hazard-based tsunami early warning models.

Tsunami Loss for Iwanuma
The AIC score (Figure 13a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and decreases slowly to the minimum point (around 20 to 42 sensors). Figure 13b shows the MSE of M3 (base plus 33 sensors infor-

Tsunami Loss for Iwanuma
The AIC score (Figure 13a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and decreases slowly to the minimum point (around 20 to 42 sensors). Figure 13b shows the MSE of M3 (base plus 33 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 5 min for both models, similar to the optimal waiting time of 3 min using the hazard-based algorithm. Given the number of sensors added to the model and the waiting time, the selected sensors for M3 and M4 are shown in Figure 13c,d. The selected sensors for M3 cover the whole region of the S-net, whereas the selected sensors for M4 are located more evenly (one near the coastline, three in the middle, and two far from the coastline) compared to Section 3.1. evenly (one near the coastline, three in the middle, and two far from the coastline) compared to Section 3.1. Figure 14 shows how the spatial locations of the selected sensors for M4 change when the waiting time changes in forecasting the tsunami loss in Iwanuma. Table 1 lists the mean values of the MSE with regard to forecasting the tsunami loss in Iwanuma. The smaller values of MSE for M2, M3 and M4, compared to M1, clearly suggest that S-net sensors provide useful information for forecasting the tsunami risk in Iwanuma. The performance of M3 in forecasting the tsunami amplitude in Iwanuma is better than M4 (see Figure 13b and Table 1). The result indicates that a more comprehensive coverage of the seismic rupture areas facilitates the improved forecasting of the regional tsunami consequences, which tend to be more variable and exhibit more nonlinear increasing trends compared with the maximum on-shore tsunami amplitude.     Table 1 lists the mean values of the MSE with regard to forecasting the tsunami loss in Iwanuma. The smaller values of MSE for M2, M3 and M4, compared to M1, clearly suggest that S-net sensors provide useful information for forecasting the tsunami risk in Iwanuma. The performance of M3 in forecasting the tsunami amplitude in Iwanuma is better than M4 (see Figure 13b and Table 1). The result indicates that a more comprehensive coverage of the seismic rupture areas facilitates the improved forecasting of the regional tsunami consequences, which tend to be more variable and exhibit more nonlinear increasing trends compared with the maximum on-shore tsunami amplitude. Figure 13. Results of the risk-based early warning model for Iwanuma: (a) AIC score plot, (b) MSEwaiting time plot, (c) optimal sensor map for M3, and (d) optimal sensor map for M4. The pink star represents the target location (Iwanuma). The red and empty circles represent the selected and unselected sensors, respectively.

Tsunami Loss for Onagawa
The AIC score (Figure 15a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and decreases slowly to the minimum point (around 18 to 28 sensors). Figure 15b shows the MSE for M3 (base plus 24 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 10 min for both models. Given the number of sensors added to the model and the waiting time, the selected sensors for M3 and M4 are shown in Figure 15c,d. The selected sensors for M3 cover the S-net, whereas many of the selected sensors for M4 are located far from the coastline only with one close to Onagawa, and they are widely located.

Tsunami Loss for Onagawa
The AIC score (Figure 15a) decreases significantly from the base model to the model with the knee-point (around 5 to 10 sensors) and decreases slowly to the minimum point (around 18 to 28 sensors). Figure 15b shows the MSE for M3 (base plus 24 sensors information) and M4 (base plus 6 sensors information), and the suggested waiting time is 10 min for both models. Given the number of sensors added to the model and the waiting time, the selected sensors for M3 and M4 are shown in Figure 15c,d. The selected sensors for M3 cover the S-net, whereas many of the selected sensors for M4 are located far from the coastline only with one close to Onagawa, and they are widely located. Figure 16 shows how the spatial locations of the selected sensors in M4 change when the waiting time changes in forecasting the tsunami loss in Onagawa. As the waiting time increases, the locations of the off-shore stations tend to be close to Onagawa (same as the previous three analyses). The summary MSE (based on tsunami loss) values for Onagawa, summarized in Table 1, also indicate that the more spacious S-net off-shore sensors that are distributed, the higher the accuracy of the loss estimation can be. The consistency of longer waiting time in the risk-based results suggests that forecasting the tsunami loss can be an effective way to update the tsunami early warning.    Table 1, also indicate that the more spacious S-net off-shore sensors that are distributed, the higher the accuracy of the loss estimation can be. The consistency of longer waiting time in the risk-based results suggests that forecasting the tsunami loss can be an effective way to update the tsunami early warning.

Challenges in Practical Applications
The results shown in Section 3 are insightful in developing early warning systems based on an extensive OBS network. For the hazard-based models, the AIC-based model (M3) suggests that a set of 25-35 sensors that span across the rupture area is needed, while the knee-point-based model (M4) indicates that a smaller set of six sensors can be enough for establishing an effective early warning system. For the hazard-based model, a short waiting time of 3-5 min is sufficient, while for the risk-based model, a longer waiting time of 5-10 min will ensure stable forecasting. These results indicate that the warning system can be designed to release multiple updates at certain times. Early-phase warning messages could be focused on maximum tsunami wave amplitude, while at later times, the warning could be updated by focusing on possible tsunami consequences. Although issuing multiple warnings can cause some confusion, these warning procedures together with how to interpret these warning messages should also be clearly communicated with evacuees, since any warning messages cannot be free from false warning messages.
The short waiting time of 3 to 5 min for the hazard-based forecasting is only achieved by the extensive spatial coverage of the S-net. A conventional wave-monitoring network tends to have stations along the coastline, lacking the coverage in the off-shore rupture area. When the monitoring sensors are relatively far from the rupture area, a short waiting time is difficult to achieve. To demonstrate this more clearly, an additional case is investigated by restricting the available S-net stations to those close to the coastline. More specifically, since the number of selected sensors based on knee-point is six, the same number of sensors near the coastline are selected, which are shown in Figure 17a. These six sensors are the closest to Iwanuma and Onagawa. With regard to the hazard-based estimation, the model that relies on six nearshore sensors cannot achieve a forecasting performance of maximum on-shore tsunami amplitude in Iwanuma comparable to the AIC-based model until the waiting time of 20 min (Figure 17b). The same trend can be seen in the forecasting of the maximum on-shore tsunami amplitude in Onagawa (Figure 17c). The main reason for the longer waiting time is that the nearshore sensors simply cannot detect any meaningful early arrivals of potentially large tsunamis until the major tsunami waves reach these sensor locations. Note that such late detection could be fatal in issuing suitable tsunami early warnings for residents to achieve safe and timely evacuations.
In terms of risk-based forecasting, the MSE of the model with AIC-selected sensors decreases significantly in the first 10 min of waiting time and then reaches a flat trend.

Challenges in Practical Applications
The results shown in Section 3 are insightful in developing early warning systems based on an extensive OBS network. For the hazard-based models, the AIC-based model (M3) suggests that a set of 25-35 sensors that span across the rupture area is needed, while the knee-point-based model (M4) indicates that a smaller set of six sensors can be enough for establishing an effective early warning system. For the hazard-based model, a short waiting time of 3-5 min is sufficient, while for the risk-based model, a longer waiting time of 5-10 min will ensure stable forecasting. These results indicate that the warning system can be designed to release multiple updates at certain times. Early-phase warning messages could be focused on maximum tsunami wave amplitude, while at later times, the warning could be updated by focusing on possible tsunami consequences. Although issuing multiple warnings can cause some confusion, these warning procedures together with how to interpret these warning messages should also be clearly communicated with evacuees, since any warning messages cannot be free from false warning messages.
The short waiting time of 3 to 5 min for the hazard-based forecasting is only achieved by the extensive spatial coverage of the S-net. A conventional wave-monitoring network tends to have stations along the coastline, lacking the coverage in the off-shore rupture area. When the monitoring sensors are relatively far from the rupture area, a short waiting time is difficult to achieve. To demonstrate this more clearly, an additional case is investigated by restricting the available S-net stations to those close to the coastline. More specifically, since the number of selected sensors based on knee-point is six, the same number of sensors near the coastline are selected, which are shown in Figure 17a. These six sensors are the closest to Iwanuma and Onagawa. With regard to the hazard-based estimation, the model that relies on six nearshore sensors cannot achieve a forecasting performance of maximum on-shore tsunami amplitude in Iwanuma comparable to the AIC-based model until the waiting time of 20 min (Figure 17b). The same trend can be seen in the forecasting of the maximum on-shore tsunami amplitude in Onagawa (Figure 17c). The main reason for the longer waiting time is that the nearshore sensors simply cannot detect any meaningful early arrivals of potentially large tsunamis until the major tsunami waves reach these sensor locations. Note that such late detection could be fatal in issuing suitable tsunami early warnings for residents to achieve safe and timely evacuations.
will be several challenges related to missing data and large noises due to malfunctioning or maintenance of OBS sensors at the time of tsunamis. To ensure robust operation of the early warning system, multiple forecasting models can be developed for a given waiting time. As shown in Section 3, there are numerous configurations of sensors that achieve similar forecasting performance. The actual decisions for issuing warning messages can be based on an ensemble of tsunami hazard and risk forecasts based on different models, which are validated to perform equally well. In this way, the operator can minimize the chance of issuing false warning messages and avoid the failure of issuing the warning messages due to missing data and large noises. These practical aspects of the tsunami early warning systems should also be investigated in the future.

Limitations and Future Improvements
There are several limitations in the study. The first limitation is that the datasets are hypothetical. There are not enough large historical events of earthquakes and tsunamis that have occurred in a region with numerous off-shore sensors. Some bias in the results may be caused by utilizing such synthetic data. Vector quantities of the maximum onshore tsunami amplitude and tsunami loss can be forecasted, and the algorithm can be tested more comprehensively and globally in the future [41]. Regarding the statistical method, other advanced statistical methods, such as random forest and neural network [42], can be applied to forecast the tsunami impact. Moreover, multiple sources of tsunamis, such as combination of fault displacement and underwater landslide movement [43,44], can be considered in the future. In terms of risk-based forecasting, the MSE of the model with AIC-selected sensors decreases significantly in the first 10 min of waiting time and then reaches a flat trend. Whereas the MSE of the model with six near-coast sensors decreases significantly after waiting 20 min. The difference of the MSE values for the two models estimating the tsunami loss is more than estimating the tsunami amplitude (Figure 17d,e). The sensors measure wave amplitude that is more relevant to the hazard aspect of the tsunami (especially when the wave is near the coastline, it is closer to the on-shore tsunami amplitude), whereas the tsunami loss has more variability, which is difficult to forecast.
In implementing the developed tsunami early warning algorithms in reality, there will be several challenges related to missing data and large noises due to malfunctioning or maintenance of OBS sensors at the time of tsunamis. To ensure robust operation of the early warning system, multiple forecasting models can be developed for a given waiting time. As shown in Section 3, there are numerous configurations of sensors that achieve similar forecasting performance. The actual decisions for issuing warning messages can be based on an ensemble of tsunami hazard and risk forecasts based on different models, which are validated to perform equally well. In this way, the operator can minimize the chance of issuing false warning messages and avoid the failure of issuing the warning messages due to missing data and large noises. These practical aspects of the tsunami early warning systems should also be investigated in the future.

Limitations and Future Improvements
There are several limitations in the study. The first limitation is that the datasets are hypothetical. There are not enough large historical events of earthquakes and tsunamis that have occurred in a region with numerous off-shore sensors. Some bias in the results may be caused by utilizing such synthetic data. Vector quantities of the maximum on-shore tsunami amplitude and tsunami loss can be forecasted, and the algorithm can be tested more comprehensively and globally in the future [41]. Regarding the statistical method, other advanced statistical methods, such as random forest and neural network [42], can be applied to forecast the tsunami impact. Moreover, multiple sources of tsunamis, such as combination of fault displacement and underwater landslide movement [43,44], can be considered in the future.

Conclusions
This study developed an effective algorithm for tsunami early warning considering both hazard and risk aspects using sequential MLR. Utilizing off-shore sensors can improve the accuracy of forecasting the tsunami in terms of both tsunami amplitude and tsunami loss. The algorithm should use all off-shore sensors when they are available. Using the variable selection method implemented in this study, six off-shore sensors at specific locations improve the accuracy of forecasting tsunami most significantly. During the first 3 min after an earthquake occurs, the wave amplitude measured by the off-shore sensors contains enough information to estimate the maximum on-shore tsunami amplitude and the tsunami loss in Iwanuma A longer waiting time (5 to 10 min) can improve the forecasting of tsunami loss in Onagawa.