Tsunami Squares Simulation of Megathrust-Generated Waves: Application to the 2011 Tohoku Tsunami

Large subduction zone earthquakes often cause tsunamis, but observational data for hazard analysis is limited. Synthetic catalogs of seismically-generated tsunami sce-narios can be created by pairing earthquake and wave simulations. Tsunami Squares is one such wave simulator, explicitly tracking water mass and momentum, allowing simulation of dry land and the inundation process. We demonstrate a C++ port of Tsunami Squares paired with the Virtual Quake simulator by replicating the 2011 Great East Japan Earthquake and Tsunami. Comparisons of coastal run-up and wave heights with observations ﬁnds good agreement, with future improvements coming from tsunami source time dependence.


Introduction
Tsunamis present unique threats as natural disasters, because they are often accompanied by damage from the earthquakes or landslides that caused them. Seismic tsunamis are most often associated with large earthquakes in subduction zones. During such events, shifts in overhanging crustal slabs cause large seafloor displacements, resulting in corresponding displacements of water. As the ocean surface returns to equilibrium, these displacements propagate outward as tsunamis.

J o u r n a l P r e -p r o o f
There are many examples of massive loss of life and property as a result of seismic tsunamis in the Indo-Pacific region, with many of the highest-risk areas located in the developing world. Indonesia Lay et al., 2005), Cascadia (Satake et al., 1996;Nelson et al., 2006), Alaska (Shennan et al., 2009), and Chile  all feature subduction regions capable of producing tsunamis at a scale devastating to humans.
Numerical wave simulations are powerful tools for understanding tsunami hazard in at-risk areas. To reduce the impact of future disasters, physics-based simulations will play a vital role in statistical hazard analysis and tsunami early warning systems. Precomputed catalogs of synthetic tsunamigenic earthquakes and their resulting tsunamis, their inundation risks, and their observable warning signs can be used to identify hazards in time to issue warnings to threatened communities. Towards this goal, we introduce a simulator pipeline for tsunamigenic earthquakes and the coastal inundations they produce. We pair the Virtual Quake (Sachs et al., 2012;Schultz et al., 2018) longterm earthquake simulator with the Tsunami Squares wave simulator (Xiao et al., 2015;Wang et al., 2015b,a).
Typical far-field tsunami simulations employ the shallow water approximation of the Navier-Stokes fluid flow equations, the most common approach being some form of finite-difference method (Imamura, 1996). Tsunami Squares is an alternative numerical simulation algorithm. Originally developed in Fortran as a method for simulating landslide-induced tsunamis, Tsunami Squares is capable of seamlessly simulating dry land and inundation. With the goal of greater interoperability with other software, we introduce a C++ port of Tsunami Squares.
For proof of concept and validation purposes we recreate the 2011 Tohoku Earthquake and Tsunami. While Virtual Quake is not used to generate synthetic earthquakes in this study, the inferred subduction slab slip from the 2011 event is fed into Virtual Quake to recreate seafloor displacements, which are then used as initial conditions for a Tsunami Squares wave simulation. We compare waveforms and inundation to observations. We also compare the Virtual Quake/Tsunami Squares results to waveforms produced using an established simulator, the Regional Ocean Modeling System (https://www.myroms.org).
J o u r n a l P r e -p r o o f

Great East Japan Earthquake and Tsunami
On March 11, 2011, at 2:46 pm local time, a large slip event occurred in the Japan Trench roughly 72 km east of the Tohoku coast, where the Pacific plate subducts under the North American plate. The event lasted for six minutes, with maximum slips of up to 60 m over an area of 60, 000km 2 , and would eventually be categorized as a magnitude 9.1 earthquake. (Hayes et al., 2017). The quake had been preceded by several M > 6.0 foreshocks in the prior three days.
This earthquake, which represented the largest in recorded Japanese history, had a relatively shallow hypocentral depth of 25 km, and resulted in tsunamigenic sea floor displacement. The initial inundations resulting from the tsunami occurred roughly 10 minutes after the earthquake, with the most damaging inundation of the Sendai Plain continuing to occur an hour after the earthquake. The wave heights reached up to 35 meters, overwhelming the country's safeguards. The disaster resulted in an estimated death toll of 16,000 and caused approximately $210 billion in damages (Able M, 2012;Mori et al., 2012).
The Great East Japan Earthquake and Tsunami in many ways represents a worstcase natural hazard occurring in a country with best-case disaster risk reduction practices. The nonetheless great losses of life and property were in part accountable to a lack of preparedness for the least-likely scenario. Along with the wealth of observational data from -and analyses of -the event, this becomes an attractive case study for developing a pipeline for statistical hazard analysis and early warning.

Shallow Water Model
Typical tsunami flow calculations are performed by solving nonlinear long wave continuity equations. For a two-dimensional model with wave propagation through locations r r r = (x, y), in which water height is represented by the scalar H(r r r,t) and height-averaged velocity as v v v(r r r,t): J o u r n a l P r e -p r o o f and ∂ H(r r r,t)v v v(r r r,t) ∂t = −∇ · v v v(r r r,t)v v v(r r r,t)H(r r r,t) − gH(r r r,t)∇ζ (r r r,t) Here g is the gravitational acceleration, ζ (r r r,t) is the elevation of the water, and t is time. For a small, discreet time step dt, Equations 1 and 2 can be rewritten as − gH(r r r,t)∇ζ (r r r,t)dt (4)

Tsunami Squares Wave Model
The Tsunami Squares method solves equivalent equations by explicitly propagating mass and momentum across a grid of N cells with, for some cell i, center points r r r i , area A i , and topographic altitude T i . At time t, each cell contains a column of water of height H i (t), with vertically-averaged horizontal velocity v v v i (t)) and mean horizontal acceleration of a a a i (t). The mass associated with each cell is m i (t) = ρ w A i H i (t) (where ρ w is the density of water), and the momentum associated with each cell is p p p i (t) = Dry land is simply represented as cell with H i (t) = 0.
To begin a time step, the acceleration is calculated for each cell. We consider acceleration due to gravity and internal fluid drag. Gravitational acceleration is proportional to the gradient of the altitude of the top surface of the water, ζ i (t): where ν is a tuning parameter. Working against the direction of motion is frictional acceleration resulting from bulk material motion. This scales with the square of velocity: where µ d is a dynamic frictional coefficient for the material representing velocitydependent particle interaction. Frictional acceleration acts in opposition to velocity, J o u r n a l P r e -p r o o f but is never allowed to reverse the direction of flow. The total acceleration acting on a simulation cell is a a a i (t) = a a a i grav (t) + a a a i frict (t) Wave propagation simulation proceeds as illustrated in Figure 1. The water in each square is moved according to its velocity and acceleration, such that the water's new center point is located atr while its velocity and momentum are updated tõ After it moves, water from source-cell i will overlap with other simulation cells, which will receive whatever water overlaps their area. For receiving-cell j, the overlap proportion is denoted as δ A i j (t). Mass and momentum is distributed to those cells according to the proportion of overlap.
The momentum is similarly distributed: and the new velocity associated with each cell is calculated Because this method tracks volume and momentum from each time step to the next, these quantities are naturally conserved during simulation, enforcing J o u r n a l P r e -p r o o f Other fluids may be added to the simulation, for example mud in the case of a landslide. This requires frictional interactions with topography and other simulated materials, along with unique densities and internal frictional losses. Each time step, each material goes through an identical redistribution process to that outlined above, with momentum transferring between layers of materials. The C++ implementation of Tsunami Squares currently only supports simulation of water over topography.

