Modeling and Simulation of Tsunami Impact: A Short Review of Recent Advances and Future Challenges

: Tsunami modeling and simulation has changed in the past few years more than it has in decades, especially with respect to coastal inundation. Among other things, this change is supported by the approaching era of exa-scale computing, whether via GPU or more likely forms of hybrid computing whose presence is growing across the geosciences. For reasons identiﬁed in this review, exa-scale computing efforts will impact the on-shore, highly turbulent régime to a higher degree than the 2D shallow water equations used to model tsunami propagation in the open ocean. This short review describes the different approaches to tsunami modeling from generation to impact and underlines the limits of each model based on the ﬂow régime. Moreover, from the perspective of a future comprehensive multi-scale modeling infrastructure to simulate a full tsunami, we underline the current challenges associated with this approach and review the few efforts that are currently underway to achieve this goal. A table of existing tsunami software packages is provided along with an open Github repository to allow developers and model users to update the table with additional models as they are published and help with model discoverability.


Introduction
The failure of a multi-billion dollar wall designed to protect the Tohoku coasts of Japan ( Figure 1) from a level-2 tsunami (Level-2 tsunami: infrequent but highly destructive [1]) in 2011 triggered an important debate about alternative approaches to tsunami risk reduction. This debate is ongoing 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 guarantee protection from big tsunamis. Even when cost is not the main constraint-consider Japan whose GDP accounts for 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 [7][8][9], and shoreline stability [10,11]. 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, three examples of which are shown in Figure 2 along the Ring of Fire, are usually human-made hillscapes erected on the shoreline to protect the communities behind them by partially dissipating and partially reflecting the tsunami energy [12]. The question remains whether these alternatives truly protect the people and properties behind it. Robust computational simulation capabilities are required that can model such diverse types of infrastructure including the forces upon them. 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 taken from in [13]).

Figure 2.
Map of the Ring of Fire, a long coastal stretch that is most likely impacted by large tsunamis. Some tsunami mitigation parks are being constructed in South Java, Indonesia (Image: Indonesia Ministry of Marine Affairs and Fisheries), Miyagi Prefecture, Japan (Image: the Morino Project), and Constitución, Chile (Image: Architect Magazine). Adapted from the work in [12].
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, this 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, especially on complex structures. 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, especially with respect to the fine-scale simulation of tsunami run-up and tsunami-shore interaction (see, e.g., in [15][16][17][18]), contributing to further advancing the forecast of tsunami impact as envisioned by Bernard et al. [19] ten years earlier.
In light of the development of a comprehensive, multi-scale tsunami modeling infrastructure as proposed in the latest 10-year science plan of the Natural Hazards Engineering Research Infrastructure [20], this article reviews models of tsunami propagation from generation to impact with respect to their régimes of validity and model-to-model interaction. As an aide to this, a basic derivation of the assumptions that lead to the basic two-dimensional shallow water equations are provided along with some of the most direct extensions and their limitations. Connecting these assumptions and their limitations to the much more complex three-dimensional and non-hydrostatic models follows and demonstrates the complexity the current state-of-the-art faces.
The article is organized as follows. Section 2 reviews the state of the field in tsunami modeling and simulation. The different mathematical models of classical use are described in Section 3, followed by a review of the different numerical methods to solve them given in Section 4. The idea of a future comprehensive and multi-scale simulation framework is discussed in Section 5. Conclusions are given in Section 6. A table of some available tsunami software packages is given in the Appendix A.

