Should tsunami simulations include a nonzero initial horizontal velocity?

Abstract

Tsunami propagation in the open ocean is most commonly modeled by solving the shallow water wave equations. These equations require initial conditions on sea surface height and depth-averaged horizontal particle velocity or, equivalently, horizontal momentum. While most modelers assume that initial velocity is zero, Y.T. Song and collaborators have argued for nonzero initial velocity, claiming that horizontal displacement of a sloping seafloor imparts significant horizontal momentum to the ocean. They show examples in which this effect increases the resulting tsunami height by a factor of two or more relative to models in which initial velocity is zero. We test this claim with a “full-physics” integrated dynamic rupture and tsunami model that couples the elastic response of the Earth to the linearized acoustic-gravitational response of a compressible ocean with gravity; the model self-consistently accounts for seismic waves in the solid Earth, acoustic waves in the ocean, and tsunamis (with dispersion at short wavelengths). Full-physics simulations of subduction zone megathrust ruptures and tsunamis in geometries with a sloping seafloor confirm that substantial horizontal momentum is imparted to the ocean. However, almost all of that initial momentum is carried away by ocean acoustic waves, with negligible momentum imparted to the tsunami. We also compare tsunami propagation in each simulation to that predicted by an equivalent shallow water wave simulation with varying assumptions regarding initial velocity. We find that the initial horizontal velocity conditions proposed by Song and collaborators consistently overestimate the tsunami amplitude and predict an inconsistent wave profile. Finally, we determine tsunami initial conditions that are rigorously consistent with our full-physics simulations by isolating the tsunami waves from ocean acoustic and seismic waves at some final time, and backpropagating the tsunami waves to their initial state by solving the adjoint problem. The resulting initial conditions have negligible horizontal velocity.

Introduction

Tsunamis from earthquake sources pose a significant threat to coastlines around the world and have killed hundreds of thousands of people in the last few decades. The Earth science community often relies on numerical models to better understand tsunami physics (e.g., Titov and Synolakis 1998; Tanioka and Satake 1996; Kowalik 2003), perform hazard assessments (e.g., Geist and Parsons 2006; Borrero et al. 2006; González et al. 2009), and plan for early warning (e.g., Titov et al. 2005; Melgar and Bock 2013; Behrens et al. 2010).

For numerical modelers, the tsunami life cycle is often divided into three parts: generation, propagation, and inundation (Bernard et al. 2006; Mori et al. 2011; Titov and Gonzalez 1997; Watts et al. 2003; Song and Han 2011). This paper is entirely focused on the first two parts, though modeling inundation poses its own unique set of challenges and is far from being completely understood.

It is fairly straightforward to model tsunami propagation in the open ocean, as reviewed by Satake (2002). The simplest and most common way to do so is by using the linearized shallow water wave equations, based on assumptions of inviscid, incompressible water; small amplitude waves; and horizontal wavelengths greater than the ocean depth. In one horizontal dimension, the equations are

$$\frac{\partial \eta }{\partial t} + \frac{\partial (Hv)}{\partial x}= 0,$$

(1)

$$\frac{\partial v}{\partial t} + g\frac{\partial \eta }{\partial x}= 0.$$

(2)

Here, x and y the horizontal and vertical spatial coordinates, respectively, t is time, $h = h(x,t)$ is the depth of the ocean, $v = v(x,t)$ is the depth-averaged horizontal velocity, and g is the gravitational acceleration. As in Fig. 1, we define the sea surface perturbation about $y = 0$ as $y = \eta (x,t)$ and the position of the seafloor as $y = -H(x)$, so that $h(x,t) = \eta (x,t) + H(x)$. Linearization is justified assuming $|\eta (x,t)| \ll |H(x)|$. The equations describe nondispersive wave propagation at speed $\sqrt{gH}$.

Fig. 1

Geometry described by the 1-D shallow water wave Eqs. (1) and (2). Ocean depth is $h = \eta (x,t) + H(x)$, and depth-averaged horizontal velocity is v(x,t)

Time-dependent seafloor displacement, which generates the tsunami, can be accounted for through an additional term in the mass balance Eq. (1). However, it is more common to account for seafloor displacement indirectly through initial conditions. This is justified when the tsunami propagation distance over the rupture duration is much smaller than horizontal scales of interest. This condition is almost always satisfied because earthquake ruptures propagate at elastic wave speeds that are an order of magnitude faster than the tsunami propagation speed (Kajiura 1970).

Solving Eqs. (1) and (2) requires two initial conditions,

$$\eta (x,0) = \eta _0(x),$$

(3)

$$v(x,0) = v_0(x).$$

(4)

These initial conditions are, of course, determined by the tsunami generation process. For a tsunami generated by an undersea earthquake on an initially flat seafloor (i.e., $H(x)={\mathrm {constant}}$), the simplest choice for $\eta _0(x)$ is to assume that the vertical static displacement of the seafloor following the earthquake, $u_y(x)$, is mirrored at the sea surface, so that

$$\eta _0(x) = u_y(x).$$

