Speeding up and boosting tsunami warning in Chile

Chile host a great tsunamigenic potential along its coast, even with the large earthquakes occurred during the last decade, there is still a large amount of seismic energy to release. This permanent feature and the fact that the distance between the trench and the coast is just 100 km creates a difficult environment to do real time tsunami forecast. In Chile tsunami warnings are based on reports of the seismic events (hypocenter and magnitude) and a database of precomputed tsunami scenarios. However, because yet there is no answer to image the finite fault model within first minutes (before the first tsunami 5 wave arrival), the precomputed scenarios consider uniform slip distributions. Here, we propose a scheme of processes to fill the gaps in-between blind zones due to waiting of demanding computational stages. The linear shallow water equations are solved to obtain a rapid estimation of the run-up distribution in the near field. Our results show that this linear method captures most of the complexity of the run-up heights in terms of shape and amplitude when compared with a fully non-linear tsunami code. Also, the run-up distribution is obtained in quasi real-time as soon as the seismic finite fault model is produced. 10 Copyright statement. TEXT

This problem is separated in three parts.The determination of a seismic source model, the generation of an initial condition and the corresponding tsunami simulation.We define a computation domain around the earthquake source and the coastal areas in the near field.The bathymetric data used is the SRTM15 with 15 arcsec of resolution, based on the STM30 (Becker et al., 2009).
The idea here is to trade off accuracy for rapidness.In a near field earthquake tsunami context we care about the maximum 10 inundation place, the extension of the inundation until decreases to 0.5 to 1 m, and the average run-up.This model is not trying to be as accurate as possible and to determine a detailed inundation map but what intends to do is to give an idea of the main area that is going to be affected using the W-phase CMT, which is currently one of the fastest methodology to characterize large earthquake parameters in real time (Kanamori and Rivera, 2008).

The slip distribution model
Once an earthquake is firstly described by the W-phase solution, we use an elliptical slip distribution (Dmowska and Rice, 1986) over a region determined with the scaling laws obtained by Blaser et al. (2010).This serves as a first estimation while seismic waves are still traveling, and subsequent finite fault solutions are computed.Thus, we can model the near field tsunami for every finite fault model update.

The tsunami initial condition
Despite there are evidences of influence in the tsunami generation process with the source time components, for quickness purposes, we model static seafloor deformation induced by a non-uniform slip distribution including the horizontal components, as suggested in Tanioka and Satake (1996).This is obtained with that Okada's equations (Okada, 1985).

The tsunami modeling
The last part of this methodology is to obtain the tsunami heights along the coast.Usually, tsunami modeling involves complex codes to solve the fully coupled non-linear shallow water equations.Depending on the domain size and resolution, a full tsunami run can take several hours, which make real-time forecast 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 PDE with initial conditions: where η(x, y, t) is the water surface, g is the gravity acceleration, h(x, y) is the bathymetry and η 0 (x, y) is 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 an isobath of 100 m, before to reach the non-linear zone.Here, Neumann boundary condition is applied: ∂η ∂ n where n denotes the exterior unit normal vector.
The linear method (hereafter, LM) allows to obtain a quicker estimation than a full tsunami code since second-order terms are neglected but keeping the same main features.Also, this approach does not require to compute the velocity field, making the computation programs even faster.
Each simulation is compared with a fully linear water shallow water equation propagation.We select the JAGURS code (Baba et al., 2014) which runs in Fortran90 in a parallel array with OpenMP and OpenMP + MPI.This code is based in the classical finite difference method of Satake (2002).For each scenario, the tsunami is propagated during two hours of tsunami-

Tests and results
For evaluating the performance of this simple approach, we have modeled nearly all of the great tsunamis in the last two 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 exhibited in the supporting information.
For these examples we used the finite fault models obtained from USGS (Hayes, 2017), because they have proven to be operational and robust for real time operations in a global context monitoring.
All the computations were performed in a Dell Precision 7920, with two Intel Xeon Gold 6136 processor, 12 physical core each, for a total of 24 physical core, 2 threads each.For each time iteration, the domain is divided in 48 subdomains that are computed in different threads, for a parallel array (Fig. 1).
To compute the tsunami initial condition, the Okada's equations were implemented in C-language using threading, as well as the finite difference scheme for the LM.The C-code uses pthread library.With this tool, a data struct containing a pthread_t was defined, and then a routine that send every grid's subdomain to each different core's thread for computation.This method is such reliable as any other linear scheme method, because it solves the same equations.The only significant difference is the threads 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, with the used machine, in a regular grid of 4 million points with a FFM of 300 subfaults, the vertical and horizontal seafloor displacements can be calculated almost instantly (less than 5 seconds) and two hours of tsunami propagation for 2011 Tohoku, Japan earthquake can be solved in 60 seconds.

