Mathematical Models of Tsunami Propagation

1. Introduction

Compared with other types of geohazards, such as hurricanes and earthquakes, tsunamis are relatively infrequent. According to the Global Historical Tsunami Database [1], tsunamis that cause economic damage or deaths near their sources occur approximately twice per year. Tsunamis that cause damage or deaths on distant shores (i.e., over 1,000 km away) occur about twice per decade. However, the impact of a destructive tsunami can be devastating. An example in recorded history was the 2004 Indian Ocean tsunami, which was caused by a magnitude 9.1 earthquake off the coast of Sumatra in Indonesia [2]. It killed at least 225,000 people, mostly in Indonesia, Thailand, India, and Sri Lanka, and caused extensive economic and infrastructure damage in the affected areas.

Because tsunamis are rare events, it is often difficult or impossible to compile sufficient observational evidence to study and forecast tsunami hazards for a specific region. A natural alternative is to use computer models to study past tsunami events and to provide forecasts that can help communities prepare their responses before the next event strikes. With the rapid advances in computing technology, numerical modeling of tsunamis is quickly becoming a powerful tool to further our understanding of the physics of tsunamis and to predict their impacts on coastal communities.

The evolution of a typical tsunami event usually consists of three phases: the generation phase, in which a tsunamigenic source drives the movement of the water column and generate tsunami waves, the propagation phase, in which the tsunami waves travel from the source location across the open ocean, and the inundation phase, which deals with the breaking of the waves as they approach the shores. The inundation phase is usually the most difficult to model, as the dynamics of the waves during this phase are highly sensitive to small changes in the topography/bathymetry and to the patterns of the structural damages caused by the inundation. Early attempts of tsunami numerical modeling started about fifty years ago based upon simplified theories of water waves (e.g., [3,4]). At that time, computer hardware had limited capacities, and the spatial resolution of the simulations was poor (e.g., grid size of 10 km was common even in near-field simulations). With the rapid advances in computer hardware through the years, more rigorous and sophisticated water wave theories have been developed and implemented with much improved numerical accuracy and spatial resolution. In this review, we focus on some of the major mathematical models of tsunami propagation, which have been the foundation of many numerical simulations.

2. Tsunamigenic Sources

Tsunamis can be generated by a variety of sources. Some common sources include tectonic earthquakes (e.g., [5,6]), landslides (e.g., [7,8,9]), submarine volcanic eruptions (e.g., [10,11]) and certain unconventional mechanisms (e.g., [12,13]). Each different type of source mechanism requires a different numerical treatment and the study of numerical simulation techniques for tsunamigenic sources has become an important research topic by itself, especially in near-field tsunami modeling (e.g., [14]). A detailed discussion of tsunamigenic source modeling techniques is beyond the scope of this review. However, we give below a very brief discussion on the mathematical models of two most common types of tsunamigenic sources: tectonic earthquake ruptures and submarine mass failures. We note that in the far field, i.e., over two average wavelengths away from the source, the details of the source mechanisms have less impact on inundation and other important wave characteristics, which opens the possibility of adopting simplified source parameterizations (e.g., [15,16]) for far-field tsunami modeling.

Early attempts in modeling tsunamigenic earthquake ruptures usually assumed simple rectangular rupture geometries and an elastic, homogeneous, half-space model for the Earth’s crust (e.g., [17]). The seafloor deformation caused by the rupture can then be computed analytically from the well-known Okada parameter set [18] using this type of rupture model. Deformations from ruptures with more complex geometries can be approximated through superposition of the analytic solutions from small rectangular patches covering the rupture surface with certain time-dependence (e.g., [19]). Deformations from complex ruptures embedded in a 3D heterogeneous crust model can be computed using the finite-element method (FEM) (e.g., [20,21]). Such FEM rupture models have been integrated into tsunami simulations and provided improved predictions of observations, especially in the near field, when compared with those obtained using semi-analytic rupture models based on Okada’s solution (e.g., [22]). More sophisticated dynamic rupture models based on rate and state friction laws have been developed (e.g., [23,24]). However, the computational cost for these rupture models can be substantial.

Subaerial and submarine mass failures (SMFs) constitute the second most frequent type of tsunamigenic sources (e.g., [25,26,27]). SMF sources can be modelled as rigid sliding blocks by applying laws of motion developed for translational or rotational center-of-mass motions (e.g., [28,29,30]). These models can be integrated into tsunami simulations as bottom boundary conditions [29,31] or as inputs to tsunami generation models (e.g., [30,32,33,34]), such as TOPICS (Tsunami Open and Progressive Initial Conditions System) [33], to construct water wave initial conditions. Aside from the simplicity of these models, an important benefit is that many of the parameters of these models can be determined in (near-)real-time using seismic waveform observations [9], which opens the possibility of more reliable early warnings for SMF-generated tsunamis. More sophisticated SMF source models that can account for complex dynamics and material behavior have also been developed. For example, deformable slides can be modelled as (non-)Newtonian dense viscous fluids (e.g., [35,36]) or as saturated granular debris flow (e.g., [37,38,39]) or as one component in a multi-phase (multi-material) flow (e.g., [40,41]).

3. Wave Propagation

Tsunami propagation modeling can be treated as a special case of the general theory for water waves (e.g., [42]). Mathematical models for water waves can be classified based on whether they study the evolution of individual waves in space-time (i.e., phase-resolving) or the energy of individual frequency components (i.e., phase-averaging). Phase-averaging models do not treat waves individually, but instead simulate the spectral evolution of dispersive waves based on the wave energy balance equation or wave action balance equation (e.g., [43,44,45,46]). In the past three decades, various phase-averaging models, such as SWAN (Simulating Waves Nearshore) [47,48] and WAVEWATCH III [49,50,51], have been developed and widely applied to simulate and predict wind-generated waves in oceans, coastal and inland waters.