(5)

A more nuanced approach involves filtering out short wavelength variations to account for the nonhydrostatic ocean response at horizontal wavelengths smaller than H when translating the seafloor uplift to the sea surface (Kajiura 1963; Saito and Furumura 2009; Saito 2013).

Fig. 2

Snapshots from a shallow water wave model for wave height $\eta (x,t)$ (a) and (b) and horizontal velocity v(x,t) (c) and (d). The blue and red simulations both have the same initial condition on $\eta _0(x)$, as in (a), and both have oceans with constant 4 km depth. Blue has $v_0(x) = 0$ and red has $v_0(x) \ne 0$, as seen in (c). The nonzero initial condition on horizontal velocity leads to a tsunami with an asymmetric sea surface profile, whereas $v_0(x) = 0$ produces a symmetric tsunami for an ocean of constant depth

Fig. 3

Geometry of full-physics simulations. A shallow dipping thrust fault runs through the Earth, intersecting the seafloor at the trench. Material properties are piecewise homogeneous in the Earth and ocean, and $\theta$ defines the seafloor slope near the trench

If the seafloor is not flat, horizontal displacements also uplift the sea surface, as recognized by Tanioka and Satake (1996). For horizontal displacement $u_x$ along a seafloor with slope $\mathrm{{d}}H/\mathrm{{d}}x$, the contribution to $\eta _0(x)$ from horizontal displacements is $u_x \,\mathrm{{d}}H/\mathrm{{d}}x$, and thus the appropriate initial condition on sea surface displacement is

$$\eta _0(x) = u_y(x) + u_x \frac{\mathrm{{d}}H}{\mathrm{{d}}x}.$$

(6)

Because subduction zone plate boundary faults tend to dip at shallow angles of $\sim$10$^{\circ }$, $u_x$ is often significantly larger than $u_y$, making the contribution from horizontal seafloor displacement significant even for gently sloping continental shelves. The role of horizontal seafloor displacement, as described by Tanioka and Satake, has gained widespread acceptance in the tsunami modeling community, and Eq. (6) has been widely used to match tsunami data and to make probabilistic estimates of tsunami hazard (e.g., Fujii et al. 2011; Geist and Parsons 2006; Sato et al. 2011).

With regard to initial velocity, the standard approach assumes that the initial depth-averaged horizontal velocity is negligible, i.e.,

$$v_0(x) = 0.$$

(7)

This assumption is exactly correct for an inviscid ocean with a flat seafloor (e.g., Saito 2013) since in that case shear traction vanishes at the seafloor and the normal traction there is exactly vertical. But it is not immediately obvious that the initial horizontal velocity should be zero when the seafloor is sloping and is displaced horizontally. In this case, normal traction on the seafloor has a horizontal component that will impart horizontal momentum from the solid Earth to the ocean. This momentum exchange will certainly occur, as pointed out by Song et al. (2008) based on the impulse-momentum theory, but the relevant question is whether this momentum plays a non-negligible role in tsunami generation. In a compressible ocean, it is also possible for this momentum to be carried away from the source region by ocean acoustic waves.

Fig. 4

Initial condition for depth-averaged horizontal velocity used for shallow water wave code, based on theory of Song et al. (2008). For the region over the sloping seafloor, final seafloor displacement from the full-physics model is divided by a rise time and depth integrated over the height of the ocean column

Fig. 5

Comparisons of shallow water wave simulations with full-physics simulations for $\theta = 10^{\circ }$, where each row uses a different value of rise time $\tau$ when setting initial velocity in the shallow water model. Left column Initial condition on depth-averaged velocity used for the shallow water code, as per Eqs. (8) and (9). We run the shallow water code with two different values of $L_H$ (see legend). Right column Tsunami profiles at the final simulation time, comparing the full-physics result to the two different shallow water simulations

We briefly note an important consequence of choosing a nonzero value for $v_0(x)$, which might be testable with tsunami waveform observations if there is sufficient uniqueness in the solution of the inverse problem of determining initial conditions from data. A nonzero initial velocity introduces asymmetries in the sea surface profile as the tsunami propagates, as illustrated by the simple example shown in Fig. 2. A tsunami model with an initial condition other that (7) must justify this asymmetry. To the best of our knowledge, this asymmetry is not seen in tsunami observations, though it might be challenging to detect if the initial $v_0(x)$ is rather small and bathymetric effects complicate the tsunami waveform. We return to this point subsequently.

In recent years, Song et al. (2008, 2017), Song and Han (2011), Titov et al. (2016), and Xu and Song (2013) have promoted an alternative theory of tsunami generation that includes, in addition to the processes encapsulated by Eq. (6), an additional horizontal momentum transfer from the solid Earth to tsunamis in the ocean that leads to non-negligible initial velocity $v_0(x)$. They suggest that this momentum transfer increases the kinetic energy of the ocean and substantially contributes to tsunami wave heights. For example, they claim this effect accounts for two-thirds of the height of the 2004 Indian Ocean tsunami (Song et al. 2008) and half of the 2011 Tohoku, Japan, tsunami (Xu and Song 2013). If the Song et al. theory is correct and its effect on tsunami height is of the magnitude they claim, it would highlight a major gap in our previous understanding of tsunami generation. The theory therefore merits a thorough investigation.