Examples
All the earthquakes presented here, have produced tsunamis.The range of magnitude varies from 7.7 to 9.1.For each event we apply the methodology described before.By using the W-phase centroid moment tensor, a scaling law, and a elliptic slip distribution the first source is defined.Then, the linear and non-linear 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 but instead considering a FFM solution.
Table 1 shows the correlation of the runup distributions between the code JAGURS and the presented here (the linear method).Table 2 summaries the CPU times for each stage of the process for each simulation.There is a high degree of agreement within a short time.
Detailed figures of the 24 simulations are provided in the supplementary material, where maximum amplitudes, runup distribution and field measurements are displayed.
4 Application to compliment tsunami alert.Case study: The 2015 Illapel Earthquake The September 16th a 8.3 Mw earthquake took place in the Coquimbo region in Chile (Melgar et al., 2016;Fuentes et al., 2017).
The features of the event were optimals for a tsunami generation.The national agencies deployed the protocols for evacuating  the whole chilean coast, even in distant insular territories (SNAM, bulletin #1, Sept. 16th, 23:02 UTC).The decision has to be made the very minutes after the origin time.In general, an accurate prediction of the tsunami runup heights requires a precise image of the seismic source, which nowadays is not available within 5 minutes and worst for real-time by adding the tsunami calculation times.Nevertheless, we can stay close of and quasi real time approach and trigger a first estimation based on the previous elliptical source.Because this can be done in a few seconds, it can be done at the moment instead searching in a precomputed data base of scenarios that are usually restricted.For monitoring purposes, the results can be updated every time a seismic source imaging is received, for both, the near field (at 15 arcsec) and regional field (at 60 arcsec).The whole information is summarized in a coded-color map, following the official coding from the Chilean institutions (Melgar et al., 2016).
Color coded maps are self-explicative, which make them easy to interpret (Figures 2 and 3).Each region can be rapidly assigned with a color which deploy an specific protocol for evacuation.All the simulations were performed with two hours of propagation, where the main energy content plays a key role on the inundation process.Figure 4 exhibits the normalized energy rate that generates the runup history along the coast, showing that majority of the global energy is concentrated within the first hour.We can also observe that the first estimation made with the elliptic fault predicts the same levels of inundation as the further finite fault model in the near field, and with minor differences in the regional field.This agrees the common sense where finite fault models act in the monitory stage and time is not as critic as in the very first minutes.
For completeness purposes, travel times isochrones are computed along the pacific basin (Figure 5).This is performed by computing a dense set of rays following Sandanbata et al. (2018), which allows to include dispersive effects.Also, we have included the effect of the earth elasticity as shown in An and Liu (2016).This kind of maps can be computed instantly with the very first estimation of the moment tensor and then update.

Discussion and Conclusions
Despite linear theory has been used in other tsunami warning systems, for instance, in the PTWC (http://unesdoc.unesco.org/images/0022/002203/220368e.pdf), now it is used in combination with more complex sources and faster algorithms with an 5 unique and simple product easy to interpret.
The non-complexity of the source, does not seem to be a problem when it comes to do a fast run-up estimation regarding emergency response.Establishing levels of hazard in real time, tries to generate a more accurate extension of the area affected by the tsunami, the maximum inundation that can reach, and how many people will be exposed to this hazard along the chilean coast.

10
This method intends to rapidly predict the run-up distribution.Using some simplifications in the tsunami equations it is possible to model rapidly the observed run-up.Obviously, this is not a high performance computing code, but even with many details outside the modeling such as complexity of the source, fine bathymetry and simplified equations, it can predict an important percentage of the run-up.The idea of extension even without the mathematical rigorosity we would like, is not inexact when we think of the idea of an emergency response system that needs to trigger actions that can save lives and reduce economic losses after the occurrence of a large earthquake.Using the methodology of Sandabata el al. it is possible to instantaneously calculate the tsunami arrival times from sources generated in the far field.This can also be done using tsunami modeling, however this takes a longer time of calculation then this numerical methodology is accurate enough to calculate tsunami travel times from far field events.
When we compare it with tsunami modeling codes, such as JAGURS, this method includes more than 80% of the predicted run-up in a 15 arcsec bathymetry at least 20 times faster.The main purpose of this study is to produce fast run-up estimations 5 with reasonable accuracy.
The idea here is to compliment the tsunami warning information while full and long calculations are being done (Figure 6).
This simple method provide simple and reliable information about the tsunami threat, which allows to entities to take decisions properly.

Figure 1 .
Figure 1.Sketch of the partition of the computation domain for parallel running.
Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2019-9Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License. Figure 2. Near Field Simulation of the 2015 Illapel Earthquake for elliptic source (left) and finite fault model (right).Red color represents the run-ups larger than 3 meters.Orange color is between 1 m and 3 m.Yellow color ranges between 0.3 m and 1 m.Green color codes run-ups smaller than 0.3 m.

1
Figure 3. Regional Field Simulation of the 2015 Illapel Earthquake for elliptic source (left) and finite fault model (right).Red color represents the run-ups larger than 3 m.Orange color is between 1 m and 3 m.Yellow color ranges between 0.3 m and 1 m.Green color codes run-ups smaller than 0.3 m. 8 Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2019-9Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 14 February 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 4 .Figure 5 .
Figure 4. Normalized run-up energy rate during the first two hours of tsunami simulation.Upper left panel display the run-up rate along latitude and time.Upper right panel shows the final maximum run-up and bottom left panel exhibits the normalized energy rate of the whole process as a time series.

Figure 6 .
Figure 6.Flow chart of the proposed methodology.

Table 1 .
Correlation of the runup distribution between the linear solution and the JAGURS code.
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 Mw 7.7 Nicaragua Earthquake, The 2006 Mw 7.6 Java earthquake.The extension of the earthquakes varies from L = 150 km to L = 500 km, the range of peak displacement at the source varies from 3 m to 40 m.Therefore, we have tested as many earthquakes and as many source features as possible for this study.

Table 2 .
Results of the CPU time of each event, in seconds.tIC is the time for computing the initial condition, tP r the time of processing, tT P the time of the tsunami propagation and tT is the total time.