Phase-resolving wave models are based on solutions of the conservation laws of mass and momentum for fluid motions. Tsunami waves can, in their most fundamental form, be described by the Navier-Stokes equations for incompressible fluids. Numerical solutions of these equations form the basis of computational fluid dynamics (e.g., [52,53,54,55,56]). For many applications, however, simplified approximate equations can be adopted by neglecting the effects of e.g., viscosity, turbulence, non-linearity.

3.1. Conservation Laws

The fundamental principle of mass conservation of the fluid leads to the continuity equation, which can be expressed as

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathbf{u}\right)=0,$$

(1)

where $\rho$ is the fluid density, $\mathbf{u}$ is the fluid particle velocity, t is time, $\nabla ={\left(\partial /\partial x,\partial /\partial y,\partial /\partial z\right)}^T$ is the spatial gradient operator, with the vertical coordinate z measured upwards from the still-water level $z=0$. The first term on the left-hand-side (LHS) represents the temporal change in fluid density, and the second term represents the mass flow, which can be expanded as,

$$\nabla ·\left(\rho u\right)=\mathbf{u}·\nabla \rho +\rho \nabla ·\mathbf{u}.$$

(2)

Water can be approximately considered to be incompressible, especially under normal conditions. In the deep oceans at 4 km depth, where pressures are 40 MPa ( 395 standard atmospheric pressure), there is only a 1.8% decrease in volume [57]. If we approximate water as incompressible, we obtain the constraint

$$\nabla ·\mathbf{u}=0.$$

(3)

With this constraint, the general continuity equation 1 can be simplified to the incompressible continuity equation

$${\displaystyle \frac{\mathrm{D}\rho }{\mathrm{D}t}}=0,$$

(4)

where the material derivative $\mathrm{D}/\mathrm{D}t$ is defined, for an arbitrary scalar or vector field $\mathbf{A}$, as

$${\displaystyle \frac{\mathrm{D}\mathbf{A}}{\mathrm{D}t}}\equiv {\displaystyle \frac{\partial \mathbf{A}}{\partial t}}+\mathbf{u}·\nabla \mathbf{A}.$$

(5)

The fundamental principle of momentum conservation for incompressible fluids leads to the Navier-Stokes equation

$${\displaystyle \frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t}}=-{\displaystyle \frac{\nabla p}{\rho }}+\nu {\nabla }^2\mathbf{u}+{\displaystyle \frac{\mathbf{F}}{\rho }},$$

(6)

where p is the pressure, $\nu$ is the kinematic viscosity, and $\mathbf{F}$ is the external body force. The LHS of equation 6 is the acceleration of the fluid and the right-hand-side (RHS) is the summation of the density-normalized internal and external forces. The internal force acts on the surfaces of the fluid element by its neighboring fluid elements and can be separated into its normal component (i.e., pressure p) perpendicular to the surfaces and the shear stress parallel to the surfaces. For incompressible fluids, the shear stress can be written in terms of the kinematic viscosity $\nu$. For ocean waves, the external force is usually the gravity force, i.e.,

$$\mathbf{F}=\nabla \left(-\rho gz\right),$$

(7)

where g is the gravitational acceleration.

Under normal conditions, water has a kinematic viscosity of 1,002 × 10^-6 m²/s, and can be approximated as inviscid, in which case the second term on the RHS of equation 6 disappears and the incompressible Navier-Stokes equation reduces to the incompressible Euler equation,

$${\displaystyle \frac{\mathrm{D}u}{\mathrm{D}t}}=-{\displaystyle \frac{\nabla p}{\rho }}-\mathbf{g}$$

(8)

where $\mathbf{g}={\left(0,0,g\right)}^T$. In many applications, including tsunami propagation, equation 4 is often neglected, and we only need to solve equations 3 and 8, assuming that the fluid is both incompressible and inviscid.

3.2. Potential Flow Theory

The Navier-Stokes equation and the Euler equation are vectorial equations written in terms of the primitive variable, the fluid velocity vector $\mathbf{u}$. We can convert them into scalar equations if we assume that the fluid motion is irrotational, i.e., the vorticity $\mathit{\omega }$ is zero,

$$\mathit{\omega }=\nabla \times \mathbf{u}=0.$$

(9)

In this case, we can represent the fluid velocity as the gradient of a scalar function $\varphi$, which is called the velocity potential, i.e.,

$$\mathbf{u}=\nabla \varphi ,$$

(10)

since the curl of the gradient of a scalar function is always zero.

If we represent the velocity vector as $\mathbf{u}={\left(u,v,w\right)}^T={\left(\mathbf{v},w\right)}^T$, where the horizontal velocity vector $\mathbf{v}=\left(u,v\right)$, and split the spatial gradient operator ∇ into its horizontal and vertical components, i.e., $\nabla ={\left(\nabla ,{\partial }_z\right)}^T$, where $\nabla =\left({\partial }_x,{\partial }_y\right)$, and ${\partial }_{\alpha }=\partial /\partial \alpha$, $\alpha =x,y,z$, equation 10 can be written as

$$\mathbf{v}=\nabla \varphi ,w={\partial }_z\varphi .$$

(11)

Bring equation 10 into the continuity constraint (eq. 3), we obtain the Laplace equation for the velocity potential,

$${\nabla }^2\varphi =0,$$

(12)

which can be written in terms of the horizontal gradient operator ∇ as,

$${\nabla }^2\varphi +{\partial }_z^2\varphi =0.$$

(13)

Assuming the fluid is inviscid, and the external force is only gravity, we can bring equation 10 into the Euler equation (eq. 8). After some vector calculus manipulations and integrating over spatial dimensions, we obtain the Bernoulli equation [52],

$${\partial }_t\varphi +{\displaystyle \frac{1}{2}}{\left|\nabla \varphi \right|}^2+{\displaystyle \frac{1}{2}}{\left({\partial }_z\varphi \right)}^2+gz+{\displaystyle \frac{p-p_0}{\rho }}=0,$$

(14)

where $p_0$ is the reference pressure and can be set to the atmospheric pressure for tsunami simulations.

To solve the surface wave problem, we need to impose boundary conditions at the free surface, which can be represented as