State of the Field in Tsunami Forward Modeling
The multi-scale nature of tsunami dynamics makes their study in a laboratory setting challenging [21,22]. For this reason, and supported by ever cheaper computing power and data storage, numerical modeling and simulation has become the most widely utilized tool for large-scale tsunami modeling and analysis [14]. The numerical simulation of tsunamis started in Japan sixty years ago with the work by Aida [23,24] and has shown to be effective to model their generation (see, e.g., in [25][26][27][28][29]), propagation (see, e.g., in [30][31][32]), and inundation [12,16,[33][34][35], although the problem of inundation is still partially open when it comes to its numerical treatment (see, e.g., in [36][37][38] and references therein).
When it comes to tsunami-shore interaction, 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 [12,33,[39][40][41] 1D and 2D shallow water models were used to demonstrate that the nonlinearity 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 this, their correct prediction of water level seldom translates into the accurate prediction of the forces at play [42,43], possibly indicating the limitations of 1D and 2D models in specific flow régimes.
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. On the one hand, in the open ocean it is effectively a fast moving two-dimensional long wave that travels undisturbed for thousands of kilometers and can be accurately described by the 2D nonlinear 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, structures, and other debris are scoured off of the sea bed, transported, and deposited along its path, potentially leading to important coastal morphological changes (see, e.g., in [44][45][46][47][48][49][50]). By construction, however, the shallow water equations are not apt to model such complex flows, although they have been used or enhanced in the context of erosion and sediment transport in [47,[51][52][53][54], for example.
Towards an understanding of the possible limitations of the shallow water equations for a fine-grained tsunami modeling, Qin et al. [15] compared the 3D Reynolds-Averaged Navier-Stokes (RANS) solution of a turbulent propagating bore using OpenFOAM [55] against the 2D shallow water solution of the same problem using GEOCLAW [56]; 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 [57] from 2011 and by Larsen and Fuhrman [17,18] who recently underlined the growing need for high-fidelity and multi-scale models to study tsunami risk and inland propagation. Furthermore, the RANS equations classically used in engineering were shown to overproduce turbulence beneath surface waves [58], therefore suggesting the need for high-fidelity turbulence models such as large eddy simulation (LES). LES was used in the context of plunging breakers by Christensen [59] 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 [18], all critical at improving the prediction of inundation. Modeling erosion with shallow water models can be difficult as they generally require some sort of reconstruction of vertical velocity profile to reconstruct the velocity shear in these models; however, 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 [60], which require grainresolving simulations [61,62]. In a recent inter-comparison of models to study shear and separation driven coastal currents, Lynett et al. [63] concluded that "[. . . ]In general, we find that models of increasing physical complexity provide better accuracy, and that low-order three-dimensional models are superior to high-order two-dimensional models[. . . ]".
Proper fine-grained modeling is important for studying 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 [64,65], the designs of tsunami mitigation solutions around the world often incorporate vegetation. Significant experimental and numerical work has been done to analyze the effect of vegetation on wind waves, currents and, albeit less so, tsunamis (see, e.g., in [66][67][68][69][70][71][72][73][74][75][76]), particularly in the context of steady flows [77]. After the 2004 tsunami in Sumatra, Bayas et al. [78] 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 tsunamis [12] and storm surge [79][80][81][82]. Furthermore, the presence of vegetation alters sediment-transport processes and landform evolution [83][84][85][86][87]. 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 dynamics of vegetation can be found in the recent work by Mattis et al. [88,89] and Mukherjee et al. [90], who, for the first time in the context of inundation, used fluid-structure interaction algorithms to study the effect of deformation of idealized trees on the flow (See Figure 4).

Mathematical Representations and Assumptions
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,91]. In this section, we introduce these models and underline their properties and limits of validity with the aim to identify the assumptions that go into them and where those assumptions may break down, especially in the near-shore. Based on this analysis, we justify the ever increasing use of the 3D Navier-Stokes equations to study the tsunami boundary layer in the near-shore and during run-up where 1D and 2D models were shown to lead to inaccurate solutions [15,63], leading to the idea suggested in Figure 5 of having both the 2D shallow water and 3D Navier-Stokes models working together to provide a full picture of a tsunami's impact. 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. In fact, the 2D shallow water solver can extend all the way into the shore depending on the needs of the solver and compatibility layer.

The 3D Navier-Stokes Equations
We first start with arguably the best model of water flow, the three-dimensional Navier-Stokes equations. 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, and thus 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 [17,63]. The Navier-Stokes equations are the most complete model that can capture these features. Omitting their derivation from first principles, the 3D Navier-Stokes equations of incompressible flows are written as

Depth-Averaged Models
Although (1) are the model equations for tsunamis that we would like to use, they are computationally too costly to use everywhere and more than likely not needed everywhere in the domain of interest for a tsunami. We therefore need to carefully identify where and when we can approximate (1). The primary type of approximation associated with tsunamis are characterized by averaging through the depth of the water column in the gravity centered coordinate system (see Figure 6). 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. 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 two-dimensional Navier-Stokes equations (i.e., the incompressible 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 and w in the vertical, ρ is the density of water, P its pressure-which includes non-hydrostatic components at this juncture, and g is the gravitational acceleration. The underscores t, x, and y indicate the partial derivatives with respect to them. The first step in defining and justifying the depth-averaged equations lies in the nondimensionalization 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; and T the characteristic time-scale or period of the wave. With these values defined, we can also normalize the velocities and P such that where P 0 is the pressure normalization, usually taken to be atmospheric pressure, and the non-dimensional quantities are those with·. Introducing then the traditional shallow water parameter = a/λ we can write the velocity and temporal normalizations as U = √ ga, W = U, and T = λ/U. Applying these scalings to (2) 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 where δ is surface amplitude scaling. 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 nondimensionalized equations are then (dropping the·) 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 (7)-(9). Here, we will skip many of the details other than to point out two important assumptions that are often misrepresented. Continuing on the first step is to integrate the incompressibility Equation (7) such that we find The next equation to be integrated is the horizontal momentum Equation (8). 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 the term (η − z) represents the depth from the surface η and therefore the non-dimensionalized hydrostatic component of the pressure and p(x, z, t) non-hydrostatic pressure component. Then, integrating the left-hand-side of (8), 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 (8) we similarly conclude that Finally, we also integrate the vertical momentum Equation (9) 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 the notation that indicates average quantities through the depths. We will denote these by

Remark 1.
Herein 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; however, 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 us to rewrite (11) 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 . As the non-hydrostatic pressure is defined such that p(x, η, t) = 0 we can show that p = 0 therefore implying that (11) become which are of course the traditional, non-dimensionalized shallow water equations. Note that (12) 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 w(x, z, t) = w(x, t) − u x z. The crux of these calculations is then the computation of w and p. Leaving out the tedious details of these calculations but reporting the results the final equations with these assumptions leads to where p = 2 h ∂ ∂t Note that this complexity is a result of the attempt to recreate the 3D structure that we have already integrated out. This may be counterintuitive as to why we might want to do this but in fact schemes of this nature have been quite successful, namely, the Boussinesq, Serre-Green-Naghdi, and other types of semi-2D types of equations.

Mathematical Conclusions
As we have seen, (1) can be greatly reduced in complexity by integration and appropriate assumptions. While this means that the 2D equations can easily recreate the extent of flooding, the recreation of the 3D velocity fields and therefore the appropriate forces on structures is however difficult even for schemes that have been successfully implemented such as the Boussinesq and Serre-Green-Naghdi models, leaving the need for 3D Navier-Stokes models as the gold standard for impact questions. Furthermore, as we have shown, the assumptions we have made are often not as stringent as is often thought in the literature. The importance of this is that the shallow water 2D field does not preclude there being a complex horizontal velocity field, rather that its average have the special property that it commutes as mentioned above. This critical property prescribes exactly the compatibility condition that we need to have between a 2D shallow water model and a 3D Navier-Stokes model whose field does not require this commutativity property.
Each numerical method comes with its advantages and disadvantages. FD are the least costly, but they lack the geometrical flexibility of Galerkin methods and FV. FE/DG and FV are optimal for grid refinement [32,56,100] and are geometrically flexible to handle complex coastal lines; furthermore, FE/SE/DG are inherently optimal for massive parallelism [101]. 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 with it. Regardless of the numerical method of approximation, the numerical handling of the wet-dry interface to model inundation has been improved throughout the years as demonstrated in, e.g., [36][37][38][102][103][104][105][106]. It should also be noted that new methods in high-performance computing continue to push the envelope forward and will consequently lead to advances in tsunami science. This advancement will be either hampered by the same issues encountered in general geophysical turbulence computational modeling, or enjoy the benefits of embarrassingly parallel ensemble modeling as is the case with probabilistic tsunami hazard assessment (PTHA).
While Equation (1) models turbulent incompressible flows, their solution on a computational grid using any of the methods above (except for SPH) requires some attention with respect to the treatment of the viscous stresses. If we could afford to have an infinitely fine computational grid, they could be solved directly; this approach is known as direct numerical simulation (DNS). However, even with the approaching era of exa-scale computing, DNS is much too costly to be viable for tsunami simulations (The number of grid points in a computational grid necessary for DNS scales as the (9/4)th power of the Reynolds number.) [107]. Instead, by means of either a Reynolds averaging approach that gives rise to the RANS equations, or scale separation filtering for LES, turbulence can be modeled with sufficient precision at a drastically reduced computational cost with respect to DNS. RANS was first used to simulate tsunami run-up in the software COBRAS [108,109]. However, it took approximately one more decade to become a popular choice to model tsunami generated turbulence. Examples of RANS-based tsunami models are described in [17,18,[110][111][112].
In terms of solution cost, LES lies between DNS and RANS. Due to its computational cost, it is not popular among tsunami modelers, although it was recently considered as an optimal choice to study breaking waves and run-up in highly non-steady dynamics (see, e.g., in [90]). 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 interested in the analysis of tsunami impact and coastal propagation, although it would be unreasonable to use either RANS or LES to model for off-shore propagation where, instead, the 2D shallow water equations do a perfect job.
The free surface in a model based on the solution of (1) is typically accounted for in either one of two ways: by means of a volume of fluid approach (in, e.g., [112,113]) or by a level-set tracking of the water-air interface in a two-fluid model (in, e.g., [90,114]). Both of these approaches make it natural to handle the wet-dry front as the flow reaches land; when high order methods are utilized, this may not be the case for shallow-water equations whose solution is challenged in wet-dry regions for certain numerical approximations, as demonstrated by the continuous work of, for example, [36][37][38][103][104][105].
In summary, all the numerical methods of solution as well as the mathematical models 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 choose one model over another. The most prominent metrics are listed as follows.

•
Flow scale and régime. For example, are we modeling breaking waves in the vicinity of the shore or linear waves propagating in the ocean basin? Is turbulence an important factor? • 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. • Performance on the computing architecture being considered. • Overall size of the problem considered. Is one interested in a single simulation or a large ensemble in order to account for uncertainty?
These are of course only some of the issues but they 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 ongoing research on all fronts from tsunami generation, propagation, to fine-scale run-up, that tsunami modelers may soon join forces to build a unified, all-scale framework to model a full tsunami event from source to impact (see, e.g., in [27,115]). The all-scale idea is naïvely represented in Figure 5 where the tsunami's domain is subdivided in different zones as a function of the flow régime and complexity. The tsunami assumes different shapes during its lifespan from generation to impact; with time scales that range from seconds (100-1000 s at the source [50]) to up to several hours during propagation and run-up/run-down, the flow is driven by both two-dimensional and three-dimensional dynamics across spatial scales from submeters to hundreds of kilometers. Due to such enormous scale variability in time and space, not all régimes can or should be described by the same physical model. While it may be or become technically possible, it is not sensible to solve the 3D Navier-Stokes equations for turbulent flows in the far-field where the shallow water equations are appropriate and inexpensive to solve. The wave propagates in the open ocean with a small amplitude from a few centimeters up to several tens of centimeters (Ref. [50] and citations therein), which is always much smaller than the open ocean depth which is, in turn, much smaller than the wavelength. Under these circumstances, the tsunami motion is well represented by the shallow water equations [116][117][118] for linear, quasi-linear, and quasi-dispersive waves. While the wave approaches land it decelerates, its amplitude increases, and its wavelength is drastically reduced; furthermore, it loses its linear, non-dispersive nature. In these conditions a non-linear model that describes turbulent breaking waves becomes necessary. The most complete model of this is given by the 3D Navier-Stokes equations of incompressible flows. If the only metric of interest during run-up were the inundation level, the 2D shallow water equations would still provide relatively satisfactory results, but their limits become clear when the correct hydrodynamic forces on structures are to be evaluated.
To achieve a complete multi-scale model of the tsunami, different solvers may be coupled as illustrated in such a way that the 3D Navier-Stokes solver is forced at the boundaries of the shore proximity area by the shallow water model of the tsunami moving off-shore [119,120]. At the source the shallow water model is forced by an earthquake simulator or 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., in [121][122][123][124].) 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
This article gives a concise review of the different approaches to tsunami modeling from generation to impact. From the perspective of a future comprehensive multi-scale modeling infrastructure to simulate a full tsunami, the article underlines the challenges associated with this approach and reviews the efforts that are underway to achieve this goal.
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, either by means of grid refinement or high order methods in the open ocean, or by means of a complete three-dimensional model of the flow in the near shore region. 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 earthquaketsunami modeling at the source, adaptive grid refinement, and ever faster and inexpensive hybrid high-performance computing (i.e., graphics processing units) 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 leading the way towards full tsunami forecast as envisioned by Synolakis et al. [91] 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, the advent of exa-scale computing within the next decade and the ever-increasing presence of general purpose GPUs mounted on personal computers by default led us to think that such a tool may be available within the next two decades.

Acknowledgments:
The authors would like to thank two anonymous reviewers whose comments and suggestions helped improve the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
Some of the most common tsunami simulators are listed in Table A1. The table is maintained in an open repository hosted at https://github.com/mandli/tsunami-models to allow developers and scientists to add their models.