Research article 28 Jun 2019
Research article  28 Jun 2019
Speeding up tsunami forecasting to boost tsunami warning in Chile
 ^{1}Department of Geophysics, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile
 ^{2}National Seismological Center, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile
 ^{3}Géoazur, Université de Nice Sophia Antipolis, Observatoire de la Côte d'Azur, Côte d'Azur, France
 ^{1}Department of Geophysics, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile
 ^{2}National Seismological Center, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile
 ^{3}Géoazur, Université de Nice Sophia Antipolis, Observatoire de la Côte d'Azur, Côte d'Azur, France
Correspondence: Mauricio Fuentes (mauricio@dgf.uchile.cl)
Hide author detailsCorrespondence: Mauricio Fuentes (mauricio@dgf.uchile.cl)
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 realtime 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 precomputed tsunami scenarios. Finite fault modeling, however, does not provide an estimation of the slip distribution before the first tsunami wave arrival, so all precomputed 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 timeconsuming finite fault model computations.We then solve the linear shallow water equations to obtain a rapid estimation of the runup distribution in the near field. Our results show that, at a certain water depth, our linear method captures most of the complexity of the runup heights in terms of shape and amplitude when compared with a fully nonlinear tsunami model. In addition, we can estimate the runup distribution in quasirealtime as soon as the results of seismic finite fault modeling become available.
 Article