$$\Sigma \left(x,y,z,t\right)\equiv \eta \left(x,y,t\right)-z=0$$

(15)

and the fixed solid bottom,

$$z=-h\left(x,y\right).$$

(16)

The shape of the bottom boundary $h\left(x,y\right)$ is usually given as bathymetry data, while the free surface function $\eta \left(x,y,t\right)$ must be found as part of the solution.

At the free surface, we impose a kinematic boundary condition and a dynamic boundary condition. The kinematic boundary condition requires that the fluid particle must always stay in the water, i.e., the vertical velocity of a fluid particle must equal the vertical velocity of the free surface,

$${\partial }_t\eta +\nabla \varphi ·\nabla \eta -{\partial }_z\varphi =0,z=\eta \left(x,y,t\right).$$

(17)

The dynamic boundary condition requires that the normal stresses, which are given by the differences in pressure at the free surface, must be in balance, i.e., the Bernoulli equation 14 at the free surface becomes

$${\partial }_t\varphi +{\displaystyle \frac{1}{2}}{\left|\nabla \varphi \right|}^2+{\displaystyle \frac{1}{2}}{\left({\partial }_z\varphi \right)}^2+g\eta =0,z=\eta \left(x,y,t\right).$$

(18)

At the solid bottom boundary, the fluid is required not to penetrate the boundary, i.e., the fluid velocity in the normal direction of the bottom boundary $\mathbf{n}$ must vanish, $\partial \varphi /\partial \mathbf{n}=0$. This bottom boundary condition leads to

$$\nabla \varphi ·\nabla h+{\partial }_z\varphi =0,z=-h\left(x,y\right).$$

(19)

The initial value problem specified by equations 13, 17, 18, 19 can be solved given the initial fields $\eta \left(x,y,0\right)$ and $\varphi \left(x,y,z,0\right)$ at time $t=0$. Once the fields $\eta$ and $\varphi$ have been found, the momentum conservation equation, the Euler equation 8, can be used a posteriori to solve for the pressure p.

The boundary conditions at the free surface given by equations 17, 18 are nonlinear. However, the nonlinear terms can be neglected if we assume that the waves are non-steep, i.e., $H\ll \lambda$, where H is the wave height and $\lambda$ is the wavelength. The non-steep assumption is often represented using a small parameter

$$\alpha =ka\ll 1,$$

(20)

where $k=2\pi /\lambda$ is the wave number and $a=H/2$ is the wave amplitude. With this assumption, the fields $\varphi$ and $\eta$ can be expanded as perturbations around the still-water free surface $z=0$, and the second and higher-order terms can be dropped. The nonlinear terms in equations 17 and 18 can be neglected, and the resulting linear boundary conditions are

$${\partial }_t\eta ={\partial }_z\varphi ,z=0,$$

(21)

$${\partial }_t\varphi +g\eta =0,z=0.$$

(22)

Taking the time partial derivative on both sides of equation 22 and combining with equation 21, we obtain

$${\displaystyle \frac{{\partial }^2\varphi }{\partial t^2}}+g{\displaystyle \frac{\partial \varphi }{\partial z}}=0,z=0.$$

(23)

When combined with the bottom boundary condition equation 19 and the Laplace equation 13, we obtain the linear potential flow theory, which is also known as the Airy wave theory (e.g., [58,59,60] ).

An analytic solution to the linear potential flow theory can be found through separation of variables. If we consider only one horizontal dimension x and the vertical dimension z, an ansatz for a harmonic wave propagating along the x axis can be expressed as

$$\varphi \left(x,z,t\right)=Z\left(z\right)sin\left(kx-\omega t\right),$$

(24)

where $\omega$ is the angular frequency of the wave. Requiring the ansatz to satisfy the Laplace equation 13 gives

$$Z\left(z\right)=A{\mathrm{e}}^{kz}+B{\mathrm{e}}^{-kz}$$

(25)

for some constants A and B. If we assume that the bottom boundary $h\left(x,y\right)$ varies slowly, the first term in the LHS of equation 19 can be neglected, and inserting the ansatz gives

$$A/B={\mathrm{e}}^{2kh}.$$

(26)

Combining equations 25 and 26, we obtain

$$Z\left(z\right)=Ccosh\left(k\left(z+h\right)\right),$$

(27)

where C is a new constant and is often chosen such that $Z\left(0\right)=gH/\left(2\omega \right)$. The fluid potential is obtained as

$$\varphi \left(x,z,t\right)={\displaystyle \frac{gHcosh\left(k\left(z+h\right)\right)}{2\omega cosh\left(kh\right)}}sin\left(kx-\omega t\right).$$

(28)

Finally, the ansatz must satisfy equation 23, which results in

$${\omega }^2=gktanh\left(kh\right).$$

(29)

The dispersion equation 29 relates the frequency of the wave with its wavelength and implies that the phase speed, also known as celerity, $c_p=\omega /k$, is different for waves with different wavelengths. The group speed, $c_g=\mathrm{d}\omega /\mathrm{d}k$, which defines the speed of energy propagation, is also a function of the wavelength. The surface elevation can be obtained by bringing equation 28 into equation 22,

$$\eta \left(x,t\right)={\displaystyle \frac{H}{2}}cos\left(kx-\omega t\right),$$

(30)

which is a harmonic plane wave propagating along the x axis. The linearity of the linear potential flow theory allows us to construct new solutions through linear superposition of the solutions in the form of equation 28. In general, ocean waves are composed of many such harmonic components travelling in different directions.

The dispersion relation in equation 29 can be re-written as,

$$c_p^2={\displaystyle \frac{{\omega }^2}{k^2}}=gh{\displaystyle \frac{tanh\left(kh\right)}{kh}},$$

(31)

and we have the Taylor expansion

$${\displaystyle \frac{tanh\left(kh\right)}{kh}}=1-{\displaystyle \frac{{\left(kh\right)}^2}{3}}+{\displaystyle \frac{2{\left(kh\right)}^4}{15}}-{\displaystyle \frac{17{\left(kh\right)}^6}{315}}+{\displaystyle \frac{62{\left(kh\right)}^8}{2835}}+O\left({\left(kh\right)}^{10}\right).$$

