Speeding up tsunami forecasting to boost tsunami warning in Chile

Despite the occurrence of several large earthquakes during the last decade, Chile continues to have a great tsunamigenic potential. This arises as a consequence of the large amount of strain accumulated along a subduction zone that runs parallel to its long coast, and a distance from the trench to the coast of no more than 100 km. These conditions make it difficult to implement real-time tsunami forecasting. Chile issues local tsunami warnings based on preliminary estimations of the hypocenter location and magnitude of the seismic sources, combined with a database of pre-computed tsunami scenarios. Finite fault modeling, however, does not provide an estimation of the slip distribution before the first tsunami wave arrival, so all pre-computed tsunami scenarios assume a uniform slip distribution. We implemented a processing scheme that minimizes this time gap by assuming an elliptical slip distribution, thereby not having to wait for the more time-consuming finite fault model computations.We then solve the linear shallow water equations to obtain a rapid estimation of the run-up distribution in the near field. Our results show that, at a certain water depth, our linear method captures most of the complexity of the run-up heights in terms of shape and amplitude when compared with a fully nonlinear tsunami model. In addition, we can estimate the run-up distribution in quasi-real-time as soon as the results of seismic finite fault modeling become available.

1 Introduction

For decades, countries exposed to coastal inundation have done a lot of work to develop their tsunami warning systems (Doi, 2003; Wächter et al., 2012). Most tsunamis are generated by large subduction earthquakes and landslides, which, owing to the characteristics of the tsunami source process, places a real-time tsunami forecast out of reach. Regular earthquakes follow a scaling law that links their energy release (seismic moment) to their duration (Ide et al., 2007). For instance, a regular 8.5 M_w earthquake can last for about 2 min, whereas we can consider tsunami generation nearly instantaneously after the source origin time. This implies that a robust tsunami warning system must integrate several systems that monitor different potential triggers such as earthquakes and volcanoes, among others. In the case of tsunamis generated by subduction earthquakes, it is essential to detect and characterize the seismic source. Today, the W-phase method is preferred for accounting for large earthquakes in Chile, which provides a first moment tensor solution within 5 min (Riquelme et al., 2016, 2018). As a matter of fact, the regional W-phase method is now running in real time in less than 5 min (Zhao et al., 2017). This method is based on waveform inversion theory; therefore it is necessary to have an important number of broadband seismometers in the regional field. The implementation of this method relies on robust seismic networks. This paper tries to illustrate the possibility of replication of these examples in other countries with tsunami threat produced by earthquakes in the near field. It is well-known, however, that tsunami heights are very sensitive to the spatial slip distribution of the seismic source (Geist, 2002; Ruiz et al., 2015). Even after having a finite fault model (FFM), the simulation of the tsunami propagation can take several hours depending on the desired level of resolution. This is the reason why the tsunami forecasts of many of the current warning systems are based on pre-computed scenarios (Reymond et al., 2012; Gusman et al., 2014; Mulia et al., 2018). Chile and Japan use this methodology for that purpose (https://www.jma.go.jp/jma/en/News/lists/tsunamisystem2006mar.pdf, last access: 5 December 2018). This methodology, however, ignores the complexity of the seismic source and solves only for uniform slip models. We propose a methodology applicable to near-field tsunamis triggered by earthquakes that complements the monitoring systems in operation and helps make better decisions during and after an emergency alert.

2 Methodology

We can separate this problem into three main parts: (1) the estimation of a seismic source model, (2) the generation of initial conditions, and (3) the corresponding tsunami simulation. We define a computation domain around the earthquake source and the coastal areas in the near field. We use the SRTM15 bathymetric data with a 15 arcsec resolution, based on the SRTM30 (Becker et al., 2009).

The core idea consists in trading off some accuracy to gain speed. Within the context of tectonic tsunamis generated in the near field we want to know the places with the maximum inundation, the extension of the inundation until it decreases to 0.5–1 m, and the average run-up. Our model does not aim at computing a detailed inundation map with the best possible accuracy, but rather to provide a fast estimate of the main area prone to inundation relying on the W-phase CMT, currently considered one of the fastest and more accurate methods to characterize the source of large earthquakes (Kanamori and Rivera, 2008).

2.1 Slip distribution model

Once a W-phase solution provides a characterization of an earthquake, we use an elliptical slip distribution (Dmowska and Rice, 1986) over a region determined by applying the scaling law after Blaser et al. (2010). This serves as a preliminary estimation while seismic waves are still traveling, and later finite fault solutions are computed. This in turn allows us to model the near-field tsunami for every finite fault model update. The elliptic model is discretized with n_y subfaults along-dip and $n_x=⌈\frac{L}{W}{n}_y⌉$, where L and W are the length and width of the fault plane obtained with the scaling law. After setting n_y=16, all the earthquake cases analyzed in our study have enough resolution on the source area.

Figure 1 Schematic showing the discretization of the calculation domain for parallel computation.

2.2 Tsunami initial conditions

Despite evidence of influence of the source time components in the tsunami generation process, for speed purposes we model a static seafloor deformation induced by a nonuniform slip distribution that includes the horizontal components, as suggested by Tanioka and Satake (1996). This is obtained by applying the Okada equations (Okada, 1985).

Figure 2 Near-field simulation of the 2015 Illapel earthquake with an elliptical source (a) and a finite fault model (b). The colors assigned to different areas indicate the expected run-ups in meters: (1) red for run-ups larger than 3 m, (2) orange for run-ups between 1 and 3 m, (3) yellow for run-ups between 0.3 and 1 m, and (4) green for run-ups smaller than 0.3 m.

Figure 3 Regional field simulation of the 2015 Illapel earthquake for an elliptical source (a), and a finite fault model (b). The colors assigned to different areas indicate the expected run-ups in meters: (1) red for run-ups larger than 3 m, (2) orange for run-ups between 1 and 3 m, (3) yellow for run-ups between 0.3 and 1 m, and (4) green for run-ups smaller than 0.3 m.

2.3 Tsunami modeling

The last part of this methodology is the estimation of the tsunami heights along the coast. Usually, tsunami modeling involves complex codes to solve the fully coupled nonlinear shallow water equations. Depending on the domain size and resolution, a full tsunami simulation run can take several hours, which makes real-time forecast nearly impossible. To overcome this limitation, we solve the linear shallow water equations with a forward finite difference scheme. The propagation inside the domain is governed by the second-order partial differential equation (PDE) with initial conditions:

$$\begin{array}{ll}{\displaystyle }{\mathit{\eta }}_{tt}-g\mathrm{\nabla }\left(h\mathrm{\nabla }\mathit{\eta }\right)& {\displaystyle }=\mathrm{0}\\ {\displaystyle }\mathit{\eta }(x,y,\mathrm{0})& {\displaystyle }={\mathit{\eta }}_{\mathrm{0}}(x,y)\\ \text{(1)}& {\displaystyle }{\mathit{\eta }}_t(x,y,\mathrm{0})& {\displaystyle }=\mathrm{0},\end{array}$$

where $\mathit{\eta }(x,y,t)$ denotes the water surface, g the acceleration of gravity, h(x,y) the bathymetry, and η₀(x,y) the initial condition. In the open boundaries, we set a radiation condition (Reid and Bodine, 1968), whereas in the solid boundaries (coasts) we impose full reflection in a vertical wall placed at the 100 m isobath, before reaching the nonlinear zone. Here, a Neumann boundary condition is applied: $\frac{\partial \mathit{\eta }}{\partial \widehat{\mathit{n}}}$, where $\widehat{\mathit{n}}$ denotes the exterior unit normal vector. The linear method (LM) allows us to obtain a faster estimation than a full tsunami code since second-order terms are disregarded while still accounting for the same main features. In addition, this approach does not require computation of the velocity field, an added benefit that makes the computation programs even faster. Each simulation is compared to its corresponding full nonlinear shallow water equation propagation. We use the JAGURS code (Baba et al., 2014) written in Fortran90 using parallel arrays via OpenMP and OpenMP + MPI. This code is based on the classic finite difference method of Satake (2002). For each scenario, we run the simulation for the equivalent of 2 h of tsunami travel time to obtain the main features of the run-up distributions, despite the fact that later amplification of edge waves and resonance effects can occur. The approximated run-up is obtained as the maximum from the vertical wall reflection boundary condition. The resulting run-up values are on the same order of the actual run-up for a sloping beach model (Synolakis, 1987).

Figure 4 Normalized run-up energy rate during the first 2 h of tsunami simulation. Panel (a) shows the run-up rate along latitude and time, panel (b) shows the final maximum run-up, and panel (c) shows the normalized energy rate for the whole process as a time series.

Figure 5 Tsunami travel times across the Pacific basin for the 2015 Illapel earthquake. Panel (a) shows the travel times after the shallow water equations, while the travel times in (b) include the effects of dispersion and the earth elasticity for a wave frequency of 2 mHz.

3 Implementation and benchmarking

To evaluate the performance of our approach, we modeled nearly all the largest tsunamis of the last 2 decades. Most of them were already tested with an analytical approach in (Riquelme et al., 2015). The details of the propagation and run-up distribution of the 12 events tested are presented in the supporting information. For these examples we used the finite fault models provided by the USGS (Hayes, 2017), as they have proven operationally robust for real-time operations in the context of global monitoring. All the computations were performed on a Dell Precision 7920, with two Intel Xeon Gold 6136 processors, each with 12 physical cores, for a total of 24 physical cores, and two threads each. For each time iteration, the domain is divided into 48 subdomains that are computed in different threads, for a parallel array (Fig. 1). To compute the tsunami initial condition, the Okada equations were implemented in the C programming language using threading, together with the finite difference scheme for the LM. The C code uses the pthread library to define a C data structure containing a pthread, and then calls a function that sends each grid subdomain to threads running in different cores for computation. This method is as reliable as any other linear scheme method, as it solves the same equations. The only significant difference is in the thread distribution for time optimization. When a thread finishes, it computes for a certain time step and it joins with the others in order to avoid miscomputations. For instance, on the system used to run our computations for a regular grid of 4 million points with an FFM of 300 subfaults, the vertical and horizontal seafloor displacements can be calculated almost instantly (less than 5 s), and 2 h of tsunami wave propagation for the 2011 Tohoku, Japan, earthquake can be solved in 60 s.

Table 1 Correlation of the run-up distribution obtained from our linear model solution and the JAGURS code. Correlation is computed with the standard Pearson coefficient. Details can be found in the Supplement.

4 Discussion of computational results

All the earthquakes presented here have produced tsunamis. The range of magnitude varies from 7.7 to 9.1. They occurred in different subduction zones around the world. The largest ones are Tohoku in Japan and Maule in Chile. All of them show a thrust mechanism except for the Samoa event in 2009, which is a normal event. There are a few tsunami earthquakes in this section such as the 1992 M_w 7.7 Nicaragua earthquake and the 2006 M_w 7.6 Java earthquake. The geometry of the earthquakes causative fault varies from L=150 km to L=500 km; the range of peak displacement at the source varies from 3 to 40 m. Therefore, we have tested as many earthquakes and as many source features as possible for this study:

1. the 1992 M_w 7.7 Nicaragua Tsunami earthquake

2. the 2001 M_w 8.4 Southern Peru earthquake

3. the 2003 M_w 8.3 Hokkaido earthquake

4. the 2006 M_w 7.6 Java earthquake

5. the 2007 M_w 8.1 Solomon Islands earthquake

6. the 2007 M_w 8.0 Pisco earthquake

7. the 2009 M_w 8.1 Samoa Islands region earthquake

8. the 2010 M_w 8.8 Maule earthquake

9. the 2011 M_w 9.0 Tohoku earthquake

10. the 2012 M_w 7.8 British Columbia earthquake

11. the 2014 M_w 8.2 Iquique earthquake

12. the 2015 M_w 8.3 Illapel earthquake.

For each event we apply the methodology previously described, and use the W-phase centroid moment tensor, a scaling law, and an elliptic slip distribution to define the first source. Then, the linear and nonlinear tsunami simulations are performed. The resulting run-up distributions are decomposed along latitude and longitude in order to compare both models. The same procedure is repeated, this time considering an FFM solution instead. Table 1 shows the correlation between the run-up distributions obtained with the JAGURS code (nonlinear method) and the method presented in this paper (linear method). Table 2 summarizes the CPU times in seconds for different stages of the process for each simulation. There is a high degree of agreement within a short time. Detailed figures showing the results for the 24 simulations are provided in the Supplement, where maximum amplitudes, run-up distribution, and field measurements are listed. For comparison purposes, for the event in 2014 in Chile, the DART station 32 401 registered 0.25 m of amplitude (An et al., 2014), where the linear method predicts 0.39 m for the elliptic source and 0.12 m for the FFM, whereas JAGURS gives 0.55 m for the elliptic source and 0.15 m for the FFM.

Table 2 Summary of the CPU time in seconds for the 12 events. t_IC indicates the time needed to compute the initial conditions, t_Pr the processing time, t_TP the time to compute the tsunami propagation, and t_T the total time.

Figure 6 Flow chart of the methodology proposed in this study.

5 Application to compliment tsunami alert. Case study: the 2015 Illapel earthquake

On 16 September 2017 an 8.3 M_w earthquake occurred in the Coquimbo region, Chile (Melgar et al., 2016; Fuentes et al., 2017). The characteristics of this event made it an ideal case study for tsunami generation. The national agencies implemented the established protocols for evacuating the whole Chilean coast, even the more distant insular territories (SNAM, bulletin 1, 16 September, 23:02 UTC −5). Such decisions have to be made within minutes of origin time. In general, an accurate prediction of the tsunami run-up heights requires a precise image of the seismic source, which at present is not available within 5 min for real time after adding the tsunami simulation times. Nevertheless, we can come close to a quasi-real-time approach by triggering a first estimation assuming an elliptical slip distribution. This only takes a few seconds, and can at present be performed instead of searching a pre-computed database of scenarios that are usually limited. For monitoring purposes, the results can be updated every time a seismic source image is received, for both the near field (at 15 arcsec) and regional field (at 60 arcsec). All this information is summarized in a color-coded map following the official coding used by the Chilean institutions (Melgar et al., 2016). Color-coded maps are self-explanatory, which makes them easy to interpret (Figs. 2 and 3). Each region can then be rapidly assigned a color linked to a specific evacuation protocol. All the simulations were performed for 2 h of tsunami propagation where the main energy content plays a key role on the inundation process. Figure 4 illustrates the normalized energy rate that generates the run-up history along the coast, showing that the majority of the global energy is concentrated within the first hour. We can also observe that the first estimation obtained for an elliptical fault predicts the same levels of inundation as the full finite fault model in the near field, while we can observe minor differences in the regional field. This makes sense since finite fault model results become available during the tsunami monitoring stage, when time is not as critical as in the very first minutes after origin time. Note it is possible to increase the number of warning levels, allowing us to find the optimal number of states for emitting and communicating the warning bulletin. In this study we choose the UNESCO standard. For completeness, we computed the travel time isochrones across the Pacific basin (Fig. 5). These computations use a dense set of rays following Sandanbata et al. (2018), which allows us to include dispersive effects. We have also included the effect of the earth elasticity as shown in An and Liu (2016). These kinds of maps can be computed instantly together with the very first estimation of the moment tensor and then updated.

6 Conclusions

In this study we propose a method that disregards the fine complexity of the seismic source while using fine bathymetric data and a set of simplified equations. Implementation of this method allows us to model more than 80 % of the tsunami run-ups with enough accuracy for tsunami warning purposes up to 20 times faster. Our method also aims at rapidly predicting the spatial distribution of the tsunami run-ups using some simplifications in the tsunami equations. Despite lacking the mathematical rigor that we would otherwise prefer, the method we propose is not inexact within the context of an emergency response system that needs to trigger actions that can potentially save lives and reduce economic losses after the occurrence of a large earthquake. We summarized our approach in the flowchart shown in Fig. 6. Taking into account the results of our study we can list the following as the most noteworthy results.

1. Although other tsunami warning centers use linear theory as part of their operations, for instance at the Pacific Tsunami Warning Center (PTWC) (http://unesdoc.unesco.org/images/0022/002203/220368e.pdf, last access: 20 December 2018), in this study we have combined it with the use of more complex sources and faster algorithms to generate a unique and simple product easy to interpret.

2. The non-complexity of the adopted source does not seem to significantly affect the results of a fast tsunami run-up estimation for emergency response purposes. By computing different levels of tsunami hazard in near-real time we can estimate more accurately the extent of the area potentially affected by the tsunami, the maximum level of inundation, and how many people will be exposed to this hazard along the Chilean coast.

3. Using the methodology of Sandanbata et al. (2018) it is possible to instantaneously calculate the tsunami arrival times from sources generated in the far field with enough accuracy. This can also be done via tsunami modeling, but at the expense of longer computation times.

4. When compared to other tsunami modeling codes such as JAGURS, results obtained from our method match more than 80 % of the predicted run-up for 15 arcsec bathymetry while obtaining the results at least 20 times faster.

5. The simple method proposed in this study provides a fast, reliable, and intuitive characterization of the tsunami threat, which in turn allows disaster mitigation agencies to take appropriate action.