TOWARDS THE NEXT GENERATION OF TSUNAMI IMPACT SIMULATIONS

The approach to tsunami modeling and simulation has changed in the past few years more than it had in the previous two decades. This brief review describes why this modeling shift is happening and attempts to provide some insight into the future of computational tsunami research.


Introduction
. Pieces of the Tohoku multi-billion dollar tsunami protection wall after a level-2 tsunami his the coasts of Japan in 2011. While each individual concrete section of the wall did not suffer significant damage, the erosion at the foundations of the wall happened so quickly that the concrete barrier simply fell (Picture from [1].) The failure of a multi-billion dollar wall designed to protect the Tohoku coasts of Japan ( Figure 1) from a level-2 tsunami (i.e. infrequent but highly destructive) in 2011 has triggered an important debate about alternative approaches to tsunami risk reduction. This debate is on-going and has attracted broad media attention worldwide (Reuters [2], The Guardian [3], The Economist [4], The New York Times [5], Wired [6]). The question whether a wall is the best solution to tsunami mitigation lies in the significant expense required for a wall that does not necessarily guarantee protection from big tsunamis. A massive concrete wall also suggests a misleading sense of full protection that may prevent the affected communities from prompt evacuation [7], as sadly experienced in Tohoku. Even when cost is not the main constraint-consider Japan whose G.D.P. accounts for the 4.22% of the world economy, versus Indonesia's 0.93%, or Chile's 0.24%-relying on traditional concrete based solutions alone may not be desirable or sustainable, partly because of their potential long-term negative impact on the population [2], coastal ecosystems [8][9][10], and shoreline stability [11,12]. For these reasons, decision makers and engineers are increasingly considering protection solutions that rely on green designs as sustainable and effective alternatives to seawalls. These designs, which are spreading along the Ring of Fire as shown in Figure 2, are usually man-made hillscapes erected on the shoreline to protect the communities behind them by partially dissipating and partially reflecting the tsunami energy [13]. Unlike walls, these protective solutions do not give a false expectation of full protection because they are not meant to block the flow from reaching shore; instead, they are designed to let the flow through while dissipating its energy and, hence, allow the affected population to gain precious evacuation time. It is expected that similar solutions are likely to be proposed for protection from storm surge to diminish the risk of the devastation that hurricanes as strong as the recent Sandy and Katrina may cause.
While the numerical capabilities to model tsunamis in the deep ocean are well established and mature, as it will be further pointed out throughout this review it is still not the case for the modeling of tsunami-shore interaction once the wave propagates inland. The limit lies in the level of detail necessary to correctly estimate the forces involved. In this respect, in only five years since the review by Behrens and Dias [14], important steps towards high fidelity tsunami modeling at all scales have been made, contributing to further advance the forecast of tsunami impacts as envisioned by Bernard et al.
[15] ten years earlier. This review describes the different approaches to tsunami modeling and simulation used to date. It will underline their properties and limits of validity to justify the ever increasing use of the 3D Navier-Stokes equations as the ultimate tool to study the tsunami boundary layer and run up.
The effects of isolated components such as bathymetry and vegetation have been studied for idealized one-dimensional (1D) and two-dimensional (2D) settings using numerical simulations for years. For example, in [13,27,[36][37][38], 1D and 2D shallow water models were used to demonstrate that the non-linearity of tsunami waves has a significant effect on the on-shore flow velocities and propagation. Many inundation modeling efforts aim to obtain an accurate prediction of the water elevation level. While these models typically do a good job at obtaining representative values of inundation levels, their correct prediction of water level seldom translates into the accurate prediction of the forces at play [39,40].
The evaluation of a tsunami from the point of view of wave breaking, turbulence, erosion, and sediment entrainment is computationally demanding. Why is this the case? A tsunami is a naturally multi-scale phenomenon whose characteristic length scales range from the planetary scale spanning oceans to the small scales of turbulence. In the open ocean it is, effectively, a fast moving two-dimensional long wave that can travel undisturbed for thousands of kilometers and can be accurately described by the 2D non-linear shallow water equations. On the other hand, as the flow approaches the shore and moves inland, its three-dimensionality and turbulent nature become important, with boundary layer dynamics that contribute to an important shear-driven erosion and sediment transport (See Figure 3). As the tsunami advances further inland, sand, dirt, vegetation, and other debris are scoured off of the sea bed, transported, and deposited along its path, potentially leading to important coastal morphological changes (e.g., [41][42][43][44][45][46]). By construction, however, the shallow water equations are not apt to model such complex flows, although they have been used in the context of erosion by [44,[47][48][49]. Towards an understanding of the possible limitations of the shallow water equations for a fine-grained tsunami modeling, Qin et al. [50] compared the 3D Reynolds-Averaged Navier-Stokes (RANS) solution of a turbulent propagating bore using OpenFOAM [51] against the 2D shallow water solution of the same problem using GEOCLAW [52]; Qin et al. demonstrated that a 3D model for turbulent flows is necessary to correctly predict tsunami inundation and the fluid forces involved. Among the first numerical studies that attempted to model tsunami run-up as the solution of the Navier-Stokes equations we find the 2D simulations by Hsiao and Lin [53] in 2011 and by Larsen and Fuhrman [54,55] who, in 2019, underlined the growing need for high-fidelity and multi-scale models to study tsunami risk and inland propagation. Furthermore, the RANS equations classically used for engineering turbulence modeling were shown to over-produce turbulence beneath surface waves [56], hence suggesting the need for high-fidelity turbulence models such as large eddy simulation. Large eddy simulation was used in the context of plunging breakers by Christensen [57] who underlined the numerical challenges of models in reproducing experimental results on the surf zone dynamics. Much attention has been given to run up; however, not enough work has addressed the analysis of the interaction of the tsunami with the built environment and vegetation, and of the erosion caused by large tsunamis [55] which is critical to improving the prediction of inundation. Modeling erosion with shallow water models is difficult due to the lack of velocity shear in these models; but shear is the driving mechanism of erosion and entrainment. Moreover, the modeling of the dynamics of the near-bed region and of erosion is still an open challenge in coastal engineering because it is governed by important particle-particle and particle-fluid interactions [58]; if analyzed by means of fine-grained numerical models, grain-resolving simulations are necessary [59,60].
Proper fine-grained modeling is important beyond erosion. In the context of tsunami impact mitigation analysis, the next frontier may lie in the detailed simulation of vegetation-tsunami and vegetation-morphology interactions. In light of the pioneering findings by the U.S. Army Corps of Engineers that vegetation may act as a bio-shield against flooding from storm surges [61,62], the designs of tsunami mitigation solutions around the world often incorporate vegetation. Significant experimental and numerical work has been done to analyse the effect of vegetation on wind waves, currents and, albeit less so, tsunamis (e.g., [63][64][65][66][67][68][69][70][71][72]), particularly in the context of steady flows [73]. After the 2004 tsunami in Sumatra, Bayas et al. [74] estimated that vegetation along the west coast of Aceh may have reduced casualties by 5%. The benefits of coastal vegetation are not limited to attenuation. It is well known that bathymetry, topography, and coastal geomorphology (rivers, canals, and barrier islands) have a profound effect on run-up for both tsunamis [13] and storm surge [75][76][77][78]. Furthermore, the presence of vegetation alters sediment-transport processes and landform evolution [79][80][81][82][83]. The modeling of the interaction between the tsunami and vegetation is numerically difficult because of the flexible nature of vegetation. Efforts have been made to include the effect of flexible vegetation on the flow. Examples of one-and two-way coupling of the flow with the the dynamics of vegetation can be found in the recent work by Mattis et al. [84,85] and Mukherjee et al. [86] who, for the first time, used fluid-structure interaction algorithms to study the effect of deformation of idealized trees on the flow (See Figure 4).

Mathematical representations
While the most complete model to describe a tsunami as a highly turbulent free surface flow is given by the Navier-Stokes equations of incompressible flows, simpler models of wave propagation have been successfully adopted for decades to study tsunami propagation in the open ocean and run-up [14,18]. In this section, we introduce the different models, underline their properties and limits of validity and justify the ever increasing use of the full 3D Navier-Stokes equations to study the boundary layer and run up.

Depth-Averaged Models
The first set of approximations that we will consider are the large class of equations characterized by averaging through the depth of the water column in the gravity centered coordinate system (See Fig. 5). This precludes modifications that would handle steeper terrain in favor of reducing complexity but should be noted as an alternative for simple topographical features. We will also forgo analysis of layered depth-averaged models as though they have unique properties and largely follow the same derivation. Figure 5. Gravity centered, depth averaged coordinate system for the shallow water equations. Here h is the depth of the water column, η the difference between a defined datum, commonly a given sea-level, and the sea-surface, and b the bathymetry surface as measured from the same datum.

Scaled Equations
We first start with the inviscid version of the 2D Navier-Stokes equations (i.e. the Euler equations) of incompressible flows, one vertical and horizontal dimension, where we will ignore one of the horizontal directions for simplicity of notation and without loss of generality. With this setup and assuming a free surface and non-flat bottom boundary we have the equations with the boundary conditions where the velocity is represented by u in the horizontal direction, w in the vertical, ρ the density of water, P the pressure, which includes non-hydrostatic components at this juncture, and g the gravitational acceleration. The first step in defining and justifying the depth-averaged equations lies in the non-dimensionalization assumptions made during their definition. Because of this we will be explicit about these assumptions as they will play an important role later on. The primary scalings are where λ is a length scale, often taken as the characteristic wave-length of the waves involved, a the characteristic depth that the wave is propagating through, T the characteristic time-scale or period of the wave, and P 0 the pressure normalization, usually taken to be atmospheric pressure. With these values defined we can also normalize the velocities and P such that where the non-dimensional quantities are those with·. Introducing then the traditional shallow water parameter = a λ we can write the velocity and temporal normalizations as Applying these scaling to (1) and simplifying we are lead to the system of equations where we have also assumed that P 0 ρU 2 = 1 without loss of generality. For the boundary conditions we introduce a new scaling, η δ =η and b β =b, which for the top boundary results in λ δw =η˜t +ũηx (9) λ βw =ũbx.
For the upper boundary condition, if δ λ = O( ), the boundary condition is all of the same order. This is probably not the case as δ a unless in very shallow water but also should be noted for later. The lower boundary condition is a similar situation, if β λ = O( ) the boundary condition is all of the same order. In the case of the bathymetry scale, this is probably the case as the bathymetry will vary on the scale of the depth scale a. The final non-dimensionalized equations are then (dropping thẽ ·): with the boundary conditions where P a is the appropriately scaled atmospheric pressure condition.

Depth Integration
The next step to deriving the shallow water and similar equations is to depth integrate (11). Here we will skip many of the details other than to point out two important assumptions that are often misrepresented.

of 17
Continuing on the first step is to integrate the incompressibility equation in (11) such that we find The next equation to be integrated is the horizontal momentum equation. For this we will first introduce the assumption that the total pressure P can be split up additively into an hydrostatic and non-hydrostatic component such that P(x, z, t) = P a + (η − z) + p where now p(x, z, t) is that non-hydrostatic component. Integrating the left-hand-side of this equation then we find where we have simplified the notation using our average values making sure to maintain the non-commutativity of the average and squaring operators. For the right-hand-side of the the horizontal momentum equation we similarly conclude that Finally we also integrate the vertical momentum equation through the depth. This is very similar to the horizontal equations leading to for the left-hand-side and − η b (P a + (η − z) + p) z + 1) dz ⇒ p| z=b for the right-hand side of the equation. Putting all this together leads us to the vertically integrated equations At this juncture it is useful to start to introduce notation that indicates average quantities through the depths. We will denote these by Here in lies one of the critical misrepresentations of depth averaged models, that they assume constant velocity, or even more general quantities, throughout their vertical profiles. For example, in general u 2 = u 2 . In fact we will make some assumptions about the commutativity of the averaging operator with others but this does not preclude the consideration of non-constant functions of the velocity for instance. This leads to the following system of equations, One last equation for the pressure will be useful as it will give us a means to calculate p depending on the approximations that we will make in the next section. We can come to this equation by reconsidering the scaled vertical momentum equation and integrating from η to a vertical level b + αh where α ∈ [0, 1] leading to

Approximations
At this point we can use the equations from the previous section to derive a number of approximations commonly used in tsunami numerical modeling. This will not be exhaustive but rather suggestive to where this basic framework can be taken to derive and analyze the validity of these approximations.
The first of these approximations will assume that the averaging operator will commute with multiplication, e.g. u 2 ≈ u 2 . The terms that this ignores are often termed as differential advection; but note that this approximation does not imply that the velocity profiles are constant but rather that the averages commute. This can be important when comparing boundary layer physics for instance. Moving forward, this allows use to rewrite (17) and simplify notation so that where we have also now dropped the · notation.

Shallow Water
Finally for the shallow water equations we simply need to assume that the shallow water parameter 1 implying that 0 = p| z=b . Since the non-hydrostatic pressure is defined such that p(x, η, t) = 0 we can show that p = 0 therefore implying that (17) become which are of course the traditional, non-dimensionalized shallow water equations. Note that (2.1.2) can be used here to show that p = 0 is in fact true here.

Not So Shallow Equations
The next perhaps logical step is to maintain the assumption that differential advection is still small and therefore averages commute but that vertical momentum is not as ignorable as before, i.e.
≥ O(1). Instead we might make some assumptions about the depth profile of w and use that to derive a form of a closure expression for the non-hydrostatic pressure. Here we will give one example of an analysis where we assume that w has a linear profile through the depth along with the equations that arise due to this assumption.
First if w is linear in z we can use the incompressibility condition to find The crux of these calculations are then to compute the expressions for w and p. Leaving out the tedious details of these calculations but reporting the results the final equations with these assumptions leads to where

Navier-Stokes equations
All of the models above do a very good job at solving wave propagation and, when enhanced with an appropriate wetting-and-drying scheme, at simulating inundation. Nonetheless, they lack the three-dimensionality that allows for turbulence to be accounted for. During inundation, the interaction of the flow with the coastal features, erodible terrain, sediments, vegetation, and structures is such that the flow is fully turbulent, hence three-dimensional and characterized by shear. The correct estimation of the hydrodynamic forces responsible for the damages during run up requires a model that can capture these flow characteristics. The Navier-Stokes equations are the most complete model available to us if these flow features are to be considered. The 3D continuity and Navier-Stokes equations of incompressible viscous flows can be easily obtained by making the Euler equations 1, time dependent, by extending them to three spatial dimensions, and by including molecular viscosity. We write them as follows: While Equations (21) do model turbulent incompressible flows, its solution on a computational grid requires some attention with respect to the treatment of the viscous stresses. If we could afford to have an infinitely fine computational grid, we could solve them directly; an approach, this one,  [87,88]. But it took approximately one more decade to become a popular choice to model tsunami generated turbulence. Examples of RANS tsunami models are described in [54,55,[89][90][91]. In terms of solution cost, LES lies between DNS and RANS. Due to its computational cost, it is not popular among tsunami modelers yet although it was recently considered as an optimal choice to study breaking waves and run up in highly non-steady dynamics (e.g. [86]). As the era of exa-scale computing on hybrid CPU-GPU architectures is fast approaching, we expect LES to become the model of choice among tsunami inundation modelers. Some of the most common tsunami simulators are listed in Table 1.
The free surface in a model based on the solution of (21) is typically accounted for in either one of two ways: by means of a volume of fluid approach (in, e.g. [91,92] ) or by a level-set tracking of the water-air interface in a two-fluid model (in, e.g., [86,93,94]). Both of these approaches make it natural to handle the wet-dry front as the flow reaches land; this is not the case for shallow-water equations whose numerical solution is challenged in wet-dry regions.

Numerical Solution
The equations described across this paper have been solved for decades using a wide range of numerical approximations. From the classical finite difference method (FD) used in, e.g. MOST [99], to 1 The number of grid points in a computational grid necessary for DNS scales as the (9/4) th power of the Reynolds number.  [25,52,95,107,118]), to the finite element method (FEM; e.g. [86,119]) and high-order continuous and discontinuous Galerkin methods (CG/DG; e.g. [13,26,33,94,106]), to the less common smoothed particle hydrodynamics method (SPH; e.g. used in GPUSPH [116]) that has been gaining popularity to model free surface flows in engineering and hydraulics from its origins in astrophysics.
Each numerical method comes with its advantages and disadvantages. FD are less costly than CG/DG, but they lack the geometrical flexibility of Galerkin methods and FV. CG/DG and FV are optimal for grid refinement [26,52,120] and are geometrically flexible to handle complex coastal lines; furthermore, CG/DG are inherently optimal for massive parallelism [121,122]. SPH is good at tracking the free surface of the three dimensional flow, but it is expensive and it is still unclear how to treat boundary layers. Regardless of the numerical method of approximation, the handling of the wet-dry interface to model inundation is still a difficult problem to solve, especially so for high-order methods as demonstrated by the extensive work on wetting and drying presented in, e.g., [31][32][33][34][35]123].
In summary, all these methods have their use case and there is no one that can rule them all. Similarly the adage "using the right tool for the job" holds here as much as it does when doing carpentry. It is important therefore to consider what metrics may lead to the decisions to chose one model over another. Among the most prominent are: • Performance on the computing architecture being considered.
• Complexity of the physics needed. Here the difference between using a wall boundary condition at the shore and doing true wetting and drying may be significant as does the representation of true turbulent flow. • Overall size of the problem considered. Are you interested in a single simulation or a large ensemble in order to take into account uncertainty.
These are of course only some of the issues but highlight the difficulties faced and the need for hybrid solver techniques to bridge the gaps between the advantages and disadvantages between solvers.