Simulation
The basic program flow of Tsunami Squares is illustrated in Figure 2. Tsunami Squares begins by reading a file containing the combined topography/bathymetry of the region to be simulated. This defines the location, size, and altitude of each interacting square. Currently, the ETOPO1 combined bathymetry/topography map provided by the

National Oceanic and Atmospheric Association ([dataset] Amante and Eakins, 2009)
is used to define the pre-event sea floor at 1 arc minute resolution. Each square is filled with water up to sea level, with no water being added initially to squares above sea level. The initial conditions of the simulation are then set according to one of three options: 1. Vertically displacing the topography of each square according to a provided displacement file. Such displacements might be modeled from an earthquake.
These seafloor displacements propagate to the surface of the water, providing a nonzero water surface gradient for square accelerations.
2. Reading a previously-generated tsunami simulation output file, allowing simulations to be resumed. Simulation can be resumed on the same or different topographic grid.
3. Placing a Gaussian-distributed pile of water on top of the ocean surface. This is used for physics verification purposes, as described in further detail in Section 3.5.
J o u r n a l P r e -p r o o f After the simulation is initialized, the loop over time steps is entered. In each loop, the simulation state is written to the output file, squares are accelerated and moved, according to the physics of Section 3.2, and the heights and momenta of the squares is smoothed. When all time steps have been simulated, the final simulation state is written and the program terminates.
In methods like Tsunami Squares, smoothing is an important step to avoid nonphysical artifacts of the discrete simulation. In Tsunami Squares, this is accomplished by comparing the water heights and momenta in each square to those of its nearest neighbors. For each neighbor, the square with more water lends δ smooth * (H sel f − H neighbor ) to the neighbor with lower height, where Likewise, the square with more momentum contributes δ smooth * (p p p sel f − p p p neighbor ) to the lower-momentum square. The form of Equation 15, as well as its constants, was chosen for best performance across many different heights of water. The height dependence causes smoothing to be stronger in deeper water, and less effective near coasts, where variations in height and momentum between neightboring water columns J o u r n a l P r e -p r o o f are important for accurate inundation simulation. The min operation ensures that no more than 7.5% of the height or momentum difference is traded between squares each smoothing pass. The volume and momentum exchanges are calculated in parallel and then distributed to prevent a directional bias in smoothing. This, however, results in the occasional over-drawn square, which must be avoided to prevent simulation crashes (there's no such thing as negative amounts of water). Therefore, these "overdrafts" are filled back in, and the overall volume of water in the simulation is normalized such that total volume is conserved. Several smoothing passes each time step are usually necessary to produce artifact-free simulations, though the minimum necessary number of passes should be used.

C++ Port
The original implementation of the Tsunami Squares algorithm was written in Fortran by Ward (Xiao et al., 2015;Wang et al., 2015b,a). The algorithm has been rewritten in the widely-used language C++ for ease of use for a broader audience and interoperability with other software. Additional functionality has also been added: 1. The ability to load a previous simulation on a different-sized simulation grid, allowing waves distant from shore to be simulated on a coarser, computationally efficient grid, and then transitioned to a finer-resolution grid for accurate inundation simulation.
2. Use of the GeographicLib and Boost Geometry libraries, providing accurate distance and overlap calculations on a spherical or ellipsoidal Earth.
3. Parallelization of the most computationally intense sections of the code using OpenMP, allowing performance gains for users with multiple computational threads available.

Verification
In order to verify that the simulated waves are physically accurate, the output of the initial conditions of a Gaussian-distributed pile of water over a flat ocean bottom is compared to a semianalytic solution (Ward, 2011).

J o u r n a l P r e -p r o o f
For such Gaussian-pile initial conditions, the initial displacement from sea level of the ocean surface is written as u surf z (r r r, 0) = D c exp ( −r 2 /R 2 c ), and the displacement at time t is where J 0 is the zeroth Bessel function of the first kind, ω(k) = gk tanh(kH) is the dispersion relation for ocean depth H, and Figure 3 shows the performance of Tsunami Squares starting with the Gaussian-pile initial conditions, demonstrating good agreement between the two. Virtual Quake models the earth's crust as a homogeneous linearly elastic halfspace.
Faults are modeled as planar surfaces, representing discontinuities in that halfspace.
These are divided into smaller squares, or elements. These fault elements represent the fundamental interacting pieces of the simulation, analogous to each sliding block of the slider block model of Rundle and Jackson (1977).
Each fault element interacts with all others through elastic stress transfer, governed by quasi-static Green's functions (Okada, 1992). For each unit of slip, elements transfer shear and normal stresses to each other depending on their relative positions and orientations.
For indices running over the three Cartesian coordinates, the stress tensor change σ i j at a location x x x due to influences from all other locations x x x is given by Rundle et al. (2006) where s l (x x x ,t) is the amount of slip in direction l, and T kl i j (x x x − x x x is the Green's function tensor. Virtual Quake allows access to much of its functionality, including Okada Green's function stress and displacement equations, through the QuakeLib library. To reproduce the surface displacement caused by the Tohoku earthquake, we apply these tools to the slip distribution and fault model of Satake et al. (2013), as shown in Figure 6. The result is the vector displacement of the Earth's surface following the earthquake (Figure 7). The vertical component of this displacement constitutes the initial conditions for the ensuing tsunami simulation.
Currently, the initial seafloor uplift occurs instantaneously at the beginning of a Tsunami Squares simulation. Large subduction earthquakes actually rupture over the J o u r n a l P r e -p r o o f

Simulation Results
The simulation was initialized in the large region shown in Figure 7. The native one arc minute spatial resolution of the ETOPO1 bathymetry dataset was used, with a two second temporal resolution. The gravitational acceleration tuning parameter was set to ν = 0.5. At this resolution, four smoothing passes were used each time step. The wave was allowed to propagate for 800 simulated seconds.
After 800 seconds, the wave front was close to making landfall, so the simulation was transitioned to the smaller region shown in Figure 8 Table 1. (2017) using the Regional Ocean Modeling System (ROMS) (https://www.myroms. org). Because the smoothing algorithms used in simulation constitute a low-pass filter on the simulated waveform, we apply a fifth-order Butterworth low-pass filter on the observed height data.

Waveform Analysis and Comparison to ROMS
ROMS is an open-source tsunami model, which has been modified here with simplified ocean physics. Water is assumed to have constant density (as in Tsunami Squares). The horizontal resolution is about 5 km, with 30 levels of vertical resolution. Song et al. (2017) incorporate kinetic and potential energy imparted by the horizontal seafloor displacement, and used in-situ GNSS measurements of seafloor displacement as the initial condition of the ROMS simulation.
As a measure of agreement between simulation and observation, we follow Aida (1978). The amplitudes of the buoy observation x i for the first and second half-cycles (Leading tsunami wavefront peak and trough), as well as the wave arrival times, t 0i , are J o u r n a l P r e -p r o o f  Table 2: K, κ, t 0−c , and τ calculated for Tsunami Squares and ROMS found for each buoy i. Similar quantities are found for simulation: y i and t c . We define K i = x i /y i , representing a correction factor for the simulated amplitude. The geometric average over observation points, K, is found: A factor κ, representing fluctuation in K i over the different buoys, is found: The mean difference in observed and simulated arrival times across all four buoys, t 0−c , as well as its standard deviation, τ, is found. Table 2 shows the calculated values for both Tsunami Squares and ROMS. These are plotted in Figure 10, with κ and τ represented as error bars. Because data for buoy 205 is unavailable after the first wave peak, buoy 205 does not contribute to K or κ for the second half-cycle.
Similar performance in the amplitude correction factor K is seen for Tsunami Squares as for ROMS, with both simulators performing well. The notable difference in performance is in the arrival times, t 0−c . Figure 9 shows that, while Tsunami Squares has accurate arrival times for three of the buoys closest to the initial uplift, the wavefront is too late reaching the farthest buoy, 205. This indicates the gravitational acceleration parameter, ν, is tuned low. Stronger gravitational accelerations were experimented with (not shown here), resulting in far lower arrival time error τ, but which arrived at each station many minutes too early. This indicates that the most valuable future improvement to the C++ port of Tsunami Squares is the implementation of timedependence in the tsunami source, allowing for accurate arrival times at full-strength gravitational accelerations.
J o u r n a l P r e -p r o o f Tsunami Squares currently only uses the vertical uplift in its initial wave conditions. Song et al. (2008Song et al. ( , 2017 argues for the importance of horizontal seafloor motion during tsunamigenesis, both in additional bulk displacements of water due to bathymetric gradients, and also in kinetic energy transfer to the ocean. The incorporation of this horizontal energy transfer will also be a future improvement of the source model. Red cells indicate regions where inundation was observed but not recreated in simulation (a forecast "Miss"), green cells indicate agreement between observation and simulation (a forecast "Success"), and blue cells represent areas where inundation was not observed, but which were inundated in simulation (a forecast "False Alarm"). Most simulated inundation occurred in cells where inundation was observed. However, large amounts of coast had observed inundation where the simulated wave did not flow. Considering the locations where this was most prominent, namely in valleys along the coast in the northern simulation region, it is likely that significant inundation forecasting improvements can be made by using a higher-resolution topography dataset. While the NOAA ETOPO1 dataset is sufficient for waves distant to the coast, one arc minute resolution is insufficient for simulating the inundation and runup that occurs in narrow valleys.

Simulated Inundation
J o u r n a l P r e -p r o o f

Conclusion
This article proposes a marriage of the Virtual Quake earthquake simulator and the Tsunami Squares wave simulator for use in statistical tsunami hazard analysis and early warning. Tsunami Squares holds advantages over traditional finite difference methods of simulating the shallow water equations, including needing no special treatment for dry vs. wet cells. Tsunami Squares has been ported from Fortran to C++ in order to leverage GIS and geometry libraries, allowing for simulation on a spherical earth. The new port also allows for multiprocessing, and transitions of simulations between grids of different resolutions.
As a case study, we recreate the 2011 Tohoku earthquake and tsunami. We generate the seafloor uplift of the earthquake, which acts as initial conditions for Tsunami Squares, from the Okada greens functions and fault geometry tools packaged with the Virtual Quake earthquake simulator. We replicate the fault geometry and slip distributions outlined by Satake et al. (2013). The simulated ocean surface then reflects this seafloor displacement.
We validate our simulation against NOWPHAS buoy wave heights during the 2011 Tohoku Tsunami, as well as coastal inundation levels. In comparison to a simulation produced using the open source ROMS simulator, we find similar performance in wave amplitude, with improvement needed in wave arrival times. This will be rectified through the inclusion of time dependence in the tsunami source model.