(32)

In the long-wave limit, $h\ll \lambda$, the relative depth $kh\to 0$, $tanh\left(kh\right)/\left(kh\right)\to 1$, we have $c_p\to \sqrt{gh}$ and the group speed is no longer a function of the wavenumber, i.e., the waves are nondispersive. Beyond the leading order term in the Taylor expansion, the wavenumber dependency is retained. The leading order frequency dispersion effect can be accounted for by including the $O\left({\left(kh\right)}^2\right)$ term, which can extend the applicable relative depth range of the long wave theory to include weakly dispersive waves (e.g., [60]).

3.3. Long Wave Scaling

As shown in the potential flow theory, reasonable assumptions can be made to simplify the Navier-Stokes equations, perhaps substantially, to provide useful descriptions of the wave behavior. Among the various assumptions that have been adopted when developing mathematical models of tsunami propagation, there are three important length scales: the characteristic water depth $h_0$, the characteristic wavelength in the horizontal direction ${\lambda }_0$, and the characteristic wave amplitude $a_0$. We can introduce two dimensionless parameters, which are often assumed to be small for tsunami propagation modeling in the open ocean,

$$ϵ={\displaystyle \frac{a_0}{h_0}},$$

(33)

$$\mu ={\displaystyle \frac{h_0}{{\lambda }_0}}.$$

(34)

The scale factor $ϵ$ is often used to measure the nonlinearity of the problem and the scale factor $\mu$ defines the relative depth and ${\mu }^2$ often serves as a measure of frequency dispersion (e.g., [61,62]).

A set of dimensionless independent and dependent variables can be defined in terms of the characteristic length scales [61] and the equations written in these dimensionless variables are more amenable to scale analysis (or order-of-magnitude analysis). The new dimensionless independent variables are defined as,

$$\tilde{x}={\displaystyle \frac{x}{{\lambda }_0}},\tilde{y}={\displaystyle \frac{y}{{\lambda }_0}},\tilde{z}={\displaystyle \frac{z}{h_0}},\tilde{t}={\displaystyle \frac{c_0}{{\lambda }_0}}t,$$

(35)

where $c_0=\sqrt{g{h}_0}$ is the characteristic speed of tsunami propagation in the open ocean, ranging from 356 km/h for 1-km water depth to 712 km/h for 4-km water depth. The new dimensionless dependent variables are defined as [61,62],

$$\tilde{u}={\displaystyle \frac{u}{ϵc_0}},\tilde{v}={\displaystyle \frac{v}{ϵc_0}},\tilde{w}={\displaystyle \frac{\mu w}{ϵc_0}},\tilde{\eta }={\displaystyle \frac{\eta }{a_0}},\tilde{\varphi }={\displaystyle \frac{\varphi }{ϵc_0{\lambda }_0}},\tilde{h}={\displaystyle \frac{h}{h_0}},\tilde{p}={\displaystyle \frac{p}{\rho g{a}_0}}.$$

(36)

With these definitions of the dimensionless variables, the bottom boundary and the free surface boundary are given by

$$\tilde{z}=-\tilde{h},\tilde{z}=ϵ\tilde{\eta },$$

(37)

respectively. The partial derivatives with respect to the spatial and temporal coordinates can be obtained through the chain rule, i.e., ${\partial }_{\alpha }={\displaystyle \frac{\partial \tilde{\alpha }}{\partial \alpha }}{\partial }_{\tilde{\alpha }}$ for $\alpha =x,y,z,t$.

The governing equations can be written in terms of the dimensionless variables [61]. In the following, for the sake of convenience, we drop the tilde on the dimensionless variables. The continuity constraint, equation 3, can be written in the dimensionless variables as

$${\mu }^2\left({\partial }_{x}u+{\partial }_{y}v\right)+{\partial }_{z}w=0.$$

(38)

The irrotational constraint, equation 9, is transformed into,

$${\partial }_{y}u={\partial }_{x}v,{\partial }_{z}v={\partial }_{y}w,{\partial }_{x}w={\partial }_{z}u.$$

(39)

The momentum conservation law, the Euler equation 8, can be written as

$${\mu }^2{\partial }_{t}u+ϵ{\mu }^2\left(u{\partial }_{x}u+v{\partial }_{y}u\right)+ϵw{\partial }_{z}u+{\mu }^2{\partial }_{x}p=0,$$

(40)

$${\mu }^2{\partial }_{t}v+ϵ{\mu }^2\left(u{\partial }_{x}v+v{\partial }_{y}v\right)+ϵw{\partial }_{z}v+{\mu }^2{\partial }_{y}p=0,$$

(41)

$$ϵ{\partial }_{t}w+{ϵ}^2\left(u{\partial }_{x}w+v{\partial }_{y}w\right)+{\left({\displaystyle \frac{ϵ}{\mu }}\right)}^{2}w{\partial }_{z}w+ϵ{\partial }_{z}p+1=0.$$

(42)

The boundary condition at the fixed solid bottom, equation 19, is converted into

$$w+{\mu }^2\left(u{\partial }_{x}h+v{\partial }_{y}h\right)=0,z=-h.$$

(43)

The kinematic boundary condition, equation 17, and the dynamic boundary condition at the free surface, are expressed as

$$w={\mu }^2{\partial }_t\eta +ϵ{\mu }^2\left(u{\partial }_x\eta +v{\partial }_y\eta \right),p=0,z=ϵ\eta .$$

(44)

The governing equations written in the dimensionless variables contain the two scaling parameters $ϵ$ and $\mu$, which facilitates scale analysis of the contributions of each term to the nonlinearity and dispersion of the flow problem.

3.4. Depth Integration