Towards a multi-scale framework from source to impact
It seems plausible, given the on-going research on all fronts from tsunami generation, propagation, to fine-scale run up, that the community of tsunami modelers will soon join forces to build a unified, all-scale framework to model a full tsunami event from source to impact (see, e.g. [21,124]). The all-scale idea is naïvely represented in Figure 6. The 3D Navier-Stokes solver is forced at the boundaries by the shallow water model of the tsunami moving off-shore, which is forced by an earthquake simulator or, in the case of landslide triggered tsunamis, landslide model. Arguably the most difficult and still unclear part of a modeling infrastructure of this type lies in the coupling across models. Coupling across models is an open field of research of its own (see, e.g. [125][126][127][128].) Parallel performance, data exchange, and time scale interactions, are among the difficulties to be overcome in designing coupling algorithms across software packages. We envision a major effort to optimize software coupling and make it efficient for costly tsunami simulations towards real-time modeling of a full tsunami event.

Conclusions
The current state of the field seems to indicate that the numerical tsunami modeling community is moving towards the high-fidelity modeling of tsunamis across scales. We expect an increasing effort to couple different models and allow the interaction of an earthquake simulator with the large scale dynamics in the open sea and fine scale dynamics in the coastal region, all within the same software infrastructure. While improved earthquake-tsunami modeling at the source, adaptive grid refinement, and ever faster and inexpensive hybrid high-performance computing (i.e. graphics processing units, GPUs) have contributed to the advancement of tsunami modeling at large scales, we are now witnessing the very beginning of a new tsunami modeling era; the introduction of high-fidelity simulations of inundation and enhanced coupling of software packages across all scales is finally Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 October 2020 doi:10.20944/preprints202010.0394.v1 Figure 6. Representation of a comprehensive off-to-on-shore flow framework. The 3D Navier-Stokes solver is forced at the boundaries by the shallow water model of the tsunami off-shore, which is forced by an earthquake simulator.
leading the way towards real-time tsunami forecast as envisioned by Synolakis et al. [18] fifteen years ago.
While it is difficult to estimate how long it will be before a full scale tsunami simulation can be run in real time (or faster) on a laptop, the pace at which the necessary tools are being developed is fast. At risk of being too optimistic, current research leads us to think that such a tool may be available within the next fifteen years.