(4226 KB) 
Supplement
(28513 KB)  BibTeX
 EndNote
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 realtime 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 Wphase 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 Wphase 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 wellknown, 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 precomputed 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 nearfield tsunamis triggered by earthquakes that complements the monitoring systems in operation and helps make better decisions during and after an emergency alert.
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 runup. 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 Wphase 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 Wphase 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 nearfield tsunami for every finite fault model update. The elliptic model is discretized with n_{y} subfaults alongdip and ${n}_{x}=\u2308\frac{L}{W}{n}_{y}\u2309$, 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.
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).
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 realtime 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 secondorder partial differential equation (PDE) with initial conditions:
where $\mathit{\eta}(x,y,t)$ denotes the water surface, g the acceleration of gravity, h(x,y) the bathymetry, and η_{0}(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 secondorder 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 runup distributions, despite the fact that later amplification of edge waves and resonance effects can occur. The approximated runup is obtained as the maximum from the vertical wall reflection boundary condition. The resulting runup values are on the same order of the actual runup for a sloping beach model (Synolakis, 1987).
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 runup 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 realtime 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.
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:

the 1992 M_{w} 7.7 Nicaragua Tsunami earthquake

the 2001 M_{w} 8.4 Southern Peru earthquake

the 2003 M_{w} 8.3 Hokkaido earthquake

the 2006 M_{w} 7.6 Java earthquake

the 2007 M_{w} 8.1 Solomon Islands earthquake

the 2007 M_{w} 8.0 Pisco earthquake

the 2009 M_{w} 8.1 Samoa Islands region earthquake

the 2010 M_{w} 8.8 Maule earthquake

the 2011 M_{w} 9.0 Tohoku earthquake

the 2012 M_{w} 7.8 British Columbia earthquake

the 2014 M_{w} 8.2 Iquique earthquake

the 2015 M_{w} 8.3 Illapel earthquake.
For each event we apply the methodology previously described, and use the Wphase 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 runup 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 runup 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, runup 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.
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 runup 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 quasirealtime 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 precomputed 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 colorcoded map following the official coding used by the Chilean institutions (Melgar et al., 2016). Colorcoded maps are selfexplanatory, 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 runup 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.
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 runups 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 runups 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.

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.

The noncomplexity of the adopted source does not seem to significantly affect the results of a fast tsunami runup estimation for emergency response purposes. By computing different levels of tsunami hazard in nearreal 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.

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.

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

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.
No data sets were used in this article.
The supplement related to this article is available online at: https://doi.org/10.5194/nhess1912972019supplement.
MF developed the idea and primary codes and tests. SA wrote all codes in C language with parallel computing and ran the simulations. SR compiled the catalogs of earthquakes used in this study and BD provided some FFMs to test the numerical tsunami model. The manuscript was prepared by MF and SR with supervision and contribution from all authors.
The authors declare that they have no conflict of interest.
This study was enterally supported by the Programa de Riesgo Sśmico.
This paper was edited by Maria Ana Baptista and reviewed by Victor Sardina and three anonymous referees.
An, C., Sepúlveda, I., and Liu, P. L. F.: Tsunami source and its validation of the 2014 Iquique, Chile, earthquake,, Geophys. Res. Lett., 41, 3988–3994, 2014. a
An, C. and Liu, P. L.: Analytical solutions for estimating tsunami propagation speeds, Coast. Eng., 117, 44–56, 2016. a
Baba, T., Takahashi, N., Kaneda, Y., Inazawa, Y., and Kikkojin, M.: Tsunami Inundation Modeling of the 2011 Tohoku Earthquake Using ThreeDimensional Building Data for Sendai, Miyagi Prefecture, Japan, in: Tsunami Events and Lessons Learned, Advances in Natural and Technological Hazards Research, edited by: Kontar Y., SantiagoFandiño V., and Takahashi, T., 35, Springer, Dordrecht, 2014. a
Becker, J. J., Sandwell, D. T., Smith, W. H. F., Braud, J., Binder, B., Depner, J., Fabre, D., Factor, J., Ingalls, S., Kim, S.H., Ladner, R., Marks, K., Nelson, S., Pharaoh, A., Sharman, G., Trimmer, R., Von Rosenburg, J., Wallace G., and Weatherall, P.: Global bathymetry and elevation data at 30 arc seconds resolution: SRTM30_PLUS, Mar. Geod., 32, 355–371, 2009. a
Blaser, L., Krüger, F., Ohrnberger, M., and Scherbaum, F.: Scaling relations of earthquake source parameter estimates with special focus on subduction environment, B. Seismol. Soc. Am., 100, 2914–2926, 2010. a
Doi, K.: Tsunami Warning System in Japan, in: Zschau, J. and Küppers, A., Early Warning Systems for Natural Disaster Reduction, Springer, Berlin, Heidelberg, 2003. a
Dmowska, R. and Rice, J. R.: Fracture theory and its seismological applications, Continuum Theories in Solid Earth Physics, 3, 187–255, 1986. a
Fuentes, M., Riquelme, S., Hayes, G., Medina, M., Melgar, D., Vargas, G., Gonzalez, J., and Villalobos, A.: A study of the 2015 Mw 8.3 Illapel earthquake and tsunami: Numerical and analytical approaches, in: The Chile2015 (Illapel) Earthquake and Tsunami, 255–266, Birkhäuser, Cham, 2017. a
Geist, E. L.: Complex earthquake rupture and local tsunamis, J. Geophys. Res.Sol. Ea., 107, ESE 21–ESE 215, https://doi.org/10.1029/2000JB000139, 2002. a
Gusman, A. R., Tanioka, Y., MacInnes, B. T., and Tsushima, H.: A methodology for nearfield tsunami inundation forecasting: Application to the 2011 Tohoku tsunami, J. Geophys. Res.Sol. Ea., 119, 8186–8206, 2014. a
Hayes, G.: The finite, kinematic rupture properties of greatsized earthquakes since 1990, Earth Planet. Sci. Lett., 468, 94–100, 2017. a
Ide, S., Beroza, G. C., Shelly, D. R., and Uchide, T.: A scaling law for slow earthquakes, Nature, 447, 76–79, 2007. a
Kanamori, H. and Rivera, L.: Source inversion of Wphase: speeding up seismic tsunami warning, Geophys. J. Int., 175, 222–238, 2008. a
Melgar, D., Allen, M., Riquelme, S., Geng, J., Liang, C., Bravo, F., Baez, J. C., Parra, H., Barrientos, S., Peng, F., Bock, Y., Bevis, M., Caccamise, Dana., Vigny, C., Moreno, M., and Smalley, R.: Local tsunami warnings: Perspectives from recent large events., Geophys. Res. Lett., 43, 1109–1117, 2016. a, b
Mulia, I. E., Gusman, A. R., and Satake, K.: Alternative to nonlinear model for simulating tsunami inundation in realtime, Geophys. J. Int., 214, 2002–2013, 2018. a
Okada, Y.: Surface deformation due to shear and tensile faults in a halfspace, B. Seismol. Soc. Am., 75, 1135–1154, 1985. a
Reid, R. O. and Bodine, B. R.: Numerical model for storm surges in Galveston Bay, Journal of the Waterways and harbors Division, 94, 33–58, 1968. a
Reymond, D., Okal, E. A., Hébert, H., and Bourdet, M.: Rapid forecast of tsunami wave heights from a database of precomputed simulations, and application during the 2011 Tohoku tsunami in French Polynesia, Geophys. Res. Lett., 39, L11603, https://doi.org/10.1029/2012GL051640, 2012. a
Riquelme, S., Fuentes, M., Hayes, G. P., and Campos, J.: A rapid estimation of near field tsunami runup, J. Geophys. Res.Sol. Ea., 120, 6487–6500, 2015. a
Riquelme, S., Bravo, F., Melgar, D., Benavente, R., Geng, J., Barrientos, S., and Campos, J.: W phase source inversion using high rate regional GPS data for large earthquakes, Geophys. Res. Lett., 43, 3178–3185, 2016. a
Riquelme, S., Medina M., Bravo, F., Barrientos, S., Campos, J., and Cisternas A.: Wphase Real Time Implementation and Network Expansion: The Experience in Chile 2012–2017s, Seismol. Res. Lett., 89, 2237–2248, 2018. a
Ruiz, J. A., Fuentes, M., Riquelme, S., Campos, J., and Cisternas, A.: Numerical simulation of tsunami runup in northern Chile based on nonuniform k2 slip distributions, Nat. Hazards, 79, 1177–1198, 2015. a
Sandanbata, O., Watada, S., Satake, K., Fukao, Y., Sugioka, H., Ito, A., and Shiobara, H.: Ray Tracing for Dispersive Tsunamis and Source Amplitude Estimation Based on Green's Law: Application to the 2015 Volcanic Tsunami Earthquake Near Torishima, South of Japan, Pure Appl. Geophys., 175, 1371–1385, 2018. a, b
Satake, K.: Tsunamis, in: Lee, W. H. K., Kanamori, H., Jennings, P. C., Kisslinger, C., International handbook of earthquake and engineering seismology, vol. 81A, Academic, Amsterdam, 437–451, 2002. a
Synolakis, C. E.: The runup of solitary waves, J. Fluid Mech., 185, 523–545, https://doi.org/10.1017/S002211208700329X, 1987. a
Tanioka, Y. and Satake, K.: Tsunami generation by horizontal displacement of ocean bottom, Geophys. Res. Lett., 23, 861–864, 1996. a
Wächter, J., Babeyko, A., Fleischer, J., Häner, R., Hammitzsch, M., Kloth, A., and Lendholt, M.: Development of tsunami early warning systems and future challenges, Nat. Hazards Earth Syst. Sci., 12, 1923–1935, https://doi.org/10.5194/nhess1219232012, 2012. a
Zhao, X., Duputel, Z., and Yao, Z.: Regional W phase source inversion for moderate to large earthquakes in China and neighboring areas, J. Geophys. Res.Sol. Ea., 122, 10052–10068, https://doi.org/10.1002/2017JB014950, 2017. a
 Abstract
 Introduction
 Methodology
 Implementation and benchmarking
 Discussion of computational results
 Application to compliment tsunami alert. Case study: the 2015 Illapel earthquake
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Methodology
 Implementation and benchmarking
 Discussion of computational results
 Application to compliment tsunami alert. Case study: the 2015 Illapel earthquake
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Review statement
 References
 Supplement