It can be computationally demanding to solve the nonlinear, three-dimensional, initial boundary value problem specified in equations 38-44 for tsunami propagation spanning many wavelengths. Depth integrating the governing equations can reduce the three-dimensional problem to a two-dimensional problem, leading to substantial savings in computational cost. Such quasi-3D approaches do not seek to resolve the exact depth-dependence of the fields but allow different assumptions about the fields’ depth profiles to be adopted based on the practical application (e.g., [61]).

The nondimensional continuity constraint, equation 38, can be integrated vertically from the solid bottom to the free surface. By applying the Leibniz integral rule and the boundary conditions in equations 43, 44, we obtain

$${\partial }_t\eta +{\partial }_x{\int }_{-h}^{ϵ\eta }u\mathrm{d}z+{\partial }_y{\int }_{-h}^{ϵ\eta }v\mathrm{d}z=0.$$

(45)

The pressure field can be obtained by integrating the vertical component of the momentum conservation law, equation 42, from an arbitrary depth z inside the flow to the free surface. After applying the Leibniz integral rule and the boundary conditions at the free surface, equation 44, we obtain,

$$p=\eta -{\displaystyle \frac{z}{ϵ}}+{\partial }_t{\int }_z^{ϵ\eta }w\mathrm{d}z+ϵ{\partial }_x{\int }_z^{ϵ\eta }uw\mathrm{d}z+ϵ{\partial }_y{\int }_z^{ϵ\eta }vw\mathrm{d}z-{\displaystyle \frac{ϵ}{{\mu }^2}}{w}^2.$$

(46)

The two horizontal components of the momentum conservation law, equations 40 and 41 are integrated from the solid bottom to the free surface. After applying the Leibniz rule, the continuity constraint (eq. 38) and the boundary conditions (eq. 43-44), we obtain,

$${\partial }_t{\int }_{-h}^{ϵ\eta }u\mathrm{d}z+ϵ{\partial }_x{\int }_{-h}^{ϵ\eta }{u}^2\mathrm{d}z+ϵ{\partial }_y{\int }_{-h}^{ϵ\eta }uv\mathrm{d}z+{\partial }_x{\int }_{-h}^{ϵ\eta }{p\mathrm{d}z-p|}_{z=-h}{\partial }_{x}h=0,$$

(47)

$${\partial }_t{\int }_{-h}^{ϵ\eta }v\mathrm{d}z+ϵ{\partial }_x{\int }_{-h}^{ϵ\eta }uv\mathrm{d}z+ϵ{\partial }_y{\int }_{-h}^{ϵ\eta }{v}^2\mathrm{d}z+{\partial }_y{\int }_{-h}^{ϵ\eta }{p\mathrm{d}z-p|}_{z=-h}{\partial }_{y}h=0.$$

(48)

The vertical velocity w can be obtained by integrating the continuity constraint (eq. 38) from the solid bottom to an arbitrary depth z,

$$w=-{\mu }^2\left({\partial }_x{\int }_{-h}^{z}u\mathrm{d}z+{\partial }_y{\int }_{-h}^{z}v\mathrm{d}z\right).$$

(49)

The depth integrated governing equations 45-49 only contain differential operators with respect to the horizontal coordinates ${\partial }_x,{\partial }_y$ and time ${\partial }_t$. All the depth integrals are evaluated exactly, and no approximations have been made so far.

3.5. Shallow-Water Equations

To obtain the shallow-water equations, which have been highly successful for tsunami propagation modeling, we make the long-wave approximation, i.e., the wavelength is much larger than the depth of the flow ($\mu \ll 1$). In such a regime, the vertical components of the velocity, w, and the acceleration, ${\partial }_{t}w$, can be neglected, and the expression for dimensionless pressure, eq 46, can then be simplified as

$$p=\eta -{\displaystyle \frac{z}{ϵ}}.$$

(50)

The horizontal components of the velocity are assumed to be invariant with depth, therefore the continuity constraint, eq 45, can be expressed as

$${\partial }_t\eta +{\partial }_x\left[\left(ϵ\eta +h\right)u\right]+{\partial }_y\left[\left(ϵ\eta +h\right)v\right]=0.$$

(51)

Bringing the hydrostatic pressure, eq 50, into the two horizontal components of the momentum equations 40, 41, we obtain

$${\partial }_{t}u+ϵ\left(u{\partial }_{x}u+v{\partial }_{y}u\right)+{\partial }_x\eta =0,$$

(52)

$${\partial }_{t}v+ϵ\left(u{\partial }_{x}v+v{\partial }_{y}v\right)+{\partial }_y\eta =0.$$

(53)

Equations 51-53 are the nonlinear shallow-water equations (NSWE) written in dimensionless variables.

The NSWE can be linearized if we consider small disturbances about a fluid at rest, i.e., $ϵ\ll 1$, and the horizontal velocity components u and v are negligible. In this case, the nonlinear advection terms (i.e., second terms on the LHS) of equations 52-53 can be dropped

$${\partial }_{t}u+{\partial }_x\eta =0,$$

(54)

$${\partial }_{t}v+{\partial }_y\eta =0.$$

(55)

And the continuity constraint, eq 51, can be simplified to

$${\partial }_t\eta +{\partial }_x\left(hu\right)+{\partial }_y\left(hv\right)=0.$$

(56)

Equations 54-56 are the linear shallow-water equations (LSWE) written in dimensionless variables.

Tsunami modeling based on LSWE started about five decades ago (e.g., [3,4]). The finite-difference method (FDM) was adopted to discretize the LSWE and the simulated tsunami waveforms at tide gauges and the predicted coastal tsunami height distributions were generally consistent with the observations [3]. Such tsunami propagation simulations based on LSWE have been widely adopted in practice for research and operational purposes (e.g., [63,64,65]).