Fig. 6

Geometry of the Tohoku simulation in this study. We nucleate a dynamic rupture on a thrust fault that intersects the seafloor at the trench. Material properties in the Earth are the same as in Kozdon and Dunham (2013, 2014)

Fig. 7

Comparisons of shallow water simulations with full-physics simulations for the Tohoku geometry, where each row uses a different value of rise time $\tau$ when setting initial velocity in the shallow water model. Left column Initial condition on depth-averaged velocity used for the shallow water code, as per Eqs. (8) and (9). Right column Tsunami profile at 500 s

Drawing from impulse-momentum theory, Song et al. propose that the force exerted by the moving seafloor slope on the ocean transfers horizontal momentum to the ocean. They capture this effect in an approximate manner by estimating the accelerated horizontal particle velocity $\Delta v(x,y)$ of water in an incompressible ocean in the vicinity of a moving slope as

$$\Delta v(x,y) = {\left\{ \begin{array}{ll} u_x(x)/\tau , &{} -H(x) \le y \le -R_x \\ 0, &{} {\mathrm {otherwise}} \end{array}\right. } , \;\;\; R_x = H(x) - \min \bigg \{ H(x),L_H \frac{\mathrm{{d}}H}{\mathrm{{d}}x} \bigg \}$$

(8)

where $\tau \sim$ 5–20 s is the rise time, and $L_H$ is the effective range of horizontal motion (Xu and Song 2013; Titov et al. 2016; Song et al. 2017). In more than one publication, Song and Han (2011) and Xu and Song (2013) have used $L_H = 1.5H_\mathrm{{max}}$, where $H_\mathrm{{max}}$ is the maximum depth of the ocean. However, Song et al. (2017) have recently suggested that a position-dependent effective range, $L_H = 0.5 H(x)$, might be more appropriate. The choice of $L_H$ has a large impact on the amplitude of tsunami waves, and so we examine both versions of $L_H$ in this study. The velocity given by (8) is depth integrated, and serves as an initial condition on velocity for a shallow water wave model, such that

$$v_0(x) = \frac{1}{H(x)} \int _{-H(x)}^0 \Delta v(x,y)\mathrm{{d}}y.$$

(9)

In this paper, we rigorously test the tsunami generation theory of Song et al. by running several “full-physics” earthquake and tsunami models (details provided subsequently) and comparing them to the results of a simple shallow water model. By varying the initial horizontal velocity condition of the shallow water model, we can determine whether or not the initial condition described by Song et al. provides an accurate match to tsunami waveforms in the more realistic full-physics model.

Fig. 8

Misfit between the sea surface height of full-physics simulations and shallow water simulations, calculated at the final time step plotted in previous figures using the $L_2$ norm. a Misfits for the idealized structure in Fig. 3, calculated at $t=600$ s as in Fig. 5. b Misfits for the Tohoku geometry in Fig. 6, calculated at $t=500$ s as in Fig. 7. For both sets of simulations and both choices of $L_H$, misfit is minimized as rise time $\tau$ approaches infinity (initial velocity $v_0 \rightarrow 0$)

Fig. 9

Kinetic energy as a portion of the total energy in the ocean calculated at time $t=0$, for shallow water models with a the idealized structure and b the Tohoku structure. As $\tau$ increases, initial velocity and kinetic energy decrease while gravitational potential energy is unchanged. By comparing with Fig. 8, we observe that models with low misfit are those with negligible initial KE

Modeling approach

As mentioned above, we perform two parallel sets of tsunami simulations, one with a depth-resolved, coupled solid Earth-ocean model and the other with a depth-averaged shallow water wave model.

The fully-coupled, or full-physics, model is a 2-D depth-resolved model introduced by Lotto and Dunham (2015), building on earlier work by Maeda and Furumura (2013). We solve the elastic wave equation in the solid Earth and the acoustic wave equation in the compressible, inviscid ocean, and incorporate gravity via a modified free surface boundary condition. Small wave amplitudes, relative to ocean depth, are assumed to justify a linearization of the free surface condition. However, no assumption is made regarding horizontal wavelengths relative to ocean depth; the ocean response becomes nonhydrostatic at short wavelengths, giving rise to well-known effects such as the filtering of short wavelength features in the initial sea surface height and dispersion during tsunami propagation (Kajiura 1963, 1970). In the full-physics model, the earthquake source and all waves are simultaneously generated, unlike the typical approach of first modeling the earthquake and then using the final displacement from the earthquake model to set initial conditions for the tsunami. We use this model to simulate the full seismic, ocean acoustic, and tsunami wavefields. This model contains the full physics of tsunami generation; the seafloor moves dynamically in response to fault rupture and that motion translates directly to motion and wave propagation in the water. Forces exerted by the solid Earth on the ocean generate both ocean acoustic waves and tsunamis.

We employ two structural models in this study. The first is an idealized subduction zone geometry with a homogeneous solid Earth, and the second is based on our previous studies of the 2011 Tohoku earthquake (Kozdon and Dunham 2013, 2014). In the idealized model, an ocean with a sloping seafloor rests above a homogeneous solid Earth, with a planar thrust fault intersecting the seafloor at its deepest point (Fig. 3). The fault dips at $10^{\circ }$. This shallowly dipping fault is typical of subduction zone plate boundaries, especially near the trench (e.g., Kopp and Kukowski 2003; Gutscher et al. 2001). The seafloor bathymetry to the right of the trench is flat, for simplicity. To the left of the trench, the bathymetry is set by a hyperbolic curve that asymptotes to horizontal away from the trench and to a slope of angle $\theta$ near the trench. We vary $\theta$ across simulations, which changes the depth of the trench and the slope of the seafloor. Results with different $\theta$ are all qualitatively similar, so only simulations with $\theta = 10^{\circ }$ are reported here. The vertical position of the seafloor is defined as $y = -H(x)$; $y = -H_1 = -1$ km at the landward (left) side of the domain and $y = -10$ km at the trench for $\theta = 10^{\circ }$. The trench is slightly deeper than in most subduction zones, but this hyperbolic seafloor shape is a simple way to parameterize the geometry while avoiding too many free parameters.

Although real subduction zones often have compliant sedimentary prisms that affect the rupture process (Lotto et al. 2017), here we assume for simplicity in our idealized structural models that the Earth has homogeneous material properties, with s-wave speed $c_s = 4.0\,{\mathrm {km/s}}$, p-wave speed $c_p = 7.0\,{\mathrm {km/s}}$, and density $\rho = 2800\,{\mathrm {kg/m}}^3$. In the ocean, $c_s = 0$, $c_p = 1.5\,{\mathrm {km/s}}$, and $\rho = 1000\,{\mathrm {kg/m}}^3$. The gravitational acceleration is $g=9.8\,{\mathrm {m/s}}^2$.

Fig. 10

Initial conditions for tsunami problems as determined by the adjoint method. The dotted orange line is the final tsunami profile of a zero gravity simulation (the tsunami doesn’t propagate). The green lines are the sea surface height and depth-averaged horizontal velocity fields determined by the adjoint method. These are the most consistent choices for initial conditions $\eta _0(x)$ and $v_0(x)$

Fig. 11

Total energy in the ocean for the full-physics problem described in “Adjoint problem” and “Energy considerations” sections. An earthquake nucleates at time $t=0$, generating acoustic waves and a tsunami. At time $t=1000$ s, we switch to the adjoint (time-reversed) problem and the tsunami propagates backward. The black dashed line marks the time of the “initial conditions” in Fig. 10. At this time, virtually all of the energy in the ocean is gravitational potential energy and none is kinetic energy

The full-physics model initiates with an earthquake on the fault. For the purposes of this study, either a kinematic or dynamic source model could be used. We choose to use a dynamic rupture model where slip is controlled by the interplay between elastodynamics and a fault friction law. For simulations using the idealized structural model, we nucleate rupture on the fault at the point $x = -75\,{\mathrm {km}}$, where $x=0$ defines the position of the trench. The rupture propagates along the entire $\sim$140-km length of the fault according to the rate-and-state friction law (Kozdon and Dunham 2013):

$$\frac{\mathrm{{d}} T}{\mathrm{{d}}t} = \frac{a \bar{\sigma }}{V} \tanh \Big (\frac{T}{a \bar{\sigma }}\Big ) \frac{\mathrm{{d}}V}{\mathrm{{d}}t} - \frac{V}{l} \big [ |T| - \bar{\sigma } f_{ss}(V) \big ],$$

(10)

with the steady-state friction coefficient given by

$$f_{ss}(V) = f_0 - (b-a)\ln (V/V_0),$$

(11)

where T is the shear strength of the fault, $\bar{\sigma }$ is the effective normal stress, V is the slip velocity, $l=0.8$ m is the state evolution distance, $f_0 = 0.6$ is the friction coefficient for steady sliding at reference velocity $V_0 = 1\,\upmu {\mathrm {m/s}}$, $b=0.02$, and $a=0.016$. The dimensionless parameter $b-a$ determines whether friction increases or decreases with changes to V. Since $b-a = 0.004$ for our simulations, friction is velocity-weakening and unstable slip can occur.

In order to account for poroelastic changes in pore pressure, $\Delta p$, as a response to changes in total normal stress, $\Delta \sigma$, we use the linear relation $\Delta p = B \Delta \sigma$, with $B = 0.6$. The poroelastic effect alters effective normal stress as Kozdon and Dunham (2013)

$$\bar{\sigma } = \bar{\sigma }_0 +(1-B)\Delta \sigma,$$

(12)

which serves to partially buffer the changes in effective normal stress. Equation (12) is a limiting case of a model (Cocco and Rice 2002) where the fault is surrounded by highly damaged material; in that limit, B is Skempton’s coefficient.

Initial effective normal stress $\bar{\sigma }_0$ is calculated as the difference between normal stress in the Earth and pore pressure. Normal stress on the fault is set to lithostatic, and pore pressure is hydrostatic; effective normal stress increases with depth to a maximum of $\bar{\sigma }_0 = 40$ MPa, at which point pore pressure increases lithostatically. Initial shear stress on the fault, $T_0$, is set to $T_0 = 0.6 \bar{\sigma }_0$.

We compare the results of the full-physics model to those of a 2-D linearized shallow water wave model (i.e., Eqs. 1, 2), for which we use the same bathymetry as the full-physics model and set initial conditions based on the results of the full-physics model as follows. To set the initial sea surface height, $\eta _0(x)$, we first run a simulation using the full-physics model where we set the gravitational acceleration to $g=0$ so that surface gravity waves do not propagate. After seismic and ocean acoustic waves propagate out of the domain, the sea surface is statically perturbed in a profile that resembles the seafloor displacement but with short wavelength variations filtered out (Kajiura 1963; Lotto and Dunham 2015; Saito 2013). This static sea surface perturbation becomes the initial condition on $\eta _0(x)$ for the shallow water model.

We set the initial condition on depth-averaged velocity $v_0(x)$ using the expressions provided by Song et al. (2008) and are summarized in Eqs. (8) and (9) and Fig. 4. Here, $u_x(x)$ is taken from the final displacement of the seafloor in the zero gravity version of the full-physics simulation, $R_x$ comes from the model geometry, and the rise time $\tau$ is a free parameter that we vary to control the amplitude of the initial velocity.

Comparing full-physics to shallow water simulations

Here, we compare the results of the shallow water model to those of the full-physics model, using the idealized material structure and geometry for both models. In order to make a proper comparison, we must first account for a time delay that comes about because the full-physics model begins with an earthquake rupture that generates the tsunami over tens of seconds whereas the shallow water model is initialized at time $t=0$ with a tsunami profile already in place. We calculate the time delay by comparing $\eta (x,t)$ at each time step of the full-physics simulation to the initial profile $\eta _0(x)$ of the equivalent shallow water simulation (which, as described previously, is set by a zero gravity full-physics simulation). We calculate the misfit with respect to $\eta _0(x)$ at each time step using the $L_2$ norm, with the time of minimum misfit is taken to be the time delay between the full-physics model and the shallow water model. Alternatively, we can calculate the time delay by selecting some time step of $\eta (x,t)$ in the full-physics model after which all acoustic waves have exited the domain and matching $\eta (x,t)$ of the shallow water model, subject to some time delay, to that in the full-physics model. The time delay from this alternate approach yields basically the same results as the approach previously described. Though the time delay changes slightly between simulations, the offset is typically $\sim$85 s, meaning that it takes about 85 s (roughly, the rupture duration) for the full-physics model to generate a tsunami profile that matches the initial condition of the shallow water code.

We run two sets of shallow water simulations with initial horizontal velocity set by Eqs. (8) and (9), one using $L_H = 1.5H_\mathrm{{max}}$ and the other using $L_H = 0.5H(x)$. The only remaining free parameter is the rise time $\tau$, for which Song et al. use a value of $\tau \sim 5$–20 s for their Sumatra and Tohoku simulations (Song and Han 2011; Xu and Song 2013; Song et al. 2017). We vary $\tau$ over a much larger range, as we find that using $\tau \sim 5$–20 s leads to tsunami waves of unrealistically large amplitudes, as evident in the bottom rows of Fig. 5.

Rather, we find that the shallow water tsunami best fits the full-physics tsunami when rise time $\tau$ approaches infinity (top rows of Fig. 5), i.e., when $v_0(x) = 0$. Note that even when tsunami amplitudes match, the waveforms do not align perfectly because the tsunami experiences dispersion of short wavelengths in the full-physics model but not in the shallow water model. With $\tau$ interpreted as the local rise time, setting $\tau = \infty$ is inconsistent with the assumption that the seafloor completes its deformation over a timescale smaller than the tsunami propagation time over horizontal wavelengths characterizing the deformation profile. However, more generally $\tau$ can be viewed simply as a scalar parameter that controls the amplitude of the initial velocity, with the shape or profile still given by Eqs. (8) and (9).

Shallow water models with non-negligible initial velocity (i.e., those with finite rise times) produce the same sort of asymmetry in tsunami waveform as illustrated in the simple example in Fig. 2. In contrast, the full-physics simulation yields a tsunami profile that is nearly symmetrical in amplitude between seaward and landward waves, though this is complicated somewhat by the nonuniform ocean depth.

We can quantify the misfit between the tsunamis produced by the shallow water and full-physics models by taking the $_2$ norm of the difference between the final sea surface heights plotted in Fig. 5. Those misfits are shown in Fig. 8a. The misfit is never zero because the full-physics model accounts for the finite source duration, surface gravity wave dispersion, and ocean compressibility, which are neglected in the shallow water model. For either choice of $L_H$, the smallest misfits occur in the limit as rise time approaches infinity, i.e., $\tau \rightarrow \infty$ so $v_0 \rightarrow 0$. When using $L_H = 0.5 H(x)$, we find that misfits are nearly at their minimum for $\tau \ge 40$ s. In Fig. 9a, we plot kinetic energy (KE) as a portion of the total energy in the ocean at the initial time step of each shallow water simulation. We define energy and discuss it more thoroughly in “Energy considerations” section, but for now it is sufficient to view KE as a quantitative scalar measure of particle velocities in the ocean. By comparison with Fig. 8a, we observe that models with minimum misfit are the same as those with negligible initial KE. In other words, even though models with $L_H = 0.5 H(x)$ and $\tau \ge 40$ s produce small misfits, they only do so because they essentially have zero initial horizontal velocities.

Tohoku simulations

To supplement the idealized set of simulations discussed in the previous section, we perform a similar analysis on a 2-D simulation of the 2011 Tohoku-oki earthquake and tsunami. The Tohoku simulation, whose seafloor geometry is pictured in Fig. 6, is nearly identical to the simulation performed by Kozdon and Dunham (2013, 2014) with velocity-strengthening friction near the trench and a heterogeneous material structure. Unlike that prior work, the Tohoku simulation presented here accounts for gravity and tsunami generation and propagation, which allows us to test the Song et al. momentum transfer hypothesis for a more realistic problem.

Note that, as seen in Fig. 6, the ocean layer rises above and covers the land in an unphysical way. This is an artifact of the way we break our full-physics model into curvilinear computational “blocks,” but it does not affect our results as long as we are interested in tsunami propagation prior to wave arrival at the coast.

Just as we did for the simplified geometry in the previous section, we run two shallow water simulations setting initial horizontal velocity with Eqs. (8) and (9), one using $L_H = 1.5H_\mathrm{{max}}$ and the other using $L_H = 0.5H(x)$. The bathymetry of the full-physics simulation determines H(x) in the shallow water models, and the initial condition $\eta _0(x)$ comes from the final displacement of the sea surface in a full-physics Tohoku simulation with no gravity. As with the idealized geometry simulations, the time delay between the full-physics and shallow water simulations is calculated by finding the time step of least misfit between $\eta (x,t)$ of the full-physics model and $\eta _0(t)$ of the shallow water model. For the Tohoku geometry, that misfit is $\sim$70 s.

The results, shown in Fig. 7, demonstrate once again that the best choice for rise time using the Song et al. approach in Eq. (8) is $\tau = \infty$, i.e., $v_0(x) = 0$. We quantify the misfits between the final tsunamis produced by the shallow water and full-physics models in Fig. 8b, and we observe that misfits are lowest as $\tau \rightarrow \infty$ and $v_0 \rightarrow 0$. Just as in the previous section, misfits are at their minimum when initial kinetic energy (KE) in the ocean is negligible (Fig. 9b), i.e., when $v_0(x) = 0$. These results do not indicate that $v_0(x) = 0$ is the correct, or even most self-consistent, choice for an initial condition on horizontal velocity for shallow water models, as this experiment explores only two possible initial velocity profiles. However, the results do imply that with these velocity profiles, the Song et al. theory makes inconsistent predictions of the tsunami, especially when choosing $\tau = 5{-}20$ s. We furthermore see no evidence arguing against using $v_0(x)=0$ as the initial velocity for most tsunami modeling problems.

Adjoint problem

Next, we use an alternative approach to determine tsunami initial conditions that avoids any assumptions regarding the initial velocity profile. In concept, the optimal way to constrain the initial condition on horizontal velocity is through direct examination of the particle velocity field of our full-physics model. If we integrate horizontal velocity over the ocean depth at an appropriate “initial” time, we could simply use that horizontally variable, depth-averaged velocity profile as an initial condition for a shallow water simulation. However, this approach is unfeasible: Setting aside the question of how to select an initial time, our full-physics model produces large amplitude seismic and ocean acoustic waves that interfere with the wavefield of the tsunami for the first few hundred seconds. This makes it impossible to isolate the tsunami particle velocity from the total particle velocity for use in setting the initial condition on $v_0(x)$.

However, appropriate (i.e., self-consistent) tsunami initial conditions can be determined using the following adjoint wavefield procedure. We first run a full-physics simulation (with the same idealized structure as in Fig. 3) for 1000 s of simulation time, long enough for all seismic and acoustic waves to exit through the boundaries of the domain, leaving only the slower-propagating tsunami in the ocean. Having isolated the tsunami, we then use the stress and velocity fields at the end of this first simulation as an initial condition on a time-reversed adjoint simulation. In addition to running backward in time, the adjoint simulation involves flipping the sign of the particle velocity. Since the wave propagation problem is self-adjoint, no changes are required to the code to perform an adjoint simulation. In the adjoint simulation, the tsunami propagates backward in time, and landward and seaward propagating tsunami waves converge above the earthquake source into the initial tsunami profile. We determine the initial condition by finding the time where the sea surface profile from the adjoint simulation best matches the sea surface profile of the equivalent forward simulation with zero gravity (as in the top of Fig. 10). The sea surface height and depth-averaged horizontal velocity fields at this time are identified as $\eta _0(x)$ and $v_0(x)$, the initial conditions that are most consistent with the full-physics simulations.

Not surprisingly, we can find a time step in the adjoint problem where the tsunami profile very closely matches that of the equivalent zero gravity simulation. At this time step, the depth-averaged horizontal velocity in the ocean is of order $10^{-2}$ m/s and appears in a profile that does not match either of those profiles suggested by Song et al. The adjoint calculation suggests that the self-consistent initial condition is one that has nonzero, but very small, horizontal velocity. This result does not corroborate the shape of either of the profiles proposed by Song et al., nor is the velocity amplitude consistent with predicted initial velocity using the rise times that Song et al. recommend.

Energy considerations

Song et al. (2008, 2017), Song and Han (2011), and Titov et al. (2016) have frequently brought the concept of energy into their discussion of tsunamis, making estimates of kinetic energy (KE) and gravitational potential energy (GPE) in the ocean, immediately after the earthquake prior to tsunami propagation, and attributing the kinetic energy to momentum transfer from the horizontal motion of the seafloor during tsunami generation. In this section, we calculate the energy balance of the ocean for the full-physics model using the idealized geometry discussed in a previous section.

In addition to KE and GPE, in our full-physics framework we must additionally consider the stored compressional energy in the ocean, which we refer to as acoustic energy (AE) (Lotto and Dunham 2015). More precisely, the AE is the recoverable internal energy associated with compressing or expanding the water. We define the total energy in the ocean, E, as the sum of these three components:

$$\begin{aligned} E&= {\text {KE}}+{\text {AE}}+{\text {GPE}}\nonumber \\&=\int _{-H(x)}^0 \int _{-D}^{D}\frac{1}{2}\rho \varvec{v}^2\, {\mathrm {d}}x{\mathrm {d}}y + \int _{-H(x)}^0 \int _{-D}^{D} \frac{1}{2K}p^2 \,{\mathrm {d}}x{\mathrm {d}}y + \int _{-D}^{D} \frac{1}{2}\rho g \eta ^2 \,{\mathrm {d}}x, \end{aligned}$$

(13)

where $\varvec{v}$ is the total velocity, p is the pressure, $K = \rho c_p^2$ is the bulk modulus of water, and $D=5000$ km is the horizontal distance from the trench to either end of the domain. Energy can escape out of the side boundaries as waves pass through the edges of the domain.

Figure 11 tracks the components of E over time for the same full-physics simulation discussed in “Adjoint problem” section. An earthquake nucleates at time $t=0$, generating seismic waves and deforming the elastic solid. The solid Earth does work on the ocean through normal tractions on the seafloor, transferring energy into the ocean in the form of acoustic waves and surface gravity waves. The fast-moving acoustic waves propagate out of the domain within the first few hundreds of seconds of simulation time, carrying with them much of the KE and AE. After about 300 s, only the tsunami remains within the domain, with relatively constant KE and GPE of approximately equal amplitudes and negligible AE.

At $t=1000$ s, we begin to run the adjoint problem with reversed time. In the adjoint problem, there is no source of acoustic waves and so the AE continues to be negligible. Both GPE and KE remain relatively constant until the tsunamis that were propagating in opposite directions begin to interfere with each other in the source region. This interference establishes the tsunami initial conditions, with the dashed line in Fig. 11 corresponding to the initial conditions shown in Fig. 10. At this time, virtually all of the energy in the ocean is GPE, with negligible KE, such that KE/(KE + GPE) = 0.002. This observation argues against the claim of Song et al. that the KE of the initial tsunami is higher (or even significantly higher) than the GPE.

However, it is important to note that our energy balance calculation provides quantitative support for the impulse-momentum theory promoted by Song et al. There is indeed a considerable amount of KE imparted to the ocean by the solid Earth during the tsunami generation process. However, almost all of this KE is coupled with AE in the form of acoustic waves, rather than remaining in the vicinity of the tsunami source region to provide a non-negligible contribution to initial velocity.

Discussion and conclusions

We have taken two unique approaches in our attempt to determine the appropriate initial condition on horizontal velocity $v_0(x)$ for shallow water tsunami models. The first approach involves running comparable tsunami simulations with a full-physics model and a linearized shallow water wave model, and selecting the initial condition on horizontal velocity for the shallow water model using the expressions presented by Xu and Song (2013), Titov et al. (2016) and Song et al. (2017). We then vary the rise time parameter $\tau$, which scales the amplitude of $v_0(x)$, and compare the tsunamis predicted by the shallow water and full-physics simulations.

This comparison is made for two model geometries: an idealized subduction zone geometry (Fig. 3) and a 2011 Tohoku earthquake geometry (Fig. 6). For both geometries, our simulations unequivocally demonstrate that the initial condition on $v_0(x)$ (Eqs. 8, 9) using $\tau =$ 5–20 s, as Song et al. have suggested, overestimates tsunami amplitudes by several times. Tsunami amplitudes are nearly an order of magnitude too large when choosing $L_H = 1.5H_\mathrm{{max}}$, as Song et al. have done in more than one publication (Song and Han 2011; Xu and Song 2013). If we use larger values for $\tau$ – which may be more physically relevant because the rupture process for typical tsunamigenic earthquakes lasts several tens of seconds—we find that the shallow water model produces results closer to the full-physics simulations, especially when $L_H = 0.5 H(x)$. But when $\tau \ge 80$ s and $L_H = 0.5 H(x)$, the initial horizontal velocity profile is nearly flat, virtually indistinguishable from choosing $v_0(x) = 0$, i.e., $\tau = \infty$. Quantitative measures of misfit between the two models (Fig. 8) support the conclusion that shallow water simulations best match their full-physics counterparts as rise time $\tau$ approaches infinity. Indeed, comparing Figs. 8 and 9 reveals that any shallow water simulation that reasonably matches its full-physics counterpart does so because initial KE—and, accordingly, initial horizontal velocity—is negligible.

The second approach requires us to run a full-physics simulation long enough for all seismic and acoustic waves to leave the domain, and then to run it backward in time to solve the adjoint problem. This adjoint simulation method allows us to use the full-physics simulation to directly determine $v_0(x)$. In doing so, we again find that $v_0(x)$ is so small in magnitude so as to be practically equivalent to the null hypothesis, $v_0(x) = 0$. This method produces a velocity profile that does not qualitatively match that suggested by Song et al. Most notably, the initial velocity profile determined by the adjoint method does not reach its maximum at the trench ($x=0$) and is not uniformly positive. We also find that at the time of this adjoint-based initial condition, almost all the energy in the ocean is in the form of gravitational potential energy and virtually none is kinetic energy.

We conclude that Song et al.’s proposed choice of $v_0(x)$ for modeling tsunamis will consistently overestimate tsunami heights. Furthermore, we have found no convincing evidence of any better choice than $v_0(x) = 0$, which is the default initial condition on horizontal velocity used by most tsunami modelers.

In a study published after initial submission of our manuscript, Song et al. tested their tsunami generation theory by comparing simulations using initial condition (8 and 9) to tide gauge and DART data from the 2011 Tohoku event (Song et al. 2017). They argued that their tsunami generation theory matches the data better than the conventional tsunami generation theory. However, many other authors have successfully used the standard initial condition $v_0(x) = 0$ to fit data both landward and seaward of the rupture area (Bletery et al. 2014; Gusman et al. 2012; Hooper et al. 2013; Satake et al. 2013). The fits are equally good. This suggests that there is sufficient nonuniqueness with regard to the tsunami initial conditions given the available tsunami data from the Tohoku event. It follows that we cannot equivocally disambiguate between Song’s initial conditions and the conventional initial conditions from even this well-recorded tsunami.

Song et al. have also used laboratory wave basin tests to create an empirical scaling relationship between the ratio of initial kinetic to gravitational potential energy (KE/GPE) and a normalized version of tsunami amplitude (Song et al. 2017). Using characteristic values for ocean depth, fault slip, rise time, and tsunami amplitude, they fit the 2011 Tohoku and 2004 Sumatra tsunamis into this paradigm. A full discussion of these laboratory experiments is beyond the scope of this paper, but we note certain flaws in the scaling relations of the experiments. First, tsunami amplitudes in the open ocean are much smaller than ocean depths ($|\eta | / H \ll 1$), but in these experiments $|\eta |/H \approx 0.1 - 0.7$. Nonlinearity is therefore essential in the experiments, while negligible in the real ocean. Second, water compressibility effects are vastly less important in the experiments than in the real ocean. Two dimensionless parameters quantifying compressibility effects are $gH/c^2$ and $H/c \tau$, where c is the sound speed in the ocean and $\tau$ is again the rise time or source duration. These dimensionless parameters, and therefore compressibility effects, are several orders of magnitude more important in real tsunamis than in the experiments. These differences in scaling suggest that waves in these laboratory experiments are generated and propagate in a different manner from real tsunamis.

What accounts for the stark difference between the results in this paper and in the studies of Song and collaborators? In our full-physics simulations, we observe large amplitude acoustic waves in the ocean, likely due in part to the transfer of horizontal momentum from the solid Earth to the ocean during earthquake rupture. Our full-physics model, which includes a compressible ocean, captures these acoustic waves that a model using an incompressible ocean would miss. In an incompressible ocean with a finite domain, horizontal momentum transferred from the solid Earth to the ocean would be left in the vicinity of the source region rather than being carried away by acoustic waves. However, in a full-physics, compressible ocean model most of the kinetic energy injected into the ocean by the seafloor displacement is carried away by acoustic waves and does not significantly contribute to tsunami height. Overall, we conclude that the standard practice of setting initial horizontal velocity to zero remains the best, and most consistent, procedure for tsunami modeling.