In coastal zones, where wave amplitudes increase due to shoaling, the assumption of linearity becomes invalid and the NSWE should be adopted (e.g., [63,64,65]). Numerical solutions of NSWE based on the FDM have been implemented and applied extensively (e.g., [66,67,68]). Several operational codes (e.g., [64,65,66]) allow for solving either the LSWE or the NSWE by just changing one of the input parameters. The FDM implementations use uniform or nested grids. In a nested grid, coarser spatial resolution is used in the open ocean and the shallow sea, while finer spatial resolution is used for inundation and runup on the land. Alternatively, more flexible, unstructured grids can be utilized in the finite-element method (FEM) (e.g., [69,70]) and the finite volume method (FVM) (e.g., [71,72]) to provide seamless resolution adaptivity for resolving tsunami wave behavior at different spatial scales.

3.6. Boussinesq Equations

Under the long-wave approximation adopted in the NSWE and LSWE, the speed of wave propagation depends only on the water depth. However, tsunami waves may have a wide range of frequencies caused by the mechanisms of the tsunamigenic sources (e.g., [73]), and each frequency component propagates at a different speed (i.e., frequency dispersion). At high frequencies (i.e., short wavelengths), the long-wave approximation may become invalid (e.g., [74]) and Boussinesq equations can provide better approximations than the shallow-water equations. Past experiences have shown that frequency dispersion has significant effects on tsunami waves generated by SMFs [26] and dispersive models should be used for simulating such tsunamis whenever possible.

The classic Boussinesq equations (e.g., [75]) were derived using the depth-averaged horizontal velocities, i.e.,

$$\overline{\mathbf{v}}\left(x,y,t\right)={\displaystyle \frac{1}{ϵ\eta +h}}{\int }_{-h}^{ϵ\eta }\mathbf{v}\left(x,y,z,t\right)\mathrm{d}z,$$

(57)

as the model velocities. Bringing equation 57 into the depth-integrated continuity constraint, eq 45, we obtain

$${\partial }_t\eta +\nabla ·\left[\left(ϵ\eta +h\right)\overline{\mathbf{v}}\right]=0.$$

(58)

Substituting equation 57 into the depth-integrated horizontal moment equations 47-48 and integrating, retaining terms up to $O\left(ϵ\right)$ and $O\left({\mu }^2\right)$, we obtain

$${\partial }_t\overline{\mathbf{v}}+\nabla \eta +ϵ\left(\overline{\mathbf{v}}·\nabla \right)\overline{\mathbf{v}}+{\mu }^2\left\{{\displaystyle \frac{h^2}{6}}\nabla \left(\nabla ·{\partial }_t\overline{\mathbf{v}}\right)-{\displaystyle \frac{h}{2}}\nabla \left[\nabla ·\left(h{\partial }_t\overline{\mathbf{v}}\right)\right]\right\}=0.$$

(59)

Compared with the NWSE momentum equations 52-53, the classic Boussinesq momentum equation 59 contains third-order derivatives with respect to both spatial and temporal coordinates, which are responsible for the dispersion effects.

A major drawback of the classic Boussinesq equations is the limited range of applicability. The relative depth $\mu$ must be less than about 0.2 (i.e., shallow-water waves) to keep the error in phase speed below about 5% [61,76]. To enhance deep-water accuracy, it is better to adopt the generalized Boussinesq equations (e.g., [61]). Instead of using the depth-averaged horizontal velocities, equation 57, as the model velocities, the horizontal velocities ${\mathbf{v}}_{\alpha }$ at an arbitrary depth below the surface, $z=z_{\alpha }\left(x,y\right)$, are used as model velocities in generalized Boussinesq equations. The resulting continuity and momentum equations are [61],

$${\partial }_t\eta +\nabla ·\left[\left(h+ϵ\eta \right){\mathbf{v}}_{\alpha }\right]+{\mu }^2\nabla ·\left\{\left({\displaystyle \frac{z_{\alpha }^2}{2}}-{\displaystyle \frac{h^2}{6}}\right)h\nabla \left(\nabla ·{\mathbf{v}}_{\alpha }\right)+\left(z_{\alpha }+{\displaystyle \frac{h}{2}}\right)h\nabla \left[\nabla ·\left(h{\mathbf{v}}_{\alpha }\right)\right]\right\}=0,$$

(60)

$${\partial }_t{\mathbf{v}}_{\alpha }+\nabla \eta +ϵ\left({\mathbf{v}}_{\alpha }·\nabla \right){\mathbf{v}}_{\alpha }+{\mu }^2\left\{{\displaystyle \frac{z_{\alpha }^2}{2}}\nabla \left(\nabla ·{\partial }_t{\mathbf{v}}_{\alpha }\right)+z_{\alpha }\nabla \left[\nabla ·\left(h{\partial }_t{\mathbf{v}}_{\alpha }\right)\right]\right\}=0.$$

(61)

Compared with the continuity equation 58 in classic Boussinesq equations, the continuity equation 60 in generalized Boussinesq equations has additional dispersion terms. In practice, the reference depth $z_{\alpha }$ can be chosen such that the resulting dispersion matches the Padé [2, 2] expansion of the dispersion in the Airy wave theory.

Deep-water accuracy can be further enhanced by developing high-order Boussinesq-type models (e.g., [77,78,79,80,81]). An alternative approach is to divide the water column into several layers, use low-order polynomial approximations for the velocity profiles within each layer and ensure that the velocities are continuous on the interfaces separating neighboring layers (e.g., [82,83]).

3.7. Numerical Solutions of Incompressible Navier-Stokes Equations

The Navier-Stokes (NS) equations are based on the conservation laws of mass and momentum. Unlike the shallow-water or Boussinesq equations, NS equations do not assume a priori any vertical distributions of the velocity fields and are capable of modeling complex wave behavior, such as wave breaking, turbulences, wave-structure interaction. The computational cost for direct numerical solutions of the nonlinear, three-dimensional, NS equations can be substantially higher than that for the shallow-water or the Boussinesq equations. At the current stage, practical applications of direct numerical solutions of NS equations are limited to small-scale problems (e.g., [40,84,85,86,87]).

Solution techniques for incompressible NS equations can generally be classified into Lagrangian, Eulerian and arbitrary Lagrangian Eulerian (ALE) methods. Lagrangian methods are based on tracking the movements of fluid particles through space and time and have been implemented using numerical schemes such as the FEM (e.g., [88,89,90]), the particle FEM (e.g., [91,92]), the moving particle semi-implicit method (e.g., [93,94]) and the smoothed particle hydrodynamics (SPH) (e.g., [95,96,97]). Eulerian methods are based on tracking fluid property changes over time at specific spatial locations and have been implemented using numerical schemes such as the FDM (e.g., [98]) and the FVM ([87,99]). Lagrangian methods are superior in tracking the movements of the free surface or other material interfaces and boundaries. However, large deformations can cause severe mesh distortion and may require periodic re-meshing in grid-based methods. In Eulerian methods, the computational mesh is usually fixed in space, but the positions of the free surface are generally unknown and need to be solved explicitly from an advection equation (e.g., [100,101]). The ALE methods have been developed to combine the advantages of the Lagrangian and Eulerian methods (e.g., [102,103,104]).

3.7.1. Projection Method

A major difficulty in solving the incompressible NS equations using the Lagrangian, Eulerian or ALE methods is the calculation of the pressure field. Unlike compressible NS equations, there is no evolution equation or the equation of state in incompressible NS equations, and the pressure and velocity fields are coupled implicitly through the incompressibility constraint, equation 3. One of the most widely used approaches for coping with this difficulty is the projection method [98], which is based on the Hodge decomposition of a vector field (e.g., [105]). The original projection method is a non-incremental pressure-correction scheme [106]. In the first step, an intermediate velocity ${\mathbf{u}}^{*}$ is computed from the momentum conservation, equation 6, by ignoring the pressure term,

$${\displaystyle \frac{{\mathbf{u}}^{*}-{\mathbf{u}}^n}{\Delta t}}=-\left({\mathbf{u}}^n·\nabla \right){\mathbf{u}}^n+\nu {\nabla }^2{\mathbf{u}}^n+{\displaystyle \frac{{\mathbf{F}}^n}{\rho }},$$

(62)

where superscript n represents the $n^{\mathrm{th}}$ time step and $\Delta t$ is the time step length. In the second step, the intermediate velocity is corrected to obtain the velocity field at the next time step using the pressure field,

$${\displaystyle \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^{*}}{\Delta t}}=-{\displaystyle \frac{\nabla p^{n+1}}{\rho }},$$

(63)

where the pressure field $p^{n+1}$ is obtained by solving the Poisson equation,

$${\nabla }^2{p}^{n+1}={\displaystyle \frac{\rho }{\Delta t}}\nabla ·{\mathbf{u}}^{*},$$

(64)

which results from taking the divergence on both sides of equation 63 and imposing $\nabla ·{\mathbf{u}}^{n+1}=0$. Equation 63 can be rewritten as

$${\mathbf{u}}^{*}={\mathbf{u}}^{n+1}+{\displaystyle \frac{\Delta t}{\rho }}\nabla p^{n+1},$$

(65)

which is the Hodge decomposition of the vector field ${\mathbf{u}}^{*}$ given the boundary condition $\partial p^{n+1}/\partial \mathbf{n}=0$, where $\mathbf{n}$ is the normal direction of the domain boundary. This artificial Neumann boundary condition induces a numerical boundary layer that limits the accuracy of the scheme (e.g., [106,107]). A remedy is to use staggered grids (e.g., [108,109]) and store the pressure and velocity components at different locations. Over the past few decades, a large amount of work has been dedicated to the construction, analysis, and application of projection-type schemes for incompressible NS equations (e.g., [106]) and such schemes have been implemented in tsunami simulations codes, such as NHWAVE [87,110].

3.7.2. Direct Solution Based on Weighted Least-Squares Method

The projection method is an operator splitting approach, in which the velocity and pressure fields are computed separately. An alternative approach is to solve directly for both velocity and pressure fields together at each time step [111]. The momentum conservation in equation 6 can be discretized in time and re-organized as

$${\mathbf{u}}^{n+1}+{\displaystyle \frac{\Delta t}{\rho }}\nabla p^n={\mathbf{u}}^n-\left({\mathbf{u}}^n·\nabla \right){\mathbf{u}}^n+\Delta t\nu {\nabla }^2{\mathbf{u}}^n+{\displaystyle \frac{\Delta t}{\rho }}{\mathbf{F}}^n.$$

(66)

Together with the continuity constraint at time step $n+1$,

$$\nabla ·{\mathbf{u}}^{n+1}=0,$$

(67)

and associated boundary conditions, we can solve directly for ${\mathbf{u}}^{n+1}$ and $p^n$at each time step. In practice, equations 66-67 and the associated boundary conditions are discretized in space, and we obtain an over-determined linear system of equations, which can be solved using the least-squares method [111,112].

Compared with the projection method, this direct solution approach removes pressure corrections that may corrupt boundary conditions because of the underlying Helmholtz-Hodge theorem [105,106] and any spurious pressure oscillations in the numerical solutions due to the decoupling between the velocity field and the pressure field [111]. However, the over-determined linear system resulting from discretizing equations 66-67, together with the associated boundary conditions, can be substantially (about 4 times for 3D problems or 3 times for 2D problems) larger than the linear system resulting from discretizing the Poisson equation 64 in the projection method. To improve the convergence of the direct solution approach, optimal weights can be applied to the discretized continuity constraint, equation 67, and the associated Dirichlet and Neumann boundary conditions, during the least-squares solution of the over-determined linear system [111]. To derive these optimal weights, the L2-norm of the errors associated with the continuity constraint, equation 67, and the Dirichlet and Neumann boundary conditions are required to balance with the L2-norm of the error associated with the momentum conservation equation 66 [111].

3.7.3. Free Surface Capture

Unlike the shallow-water equations or the Boussinesq equations, in which the free surface elevation $\eta \left(x,y,t\right)$ is explicitly solved for, the incompressible NS equations do not have $\eta$ as a dependent variable. The moving free surface needs to be determined indirectly during the solution process of the incompressible NS equations, both for resolving complex topological features, and for enforcing proper boundary conditions. One approach that has been successfully applied in tsunami simulations is to solve the free surface kinematic boundary condition, equation 17, or its depth-integrated form, equation 45, together with the incompressible NS equations in boundary-fitted curvilinear coordinates (e.g., [113]).

Several other types of numerical techniques that have been developed in the computational fluid dynamics community include the boundary-integral method (e.g., [114,115,116,117]), the marker-particle-based methods (e.g., [108,118,119]), the volume-of-fluid (VOF) methods (e.g., [100,120]) and the level-set methods (e.g., [121,122]).

The boundary-integral method is based on the solution of the Laplace equation 12, the kinematic boundary condition, equation 17, and the Bernoulli’s equation 18. The fluid velocity $\nabla \varphi$ needs to be evaluated at the free surface. The evaluation of the tangential derivatives is straightforward. The normal derivative is related to the values of $\varphi$ on the free surface by applying Green’s third identity (e.g., [114]), which provides a Fredholm integral equation of the first kind. An alternative approach is to use the dipole representation (e.g., [115]), which results in a Fredholm integral equation of the second kind. The boundary-integral method reduces the spatial dimension of the flow problem by one, thereby avoiding differentiating dependent variables across discontinuities. But it suffers from difficulties in topology changes of the free surface (e.g., pinching and merging).

Surface marker methods are closely related to the front-tracking methods, which are based on the Lagrangian representation of maker particles placed on the boundary surface and advected by the flow (e.g., [123,124,125]). The marker methods avoid the geometric complexities of reconnecting all the particles as in the front-tracking methods, at the cost of lacking a well-defined, smooth surface. The surface can form gaps and holes due to uneven spatial distributions of marker particles evolving through time. Artificial smoothing near the surface may become necessary to obtain regularized geometric properties (e.g., [126]). The application of velocity boundary conditions requires careful consideration to avoid introducing artificial asymmetries that can cause the simulation to break up [127].

The free surface can also be represented using a fixed grid (i.e., the Eulerian representation). Both the VOF methods and the level-set methods are fixed-grid methods and can support highly complex topologies. VOF methods are based on the definition of a scalar field, $C_{ijk}$, known as the volume fraction or color function,

$$C_{ijk}{V}_{ijk}\approx \int \int {\int }_{ijk}\chi \left(x,y,z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z,$$

(68)

where $\chi$ is the fluid characteristic function, which equals to 1 in the fluid and 0 otherwise, $V_{ijk}$ is the volume of the spatial cell indexed at grid coordinate $\left(i,j,k\right)$ and the volume integration is over the same spatial cell. We have $0<C_{ijk}<1$ in cells cut by the surface, and $C_{ijk}=0$ or 1 for cells away from the surface. For incompressible flow, mass conservation is equivalent to the conservation of the volume and the fluid characteristic function $\chi$, which is also advected by the velocity field $\mathbf{u}$. To reconstruct the surface, we need to find an approximation to the portion of the surface in each cut cell from the volume fraction in that cell and its neighboring cells. The reconstructed surface may not be smooth or continuous, which can reduce the accuracy of the geometric information on the boundary, leading to compromises in solution accuracy eventually. Improvements in the reconstruction process have been developed (e.g., [128,129,130,131]).

The level-set method describes the free surface as the zero level-set (i.e., isocontour) of a scalar function $\xi \left(\mathbf{x},t\right)$, i.e.,

$$\Gamma \left(t\right)=\left\{\mathbf{x}:\xi \left(\mathbf{x},t\right)=0\right\}.$$

(69)

Evolving the free surface dynamically is reduced to evolving the scalar function $\xi$, which is advected by the flow field $\mathbf{u}$, i.e.,

$${\partial }_t\xi +\mathbf{u}·\nabla \xi =0.$$

(70)

The free surface may develop topological singularities (e.g., cusps, pinching, merging), while the level-set function $\xi$ remains relatively smooth (i.e., Lipschitz continuous). Because the level-set method is Eulerian, it is unnecessary to update the mesh to explicitly track the free surface during the flow calculation, the free surface is recovered only at the end of each step by locating the zero level-set. A popular choice for the level-set function is the signed distance towards the surface, i.e.,

$$\xi \left(\mathbf{x}\right)=\pm \left|\right|\mathrm{cp}\left(\mathbf{x}\right)-{\mathbf{x}\left|\right|}_2,$$

(71)

where $\mathrm{cp}(·)$ is the closest-point function providing the coordinates of the point on the surface closest to the point x measured in L2 distance. The sign in equation 71 is chosen as positive if outside and negative if inside a closed surface. This choice of the level-set function simplifies the calculation of geometric properties such as surface normal directions and local mean curvature. The level-set method has been applied successfully to solve a variety of problems in fluid mechanics, combustion, material science (e.g., [132,133,134,135]). However, this method also has its own drawbacks. For example, the signed distance property of the level-set function can be compromised after a few iterations of the advection step in equation 70 [122]. The level-set method does not conserve mass in regions under-resolved by the grid [130,136]. Remedies to these problems have been developed (e.g., [122,137,138,139]).

4. Conclusion

Mathematical models for tsunami propagation have been consistently improving through the past few decades. The models have evolved from simple descriptions of the major characteristics of regular waves to increasingly realistic, sophisticated and robust theories for accurately representing wave behaviors across many scales involved in tsunami propagation and inundation. A direct consequence of the theoretical development is the availability of numerical models implemented in a large collection of community and commercial tsunami simulation codes, which have been widely used for research and operational purposes.

An important challenge in developing new tsunami models is to couple the different scales in a seamless way, from the large scales needed to quantify tsunami propagation to the fine scales needed for representing the breaking of tsunami waves and the subsequent flooding of structures in the exposed areas. While more detailed information can be obtained through more sophisticated simulations, such information needs to be delivered in a timely manner to benefit time-sensitive applications, such as tsunami early warning. Such demanding goals can only be achieved through a combination of judicious choices of the most important features to model and the continued development of efficient computer codes that can take full advantage of the ever-increasing computational resources.