# A preliminary depth-integrated model for tsunamis propagation including water compressibility and seafloor elasticity

Gaël Loïc Richard, K Msheik, Arnaud Duran

► To cite this version:

Gaël Loïc Richard, K Msheik, Arnaud Duran. A preliminary depth-integrated model for tsunamis propagation including water compressibility and seafloor elasticity. Electronic Journal of Probability, 2023, 10.1016/j.euromechflu.2023.01.004. hal-03805611

HAL Id: hal-03805611

https://hal.inrae.fr/hal-03805611

Submitted on 7 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A new depth-integrated model for tsunamis
propagation
including water compressibility and seafloor
elasticity

G. L. Richard†, K. Msheik*, and A. Duran*

†Univ. Grenoble Alpes, INRAE, UR ETNA, 38000 Grenoble, France

*Institut Camille Jordan, CNRS UMR 5208, Université Claude Bernard Lyon 1, 69100 Villeurbanne, France

Abstract

A model for tsunamis propagation is derived by coupling a weakly
compressible liquid layer, standing for the ocean, to a viscoelastic solid
layer of constant depth, representing the upper layer of the Earth. Be-
low this layer, the Earth is assumed to be perfectly rigid. The resulting
equations are averaged over the depth of the water layer for the liq-
uid part and are integrated over the thickness of the solid layer for the
solid part. The obtained system of equations is hyperbolic and admits
an exact equation of energy conservation. The model is dispersive and
includes an elastic branch, an acoustic branch and a gravity branch.
The elasticity of the solid layer entails: 1) a reverse dispersion effect
i.e. the phase velocity of the gravity branch decreases for small wave
numbers if the wave number decreases; 2) an arrival time delay of the
tsunami and 3) a leading negative phase i.e. a negative water eleva-
tion before the arrival of the main wave. All these effects are absent if
the sea floor is assumed to be perfectly rigid. The arrival time delay
due to elasticity comes in addition to the delay due to seawater com-
pressibility. These effects are in agreement with previous works on
tsunamis propagation. This preliminary model is a promising alter-
native approach for tsunamis propagation due to its ability to capture
the main physical effects with a relatively simple mathematical struc-
ture of equations, which is easy to solve numerically.
# 1 Introduction

Models for the propagation of tsunamis commonly use the linear shallow water equations. This approach assumes an incompressible sea water and a perfectly rigid sea floor. Moreover dispersive and nonlinear effects are neglected in these models.

Nonlinear terms can be important in coastal waters, in particular for the propagation of tsunamis into bays and harbours, but their effects are small in the deep ocean (Baba et al. 2017) [6] since the amplitude of a tsunami is very small with respect to the water depth. Dispersive effects can be included in the models by adding the Boussinesq terms in the shallow water equations, which leads to the Boussinesq-type equations. Kirby et al. (2013) [14] showed the importance of dispersive effects both in the far field of large co-seismic tsunamis as well as in the near field of more concentrated seismic sources such as submarine mass failures. In the case of the 2011 Tohoku event, Grilli et al. (2013) [10] found that dispersive effects were negligible in the near-field but that they may account for 20–40 % of tsunami amplitude in deep water, hence justifying the use of a Boussinesq-type model. Baba et al. (2017) [6] found that the Boussinesq term has a strong effect by changing the shape of the tsunami waveforms in both the near and far fields and that this was caused by the inherent dispersive effect where the short-wavelength energy propagated more slowly than the large-wavelength energy. Baba et al. (2019) [7] strongly recommended the use of dispersive equations to avoid underestimating the height of tsunamis caused by submarine landslides.

Measurements of tsunamis arrival times at long distance show a systematic delay with respect to the predictions given by the linear shallow water equations or by the Boussinesq-type models. The order of magnitude of this delay is about 1% of the propagation time, giving a 10–30 min delay in the far field. In the case of the 2004 Sumatra tsunami, Rabinovich et al. (2011) [20] found a time shift of 15 min and 10 min at two bottom pressure gauges for a total travel time of 19 h 46 min and 19 h 39 min respectively, giving a difference of 1.26 % and 0.85 %. From the record of the tsunami produced by the 2010 Central Chile earthquake at a GPS buoy located near Japan, Kato et al. (2011) [13] found a 26 min delay for a travel time close to 24 h amounting to a difference of 1.8 %. Similarly, time delays of up to 15 min were found for the 2011 Japanese Tohoku-Oki tsunami with respect to numerical simulations using the linear shallow water equations (Tang et al. 2012 [23]) and of about 7 min (up to 1.5 %) with a fully nonlinear dispersive Boussinesq model (Grilli et al. 2013 [10]).

The arrival time delay is mainly attributed to the effects of the sea water compressibility, the elasticity of the earth and the gravitational potential
variation. The compressibility of water results in a density stratification and a decrease of the tsunami propagation velocity. The sea floor undergoes an elastic deformation due to tsunami loading, which also decreases the propagation velocity.

Tsai et al. (2013) [24] obtained an estimation of the effect of the sea floor elasticity and of the variable water density. Allgeyer & Cummins (2014) [4] calculated the Earth's response to a mass loading its surface with the method of Longman (1962) [16] using a Green's function which can be expressed as a sum over spherical harmonics. The values of the Green's function are tabulated in Pagiatakis (1990) [18] who used the Preliminary Reference Earth Model (PREM) of Dziewonski & Anderson (1981) [8]. The seabed displacement including the effect of elastic loading can be found by convolving the Green's function with the change in ocean depth. This step accounts for about 70% of the computational time (Allgeyer & Cummins 2014 [4]). This method was coupled to the resolution of the nonlinear shallow water equations. This work was extended by Baba et al. (2017) [6] who included the dispersive terms of the Boussinesq equations and incorporated the effect of gravitational potential change in the expression of the Green's function.

Watada et al. (2014) [28] used a synthetic method using the tsunami normal mode theory of Ward (1980) [27] and Okal (1982) [17] in which the tsunamis are regarded as propagating dispersive surface gravity waves in the top fluid layer of a self-gravitating elastic PREM (with a 4 km deep ocean instead of 3 km in the PREM). This method takes into account the elasticity and gravitational potential changes effects. The sea water density was assumed to be constant. This work was extended by Ho et al. (2017) [11] with the inclusion of ocean density stratification, actual raypath and bathymetry.

The potential theory was used to evaluate the effect of water compressibility by Abdolali & Kirby (2017) [2]. The effects of sea-floor elasticity and gravitational potential changes were added in Abdolali et al. (2019) [1].

Apart from the arrival time delay, a small negative phase is often observed before the arrival of the first tsunami positive peak. This leading negative phase is commonly observed at far-field stations but not at near-field stations indicating that the origin of this precursory phase is not related to the earthquake source but to propagation effects (Watada et al. 2014 [28]). This negative phase is related to the effect of reverse dispersion at long periods i.e. the tsunami velocity is smaller if the period increases above $10^3$ s approximately (Watada et al. 2014 [28]). The reverse dispersion effect can be explained by the Earth elastic loading and gravitational potential changes. According to Allgeyer & Cummins (2014) [4], the density stratification results only in a time delay and does not otherwise distort the waveform but
the surface loading changes the arrival time and also distort the waveform by the emergence of a long-period negative peak prior to the positive peak of the main arrival. The negative phase is the result of the loading effect and is due to the downward flexure of the seafloor in response to the weight of the tsunami (Allgeyer & Cummins 2014 [4]).

Inazu & Saito (2013) [12] modelled the elastic loading and gravitational potential changes with a simple parametrization with one scalar parameter using a method originally proposed by Accad & Pekeris (1978) [3] for the effects of tidal self-attraction and loading. This method is faster but less accurate than the convolution method. In particular, it does not account for the change in waveform and is unable to predict the initial negative phase and the reverse dispersion effect.

The importance to take into account these effects in the models was highlighted by several authors (Yoshimoto et al. 2016 [29], Baba et al. 2017 [6], Ho et al. 2017 [11]) who noted that the ability to accurately model the far-field tsunami waveforms can improve real-time tsunami forecasting and could facilitate the use of far-field tsunami data to improve our understanding of tsunamis sources, in particular if near-field data are lacking.

The goal of this work is to derive a model for tsunami propagation including the dispersive and nonlinear effects as well as the effect of water compressibility and elastic loading. This model should be kept as simple as possible while capturing these physical effects. This implies that it should be more complex that the model of Inazu & Saito (2013) [12] to be able to predict the change in waveform, the leading negative phase and the reverse dispersion effect but simpler than the convolution or the normal mode methods. The model could be an extension of the Boussinesq-type models, because of their dispersive properties, but their mixed elliptic-hyperbolic structure with high-order space-time derivatives is not easy to implement numerically and is computationally expensive. Consequently the goal is to extend the hyperbolic model of Richard (2021) [21], which takes into account the compressibility, in order to include the seafloor elasticity. The gravitational changes are not included in the present work and are kept for a future extension. To our knowledge, this is the first attempt to extend a hyperbolic compressible Boussinesq-type model to Earth elasticity. Consequently this work is an exploratory work to determine how Earth elastic effects could be added to such a model and whether the main effects, arrival delay, initial negative phase and reverse dispersion, can be correctly taken into account.

The derivation of the model is explained in §2. The analysis of the model, including its hyperbolicity and its dispersive properties, is detailed in §3. Numerical simulations are presented in §4.
<figure><img src="image_2.png" /><figcaption>Figure 1: Definition sketch.</figcaption></figure>

# 2 Derivation of the model

## 2.1 Physical description and modelling

The propagation of a tsunami is studied in a simplified problem. The flow is two-dimensional and the still water depth, denoted by $h_0$, is constant. The wave is characterized by an elevation $\eta$ of the water with respect to its still level. We assume that the solid below the ocean is not rigid and that the position of the solid/liquid interface can change from a position $z = 0$ at rest to a position $z = b$ if there is a wave, $z$ being the upward vertical coordinate (see figure 1). The free surface is at $z = h_0 + \eta$ and the total depth is $h = h_0 + \eta - b$.

We can describe the phenomenon of sea floor loading as follows: At the arrival of a wave, the extra loading applies an extra pressure on the solid bottom which is a little bit depressed as a result. The sea floor reverts to its initial position after the departure of the wave. This suggests an elastic behaviour reminiscent of a spring. Given that the displacements in the solid part are very small, the solid can be modelled with the constitutive law of linear viscoelasticity since the Earth's crust and mantle exhibit viscoelastic behaviour (see e.g. Yuen & Peltier 1982 [30], Pagiatakis 1990 [18], Wang *et al.* 2012 [26]).

The normal mode theory shows that the tsunami modes are confined to the top tens of kilometers of crust and mantle and that Earth structures deeper than this are of no consequence (Ward 1980 [27]). The characteristic penetration depth can reach 100 or 200 km for tsunami waves of low frequencies such as co-seismic tsunamis but can decrease to 10 or 20 km for tsunamis of higher frequencies as those due to submarine mass failures. We assume thereafter that only the crust and the upper mantle have an effect on tsunami propagation. The normal mode theory shows a decay of eigenfunction amplitude with increasing depth (Ward 1980 [27]). Consequently, at a depth large
enough, of the order of magnitude of the characteristic penetration depth, the displacement is zero.

We model the Earth as a viscoelastic layer with a thickness $H$. At a depth larger than $H$, the Earth is assumed to be perfectly rigid. The value of $H$ can be chosen at one of the discontinuities in the upper layers of the Earth such as the Mohorovičić discontinuity at the boundary between the crust and the mantle or the discontinuity around a 410 km depth at the limit between the upper mantle and the transition zone. These discontinuities appear in all models of the Earth. The existence of a discontinuity around 220 km, the Lehmann discontinuity, is more disputed. It seems weaker below the oceans. All these discontinuities correspond to sharp increases of the density and of the velocity of seismic waves (S and P) and thus of the Lamé coefficients and of the rigidity of the Earth. Even in the absence of a true discontinuity, a model of solid layers with uniform density and elastic properties separated by discontinuities can be seen as a simplified description of the real Earth structure. The simplified model used to explore the possibility to include Earth elasticity in a system of equations for tsunami propagation is thus a two-layer model with an upper layer of weakly compressible water of thickness $h$ and a lower layer of viscoelastic solid of thickness $H$, with constant density and Lamé coefficients, above a rigid Earth. A more accurate model could be a model with several solid layers but this is left for a future work. For this work the goal is only to determine whether a relatively simple depth-integrated model can predict the main physical effects of solid Earth elasticity on tsunami propagation, particularly the leading negative phase and the reverse dispersion effect.

Several different viscoelastic models have been used for the Earth's mantle such as the models of Maxwell or Burgers. The model of Maxwell is more appropriate for long-term processes such as glacial isostatic adjustment. In the present work, we need a simple viscoelastic model for a solid behaviour able to reversibly revert to its initial state. The model of Kelvin-Voigt is the simplest one for this goal. It can be described as an elastic spring and a viscous damper in parallel. The strain is the same in each component and the total stress $\sigma_T$ is the sum of an elastic part $\sigma$ governed by the generalized Hooke's law and of a viscous part $\sigma'$:

$$ \sigma_T = \sigma + \sigma' \qquad (1) $$

The viscous term can be written

$$ \sigma' = 2\mu'D + \eta'(\text{tr}D)I \qquad (2) $$

where $D$ is the strain rate tensor, $I$ the identity tensor, $\mu'$ the first viscosity and $\eta'$ the second viscosity. For the elastic part, we adopt the hyperbolic
formulation of Virieux (1986) [25] which is best suited to derive a hyperbolic model (see also LeVeque 2002 [15]). Denoting by **u** the displacement in the solid and by **ε** the linear strain tensor

$$ \boldsymbol{\varepsilon} = \frac{1}{2} \left[ \mathbf{grad} \, \mathbf{u} + (\mathbf{grad} \, \mathbf{u})^{\mathrm{T}} \right] \qquad (3) $$

Hooke's law can be written

$$ \boldsymbol{\sigma} = 2\mu\boldsymbol{\varepsilon} + \lambda(\mathrm{tr}\boldsymbol{\varepsilon})\boldsymbol{I} \qquad (4) $$

where $\mu$ and $\lambda$ are the Lamé coefficients. For small perturbations, the elastic part of the stress satisfies the equation (Virieux 1986 [25])

$$ \frac{\partial \boldsymbol{\sigma}}{\partial t} = 2\mu \boldsymbol{D} + \lambda(\mathrm{tr}\boldsymbol{D})\boldsymbol{I} \qquad (5) $$

In the case of linear viscoelasticity, Newton's second law can be written

$$ \rho_s \frac{\partial \dot{\mathbf{u}}}{\partial t} = \mathbf{div} \, \boldsymbol{\sigma}_T \qquad (6) $$

where $\rho_s$ and $\dot{\mathbf{u}}$ are respectively the density and the velocity field in the solid layer. Apart from the viscous part (2), which of course induces a parabolic behaviour, the system (5)–(6) is a hyperbolic system of equations.

The liquid layer is described as an inviscid, weakly compressible fluid, as in Richard (2021) [21]. The velocity field, the pressure field and the density in the liquid are denoted by $\mathbf{v}$, $p$ and $\rho_L$ respectively. The continuity equation writes

$$ \frac{\partial \rho_L}{\partial t} + \mathrm{div}(\rho_L \mathbf{v}) = 0 \qquad (7) $$

and the momentum balance equation is

$$ \frac{\partial \rho_L \mathbf{v}}{\partial t} + \mathbf{div} (\rho_L \mathbf{v} \otimes \mathbf{v}) = \rho_L g - \mathbf{grad} p \qquad (8) $$

where $g$ is the gravitational acceleration vector.

For a 2D-flow (see figure 1), the components of the various vectors are denoted as follows: $(u_x, u_z)^{\mathrm{T}}$ for the displacement in the solid, $(\dot{u}_x, \dot{u}_z)^{\mathrm{T}}$ for the velocity in the solid, $(v_x, v_z)^{\mathrm{T}}$ for the velocity in the liquid. The components of the tensors $\boldsymbol{\sigma}$ are $\sigma_{xx}$, $\sigma_{xz}$, $\sigma_{zz}$ and the components of the tensors $\boldsymbol{\sigma}'$ are $\sigma'_{xx}$, $\sigma'_{xz}$, $\sigma'_{zz}$. In the Oz-direction, Equation (6) becomes

$$ \rho_s \frac{\partial \dot{u}_z}{\partial t} - \frac{\partial \sigma_{xz}}{\partial x} - \frac{\partial \sigma'_{xz}}{\partial x} - \frac{\partial \sigma_{zz}}{\partial z} - \frac{\partial \sigma'_{zz}}{\partial z} = 0 \qquad (9) $$
while Equation (5) yields for the components $\sigma_{xz}$ and $\sigma_{zz}$

$$ \frac{\partial \sigma_{xz}}{\partial t} - \mu \frac{\partial \dot{u}_x}{\partial z} - \mu \frac{\partial \dot{u}_z}{\partial x} = 0 \qquad (10) $$

$$ \frac{\partial \sigma_{zz}}{\partial t} - \lambda \frac{\partial \dot{u}_x}{\partial x} - (2\mu + \lambda) \frac{\partial \dot{u}_z}{\partial z} = 0 \qquad (11) $$

The components $\sigma'_{xz}$ and $\sigma'_{zz}$ are given by

$$ \sigma'_{xz} = \mu' \left( \frac{\partial \dot{u}_x}{\partial z} + \frac{\partial \dot{u}_z}{\partial x} \right) \qquad (12) $$

$$ \sigma'_{zz} = 2\mu' \frac{\partial \dot{u}_z}{\partial z} + \eta' \left( \frac{\partial \dot{u}_x}{\partial x} + \frac{\partial \dot{u}_z}{\partial z} \right) \qquad (13) $$

At the liquid/solid interface, the kinematic boundary condition

$$ \dot{u}_z(b) = \frac{\partial b}{\partial t} + \dot{u}_x \frac{\partial b}{\partial x} \qquad (14) $$

and the dynamic boundary condition

$$ -[\sigma_{xx}(b) + \sigma'_{xx}(b)] \frac{\partial b}{\partial x} + \sigma_{xz}(b) + \sigma'_{xz}(b) = p(b) \frac{\partial b}{\partial x} \qquad (15) $$

$$ -[\sigma_{xz}(b) + \sigma'_{xz}(b)] \frac{\partial b}{\partial x} + \sigma_{zz}(b) + \sigma'_{zz}(b) = -p(b) \qquad (16) $$

are satisfied. On the rigid bottom, at $z = -H$, the displacement and the velocity are equal to zero. In particular,

$$ u_z(-H) = 0; \quad \dot{u}_z(-H) = 0 \qquad (17) $$

The boundary conditions for the liquid layer are also the kinematic and dynamic boundary conditions at the free surface, with a constant atmospheric pressure which could be taken equal to zero, and the no-penetration condition at the liquid/solid interface. The derivation of the equations in the liquid layer is treated as in Richard (2021) [21] and is not repeated thereafter. On the other hand the derivation of the equations for the solid layer is detailed below.

## 2.2 Orders of magnitude and scaling

The four characteristic lengths of the problem are the wavelength $L$, the still water depth $h_0$, the wave elevation $\eta$ and the liquid/solid interface elevation
b. The thickness H of the solid layer if of the same order of magnitude as the wavelength L, from a few tens of km to a few hundred km. The still water depth is equal to a few km, typically 4 to 6 km. The wave elevation is of the order of 1 m and b is of the order of a few cm. The ratio

$$ \varepsilon = \frac{h_0}{L} \qquad (18) $$

is the traditional shallow-water parameter. In the case of a tsunami, this is a small parameter: $\varepsilon \ll 1$. Let's take the case $h_0 \sim 4$ km and $L \sim 200$ km. This gives $\varepsilon \sim 2 \cdot 10^{-2}$. With $\eta \sim 1$ m, the nonlinearity parameter $\delta = \eta/h_0 \sim 2.5 \cdot 10^{-4}$ which means that $\delta \sim \varepsilon^2$. The problem is weakly nonlinear. Taking $b \sim 2$ cm gives $b/\eta \sim \varepsilon$. With respect to the length L, we have $h_0/L \sim \varepsilon$, $\eta/L \sim \varepsilon^3$ and $b/L \sim \varepsilon^4$.

The variation of the vertical displacement $u_z$ in the solid is of the order of b. The horizontal displacement on the other hand is much smaller. With the normal mode theory Ward (1980) [27] found that in contrast to the behaviour in the ocean, vertical motion dominates in the solid. We assume that $u_x/u_z \sim \varepsilon$ which gives $u_x/L \sim \varepsilon^5$.

The order of magnitude of $\sigma_{zz}$ is the pressure of the water column of depth $h$, which is close to the hydrostatic pressure $\rho_\ell gh$ where $\rho_\ell$ is a characteristic density of the ocean which is specified below. In this pressure, the main part $\rho_\ell gh_0$ is a constant and only the small part $\rho_\ell g\eta \sim \varepsilon^2 \rho_\ell gh_0$ is variable. The order of magnitude of the other components of $\sigma$ is $\varepsilon^2 \rho_\ell gh_0$.

The values of the Lamé coefficients can be found from the values of the velocities of the seismic waves P and S. Using the values of the PREM, $\mu$ is in the range $2 \cdot 10^{10}-7 \cdot 10^{10}$ Pa and $\lambda$ in the range $3 \cdot 10^{10}-9 \cdot 10^{10}$ Pa. For a depth $h_0 \sim 4$ km, the ratios $\rho_\ell gh_0/\mu$ and $\rho_\ell gh_0/\lambda$ are of the order of $\varepsilon^2$.

Consequently the various quantities are scaled as follows:

$$ \begin{array}{lclcl} b = \varepsilon^4 L \tilde{b} & u_z = \varepsilon^4 L \tilde{u}_z & u_x = \varepsilon^5 L \tilde{u}_x & \eta = \varepsilon^3 L \tilde{\eta} \\ h = h_0 \tilde{h} & z = L \tilde{z} & x = L \tilde{x} & H = L \tilde{H} \\ p = \rho_\ell g h_0 \tilde{p} & \sigma_{xx} = \varepsilon^2 \rho_\ell g h_0 \tilde{\sigma}_{xx} & \sigma_{xz} = \varepsilon^2 \rho_\ell g h_0 \tilde{\sigma}_{xz} & \sigma_{zz} = \varepsilon^2 \rho_\ell g h_0 \tilde{\sigma}_{zz} \\ \mu = \frac{\rho_\ell g h_0}{\varepsilon^2} \tilde{\mu} & \lambda = \frac{\rho_\ell g h_0}{\varepsilon^2} \tilde{\lambda} \end{array} $$

The tildes denote dimensionless quantities.

Due to the constant part of the pressure $\rho_\ell gh_0$ (see above), $\sigma_{zz}$ should be scaled as $\sigma_{zz} = \rho_\ell gh_0 \tilde{\sigma}_{zz}$ in terms that are not a derivative of $\sigma_{zz}$ with respect to $x$ or $t$.

The characteristic time is $L/c_s$ where $c_s$ is the celerity of the transverse waves (similar to S-waves) $c_s = \sqrt{\mu/\rho_s}$. Since $\rho_\ell/\rho_s = O(1)$, we can take
simply $\tau = \sqrt{\varepsilon L/g}$ as a characteristic time and scale $t$ as $t = \tilde{t}\sqrt{\varepsilon L/g}$. The velocities are then scaled as

$$ \dot{u}_z = \varepsilon^4 \frac{L}{\tau} \tilde{u}_z ; \quad \dot{u}_x = \varepsilon^5 \frac{L}{\tau} \tilde{u}_x. \qquad (19) $$

The Reynolds number, defined as $Re = \rho_s L b/(\tau \mu')$, is taken of $O(\varepsilon^4)$ and is written $\varepsilon^4/\tilde{\mu}'$. This leads to the scaling of the viscosities

$$ \mu' = \frac{\rho_s \sqrt{gL^3}}{\sqrt{\varepsilon}} \tilde{\mu}'; \quad \eta' = \frac{\rho_s \sqrt{gL^3}}{\sqrt{\varepsilon}} \tilde{\eta}'. \qquad (20) $$

The components of the viscous part of the stress tensor are then scaled as

$$ \sigma'_{xz} = \varepsilon^2 \rho_s g h_0 \tilde{\sigma}'_{xz} ; \quad \sigma'_{zz} = \varepsilon^2 \rho_s g h_0 \tilde{\sigma}'_{zz}. \qquad (21) $$

With this scaling, Equation (9) becomes in dimensionless form

$$ \frac{\rho_s}{\rho_\ell} \frac{\partial \tilde{u}_z}{\partial \tilde{t}} - \frac{\partial \tilde{\sigma}_{xz}}{\partial \tilde{x}} - \frac{\rho_s}{\rho_\ell} \frac{\partial \tilde{\sigma}'_{xz}}{\partial \tilde{x}} - \frac{1}{\varepsilon^2} \frac{\partial \tilde{\sigma}_{zz}}{\partial \tilde{z}} - \frac{\rho_s}{\rho_\ell} \frac{\partial \tilde{\sigma}'_{zz}}{\partial \tilde{z}} = 0 \qquad (22) $$

Equation (10) can be written

$$ \frac{\partial \tilde{\sigma}_{xz}}{\partial \tilde{t}} - \tilde{\mu} \frac{\partial \tilde{u}_z}{\partial \tilde{x}} = O(\varepsilon) \qquad (23) $$

and Equation (11) becomes

$$ \frac{\partial \tilde{\sigma}_{zz}}{\partial \tilde{t}} - (2\tilde{\mu} + \tilde{\lambda}) \frac{\partial \tilde{u}_z}{\partial \tilde{z}} = O(\varepsilon) \qquad (24) $$

The components of the viscous part of the stress tensor (12)–(13) become, in dimensionless form,

$$ \tilde{\sigma}'_{xz} = \tilde{\mu}' \frac{\partial \tilde{u}_z}{\partial \tilde{x}} + O(\varepsilon) \qquad (25) $$

$$ \tilde{\sigma}'_{zz} = (2\tilde{\mu}' + \tilde{\eta}') \frac{\partial \tilde{u}_z}{\partial \tilde{z}} + O(\varepsilon) \qquad (26) $$

The dimensionless kinematic boundary condition (14) at the interface is simply

$$ \tilde{u}_z(b) = \frac{\partial \tilde{b}}{\partial \tilde{t}} + O(\varepsilon^5) \qquad (27) $$

while the dynamic boundary condition (15) becomes

$$ \tilde{\sigma}_{xz}(b) + \frac{\rho_s}{\rho_\ell} \tilde{\sigma}'_{xz}(b) = O(\varepsilon^2) \qquad (28) $$
Then the condition (16) can be written

$$ \tilde{\sigma}_{zz}(b) + \varepsilon^2 \frac{\rho_s}{\rho_\ell} \tilde{\sigma}'_{zz}(b) = -\tilde{p}(b) + O(\varepsilon^6) \quad (29) $$

In this last equation, the variable part of both $p(b)$ and $\sigma_{zz}(b)$ is of $O(\varepsilon^2)$ with respect to the constant part.

## 2.3 Depth-integrated equations

Equations (22) and (23) are integrated over the depth of the solid layer from $z = -H$ to $z = b$. We define the depth-integrated quantities

$$ q_2 = \int_{-H}^{b} \dot{u}_z \, dz \quad (30) $$

and

$$ S_{12} = \int_{-H}^{b} \sigma_{xz} \, dz \quad (31) $$

The integration of (23) leads to

$$ \frac{\partial}{\partial \tilde{t}} \int_{-\tilde{H}}^{\varepsilon^4 \tilde{b}} \tilde{\sigma}_{xz} \, d\tilde{z} - \tilde{\mu} \frac{\partial}{\partial \tilde{x}} \int_{-\tilde{H}}^{\varepsilon^4 \tilde{b}} \tilde{u}_z \, d\tilde{z} = O(\varepsilon) \quad (32) $$

which yields the depth-integrated equation

$$ \frac{\partial \tilde{S}_{12}}{\partial \tilde{t}} - \tilde{\mu} \frac{\partial \tilde{q}_2}{\partial \tilde{x}} = O(\varepsilon) \quad (33) $$

Taking into account the boundary conditions, the integration of (22) leads to

$$ \frac{\rho_s}{\rho_\ell} \frac{\partial \tilde{q}_2}{\partial \tilde{t}} - \frac{\partial \tilde{S}_{12}}{\partial \tilde{x}} = -\frac{1}{\varepsilon^2} \tilde{p}(\varepsilon^4 \tilde{b}) - \frac{1}{\varepsilon^2} \tilde{\sigma}_{zz}(-\tilde{H}) + \frac{\rho_s}{\rho_\ell} \frac{\partial}{\partial \tilde{x}} \int_{-\tilde{H}}^{\varepsilon^4 \tilde{b}} \tilde{\sigma}_{xz} \, d\tilde{z} - \frac{\rho_s}{\rho_\ell} \tilde{\sigma}'_{zz}(-\tilde{H}) + O(\varepsilon^4) \quad (34) $$

The pressure at the bottom of the ocean layer is $p(b) = -\rho_\ell gh$ or, in dimensionless form, $\tilde{p}(\varepsilon^4 \tilde{b}) = -\tilde{h}$. Equation (25) implies

$$ \int_{-\tilde{H}}^{\varepsilon^4 \tilde{b}} \tilde{\sigma}'_{xz} \, d\tilde{z} = \tilde{\mu}' \frac{\partial}{\partial \tilde{x}} \int_{-\tilde{H}}^{\varepsilon^4 \tilde{b}} \tilde{u}_z \, d\tilde{z} + O(\varepsilon^4) \quad (35) $$

$$ = \tilde{\mu}' \frac{\partial \tilde{q}_2}{\partial \tilde{x}} + O(\varepsilon^4) \quad (36) $$
We obtain

$$ \frac{\rho_s}{\rho_\ell} \frac{\partial \tilde{q}_2}{\partial \tilde{t}} - \frac{\partial \tilde{S}_{12}}{\partial \tilde{x}} = -\frac{1}{\varepsilon^2} \tilde{h} - \frac{1}{\varepsilon^2} \tilde{\sigma}_{zz}(-\tilde{H}) + \frac{\rho_s}{\rho_\ell} \tilde{\mu}' \frac{\partial^2 \tilde{q}_2}{\partial \tilde{x}^2} - \frac{\rho_s}{\rho_\ell} \tilde{\sigma}_{zz}'(-\tilde{H}) + O(\varepsilon^4) \quad (37) $$

where the expressions of $\tilde{\sigma}_{zz}(-\tilde{H})$ and $\tilde{\sigma}_{zz}'(-\tilde{H})$ have to be found.

## 2.4 Displacement at equilibrium

In a situation of equilibrium, defined by $\partial/\partial t = 0$ and $\partial/\partial x = 0$, Equation (22) reduces to

$$ \frac{d^2 u_z}{dz^2} = 0 \qquad (38) $$

With the boundary condition $u_z(-H) = 0$, the expression of $u_z$ is $u_z = K(z+H)$ where K is a constant. At the equilibrium, $b = \eta = 0$ and $h = h_0$. The solid/liquid interface is at $z = 0$ and the displacement at the interface is $u_z(0) = KH$. The dynamic boundary condition at the interface gives $\sigma_{zz}(0) = -\rho_\ell g h_0$. Since $\sigma_{zz} = (2\mu + \lambda)du_z/dz$, $\sigma_{zz}$ is a constant equal to $(2\mu + \lambda)K$, which gives the expression of K. The vertical displacement at equilibrium is thus

$$ u_{\text{zeq}} = -\frac{\rho_{\ell} g h_{0}}{2\mu + \lambda} (z + H) \qquad (39) $$

## 2.5 Non-equilibrium situation

In the following, we present a formal justification in the case of a non-equilibrium situation. We define the variable $u'_z = u_z - u_{\text{zeq}}$. Equation (22) can be written, in dimensional form,

$$ \frac{\partial^2 u'_z}{\partial t^2} - c_s^2 \frac{\partial^2 u'_z}{\partial x^2} - c_p^2 \frac{\partial^2 u'_z}{\partial z^2} - \nu' \frac{\partial^3 u'_z}{\partial x^2 \partial t} - \frac{\nu'_e}{2} \frac{\partial^3 u'_z}{\partial z^2 \partial t} = 0 \qquad (40) $$

where $c_s = \sqrt{\mu/\rho_s}$ is the celerity of the transverse waves, $c_p = \sqrt{(2\mu + \lambda)/\rho_s}$ is the celerity of the longitudinal waves (similar to P-waves), $\nu' = \mu'/\rho_s$ is a kinematic viscosity and

$$ \nu'_{e} = 2 \frac{2\mu' + \eta'}{\rho_{s}} \qquad (41) $$

is an effective kinematic viscosity.
We look for solutions of the form $u'_z = \Psi(z)\Phi(x, t)$, $\Psi(z)$ representing the variation of the vertical displacement with $z$. Equation (40) becomes

$$ \frac{1}{\Phi(x,t)} \left( \frac{\partial^2 \Phi}{\partial t^2} - c_s^2 \frac{\partial^2 \Phi}{\partial x^2} - \nu' \frac{\partial^3 \Phi}{\partial x^2 \partial t} \right) - \frac{c_p^2}{\Psi(z)} \frac{d^2 \Psi}{dz^2} - \frac{\nu'_e}{2} \frac{1}{\Phi(x,t)} \frac{\partial \Phi}{\partial t} \frac{1}{\Psi(z)} \frac{d^2 \Psi}{dz^2} = 0 \quad (42) $$

This equation can be written

$$ F(x,t) - G(z) - \frac{\nu'_e}{2c_p^2} E(x,t) G(z) = 0 \quad (43) $$

with

$$ F(x,t) = \frac{1}{\Phi(x,t)} \left( \frac{\partial^2 \Phi}{\partial t^2} - c_s^2 \frac{\partial^2 \Phi}{\partial x^2} - \nu' \frac{\partial^3 \Phi}{\partial x^2 \partial t} \right) \quad (44) $$

$$ G(z) = \frac{c_p^2}{\Psi(z)} \frac{d^2 \Psi}{dz^2} \quad (45) $$

and

$$ E(x,t) = \frac{1}{\Phi(x,t)} \frac{\partial \Phi}{\partial t} \quad (46) $$

Deriving (43) with respect to $x$ or $t$ gives respectively

$$ \frac{\partial F}{\partial x} - \frac{\nu'_e}{2c_p^2} \frac{\partial E}{\partial x} G(z) = 0 \quad (47) $$

and

$$ \frac{\partial F}{\partial t} - \frac{\nu'_e}{2c_p^2} \frac{\partial E}{\partial t} G(z) = 0 \quad (48) $$

These equations imply that $G(z)$ is a constant. Deriving (43) with respect to $z$ gives $dG/dz = 0$ which implies also that $G(z)$ is a constant. Taking a negative constant $-k^2$ we obtain

$$ \frac{d^2 \Psi}{dz^2} + k^2 \Psi = 0 \quad (49) $$

Due to the boundary condition $\Psi(z = -H) = 0$, the solution is of the form

$$ \Psi(z) = A \sin[k(z + H)] \quad (50) $$

where $A$ is a constant. The general solution is a superposition of these modes

$$ u'_z = \sum_i A_i \Phi_i(x, t) \sin[k_i(z + H)] \quad (51) $$
Since $\sigma_{zz} \simeq (2\mu + \lambda)\partial u_z/\partial z$, we have

$$ \sigma_{zz} \simeq -\rho_{\ell} g h_0 + (2\mu + \lambda) \sum_i A_i k_i \Phi_i \cos[k_i (z+H)] \quad (52) $$

This gives

$$ \sigma_{zz}(-H) = -\rho_{\ell} g h_0 + (2\mu + \lambda) \sum_i A_i k_i \Phi_i(x, t) \quad (53) $$

The order of magnitude of $k$ is $k \sim \omega/c_p$ where the angular wavenumber $\omega$ is related to the characteristic time of the tsunami $L/\sqrt{gh_0}$. Consequently

$$ kH \sim \varepsilon 2\pi \sqrt{\frac{\rho_s}{\rho_{\ell}}} \frac{\tilde{H}}{\sqrt{2\tilde{\mu} + \tilde{\lambda}}} \quad (54) $$

which implies that $kH = O(\varepsilon)$. This result can also be interpreted by noticing that $kH$ is also the ratio of the duration of the propagation of the longitudinal waves in the thickness of the solid layer to the characteristic time of the tsunami. The longitudinal waves propagates very quickly compared to the tsunami ($c_p$ is in the range 5000–8000 m·s⁻¹ whereas the celerity of a tsunami wave in the ocean is of the order of 200 m·s⁻¹) which means that $kH$ is small.

The boundary conditions at the solid/liquid interface yields $u'_z(z = b) = b$. The expression (51) of $u'_z$ can be written

$$ u'_z = \sum_i A_i \Phi_i(x, t) [k_i (z + H) + O(\varepsilon^3)] \quad (55) $$

Since $b/H = O(\varepsilon^4)$, we obtain

$$ u'_z(z = b) \simeq H \sum_i A_i k_i \Phi_i \quad (56) $$

which gives

$$ \sum_i A_i k_i \Phi_i \simeq \frac{b}{H} \quad (57) $$

and consequently

$$ \sigma_{zz}(-H) \simeq -\rho_{\ell} g h_0 + (2\mu + \lambda) \frac{b}{H} \quad (58) $$

Taking the derivative of the displacement with respect to $t$ gives the velocity

$$ \dot{u}_z = \sum_i A_i \frac{\partial \Phi_i}{\partial t} \sin[k_i (z+H)] \quad (59) $$
which reduces, neglecting terms of $O(\varepsilon^3)$ in the expansion of the sine function, to

$$ \dot{u}_z \approx (z + H) \sum_i A_i k_i \frac{\partial \Phi_i}{\partial t} \qquad (60) $$

The boundary condition (27) implies

$$ \frac{\partial b}{\partial t} \approx H \sum_i A_i k_i \frac{\partial \Phi_i}{\partial t} \qquad (61) $$

given that $b/H = O(\varepsilon^4)$. Calculating $q_2$ from its definition (30) leads to

$$ q_2 \approx \frac{H^2}{2} \sum_i A_i k_i \frac{\partial \Phi_i}{\partial t} \qquad (62) $$

using again $b/H = O(\varepsilon^4)$. The comparison of (61) and (62) yields

$$ \frac{\partial b}{\partial t} \approx \frac{2q_2}{H} \qquad (63) $$

Using (26), we obtain

$$ \sigma'_{zz} = (2\mu' + \eta') \sum_i A_i k_i \frac{\partial \Phi_i}{\partial t} \cos[k_i (z + H)] \qquad (64) $$

At the bottom of the solid layer

$$ \sigma'_{zz}(-H) = (2\mu' + \eta') \sum_i A_i k_i \frac{\partial \Phi_i}{\partial t} \qquad (65) $$

With (62), we obtain

$$ \sigma'_{zz}(-H) = (2\mu' + \eta') \frac{2q_2}{H^2} \qquad (66) $$

It is now possible to write the final form of Equation (37). In dimensional form, (37) writes

$$ \frac{\partial q_2}{\partial t} - \frac{1}{\rho_s} \frac{\partial S_{12}}{\partial x} = - \frac{\rho_\ell}{\rho_s} g h - \frac{\sigma_{zz}(-H)}{\rho_s} - \frac{\sigma'_{zz}(-H)}{\rho_s} + \nu' \frac{\partial^2 q_2}{\partial x^2} \qquad (67) $$

Replacing the components $\sigma_{zz}$ and $\sigma'_{zz}$ at $z = -H$ with their expressions (58) and (66) gives

$$ \frac{\partial q_2}{\partial t} - \frac{1}{\rho_s} \frac{\partial S_{12}}{\partial x} = - \frac{\rho_\ell}{\rho_s} g \eta - c_p^2 \frac{b}{H} - \nu'_e \frac{q_2}{H^2} + \nu' \frac{\partial^2 q_2}{\partial x^2} \qquad (68) $$

It should be noted that $\tilde{h} = 1 + \varepsilon^2 \tilde{\eta} - \varepsilon^3 \tilde{b}$ which means that the term $\rho_\ell g h_0$ disappears from this equation and that the term $\rho_\ell g b$ is of $O(\varepsilon)$ with respect to the other terms of the equation and consequently negligible.
## 2.6 Equations for the ocean layer

The equations for the liquid layer were derived in Richard (2021) [21]. A depth-averaged quantity over the liquid layer is defined for any quantity A by

$$ \langle A \rangle = \frac{1}{h} \int_{b}^{b+h} A \, dz \qquad (69) $$

Denoting by $R$ the depth-averaged density divided by the density at the free surface $\rho_0$, i.e. $R = \langle \rho \rangle / \rho_0$. The depth-averaged mass conservation equation for the liquid layer can be written

$$ \frac{\partial h R}{\partial t} + \frac{\partial h R U}{\partial x} = 0 \qquad (70) $$

In practice, $R$ depends mainly on a hydrostatic contribution $\rho_H$ on the density. Assuming a constant sound velocity $a$,

$$ \rho_H = \rho_0 e^{g(h+b-z)/a^2} \qquad (71) $$

The expression of $R$ is thus

$$ R = \frac{e^{M^2} - 1}{M^2} \approx 1 + \frac{M^2}{2} \qquad (72) $$

where $M$ is a Mach number defined by

$$ M = \frac{\sqrt{gh}}{a} \qquad (73) $$

As the water is weakly compressible, this Mach number is small and can be assumed to be of $O(\epsilon)$. The expression of the hydrostatic pressure is

$$ p_H = a^2 \rho_0 \left[ e^{g(h+b-z)/a^2} - 1 \right] \qquad (74) $$

The pressure at the bottom of the ocean is $p(b) \approx p_H(b)$ since the non-hydrostatic pressure is of $O(\epsilon^2)$ and consequently

$$ p(b) \approx \rho_0 ghR \qquad (75) $$

Since we write $p(b) = \rho_\ell gh$, we have

$$ \rho_\ell \approx \rho_0 R \qquad (76) $$

Equation (70) can be written (Richard 2021 [21])

$$ \frac{\partial h}{\partial t} + \frac{\partial h U}{\partial x} = \frac{M^2}{2} Q_0 h \frac{\partial U}{\partial x} \qquad (77) $$
with

$$Q_0 = \frac{2}{M^4} (e^{-M^2} + M^2 - 1) \approx 1 - \frac{M^2}{3} \qquad (78)$$

Since $h = h_0 + \eta - b$ and since $h_0$ is a constant, we can use (63) to write this equation

$$\frac{\partial \eta}{\partial t} + \frac{\partial hU}{\partial x} = \frac{M^2}{2} Q_0 h \frac{\partial U}{\partial x} + \frac{2q_2}{H} \qquad (79)$$

The depth-averaged momentum balance equation in the $Ox$-direction writes (see Richard 2021 [21] and noticing that the term in $p_N(b)\partial b/\partial x$ is of $O(\epsilon^6)$ and thus negligible)

$$\frac{\partial hRU}{\partial t} + \frac{\partial}{\partial x} \left( hRU^2 + Q_1 \frac{gh^2}{2} + hP \right) = -ghR \frac{\partial b}{\partial x} \qquad (80)$$

where $U = \langle v_x \rangle$ is the depth-averaged horizontal velocity, $P = \langle p_N \rangle / \rho_0$ is the depth-averaged non-hydrostatic pressure divided by $\rho_0$ ($p_N$ is the non-hydrostatic pressure in the water layer) and where

$$Q_1 = \frac{2}{M^4} (e^{M^2} - M^2 - 1) \approx 1 + \frac{M^2}{3} \qquad (81)$$

Noticing that

$$\frac{\partial}{\partial x} \left( Q_1 \frac{gh^2}{2} \right) = ghR \frac{\partial h}{\partial x} \qquad (82)$$

Equation (80) can be written

$$\frac{\partial hRU}{\partial t} + \frac{\partial}{\partial x} (hRU^2 + hP) = -ghR \frac{\partial \eta}{\partial x} \qquad (83)$$

The two other equations for the ocean layer are the same as in Richard (2021) [21]. The depth-averaged vertical velocity $W = \langle v_z \rangle$ satisfies the equation

$$\frac{\partial hRW}{\partial t} + \frac{\partial hRUW}{\partial x} = \frac{3}{2} P \qquad (84)$$

and the depth-averaged non-hydrostatic pressure satisfies

$$\frac{\partial hRP}{\partial t} + \frac{\partial hRU P}{\partial x} = -a^2 \left( 2W + h \frac{\partial U}{\partial x} \right) \qquad (85)$$
# 3 Model analysis

## 3.1 Hyperbolicity

The final system includes three equations for the solid layer and four equations for the ocean layer. These two subsystems are coupled by the variables $q_2$ and $\eta$. These equations write

$$ \frac{\partial \eta}{\partial t} + \frac{\partial hU}{\partial x} = \frac{M^2}{2} Q_0 h \frac{\partial U}{\partial x} + \frac{2q_2}{H} \quad (86) $$

$$ \frac{\partial hRU}{\partial t} + \frac{\partial}{\partial x} (hRU^2 + hP) = -ghR \frac{\partial \eta}{\partial x} \quad (87) $$

$$ \frac{\partial hRW}{\partial t} + \frac{\partial hRUW}{\partial x} = \frac{3}{2} P \quad (88) $$

$$ \frac{\partial hRP}{\partial t} + \frac{\partial hRU P}{\partial x} = -a^2 \left( 2W + h \frac{\partial U}{\partial x} \right) \quad (89) $$

$$ \frac{\partial q_2}{\partial t} - \frac{1}{\rho_s} \frac{\partial S_{12}}{\partial x} = - \frac{\rho_\ell}{\rho_s} g\eta - c_p^2 \frac{b}{H} - \nu'_e \frac{q_2}{H^2} + \nu' \frac{\partial^2 q_2}{\partial x^2} \quad (90) $$

$$ \frac{\partial S_{12}}{\partial t} - \mu \frac{\partial q_2}{\partial x} = 0 \quad (91) $$

$$ \frac{\partial b}{\partial t} = \frac{2q_2}{H} \quad (92) $$

In practice, the numerical simulations show that the diffusive term $\nu' \partial^2 q_2 / \partial x^2$ has an influence only on the propagation of transverse waves but has a wholly negligible influence on the propagation of tsunamis. On the other hand, the presence of this term increases considerably the computational time. In the following, we choose $\nu' = 0$, which reduces (90) to

$$ \frac{\partial q_2}{\partial t} - \frac{1}{\rho_s} \frac{\partial S_{12}}{\partial x} = - \frac{\rho_\ell}{\rho_s} g\eta - c_p^2 \frac{b}{H} - \nu'_e \frac{q_2}{H^2} \quad (93) $$

The resulting system can be written in matrix form

$$ \frac{\partial \mathbf{V}}{\partial t} + \mathbf{A} \frac{\partial \mathbf{V}}{\partial x} = \mathbf{S} \quad (94) $$

with $\mathbf{V} = (\eta, U, W, P, q_2, S_{12}, b)^\mathrm{T}$ and

$$ \mathbf{A} = \begin{bmatrix} U & hRe^{-M^2} & 0 & 0 & 0 & 0 & -U \\ \frac{P}{hR} + g & U & 0 & \frac{1}{R} & 0 & 0 & -\frac{P}{hR} \\ 0 & 0 & U & 0 & 0 & 0 & 0 \\ 0 & \frac{a^2}{R} & 0 & U & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -\frac{1}{\rho_s} & 0 \\ 0 & 0 & 0 & 0 & -\mu & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \quad (95) $$
The seven eigenvalues of A are

$$ \lambda_1 = 0 \qquad (96) $$

$$ \lambda_{2,3} = U \qquad (97) $$

$$ \lambda_{4,5} = U \pm \sqrt{(ghR + P)e^{-M^2} + a^2/R^2} \qquad (98) $$

$$ \lambda_{6,7} = \pm \sqrt{\frac{\mu}{\rho_s}} \qquad (99) $$

The eigenvalues are real. There are seven linearly independent eigenvectors. The system is thus hyperbolic. The characteristics (98) are relative to the liquid layer and include the sound velocity. They are discussed in Richard (2021) [21]. The characteristics (97) come also from the liquid layer and correspond to transport at the fluid velocity U. The characteristics (99) are new. They represent the velocities of the transverse waves in the solid layer. The characteristic (96) comes from Equation (92).

## 3.2 Equation of conservation of energy

Multiplying Equation (91) by $2S_{12}/\rho_s\mu H$ and (90) by $2q_2/H$ and adding these equations together lead to

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \frac{q_2^2}{H} + \frac{S_{12}^2}{\rho_s \mu H} \right) - \frac{\partial}{\partial x} \left( \frac{2}{\rho_s H} q_2 S_{12} \right) \\ & = - \frac{\rho_\ell}{\rho_s} g \eta \frac{2q_2}{H} - \frac{c_p^2}{H} b \frac{2q_2}{H} - 2 \frac{\nu_e'}{H^3} q_2^2 + 2 \frac{\nu'}{H} q_2 \frac{\partial^2 q_2}{\partial x^2} \end{aligned} \qquad (100) $$

Using (92), this equation can be further written

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \frac{q_2^2}{H} + \frac{S_{12}^2}{\rho_s \mu H} + \frac{c_p^2}{2H} b^2 \right) - \frac{\partial}{\partial x} \left( \frac{2}{\rho_s H} q_2 S_{12} \right) \\ & = - \frac{\rho_\ell}{\rho_s} g \eta \frac{\partial b}{\partial t} - 2 \frac{\nu_e'}{H^3} q_2^2 + 2 \frac{\nu'}{H} \frac{\partial}{\partial x} \left( q_2 \frac{\partial q_2}{\partial x} \right) - 2 \frac{\nu'}{H} \left( \frac{\partial q_2}{\partial x} \right)^2 \end{aligned} \qquad (101) $$

From Equation (87) and (86), we can derive the equation (see Richard 2021 [21])

$$ \begin{aligned} & \frac{\partial}{\partial t} \left[ hR \left( \frac{U^2}{2} + \frac{Q_2}{R} \frac{gh}{2} + gb \right) \right] \\ & + \frac{\partial}{\partial x} \left[ hRU \left( \frac{U^2}{2} + \frac{Q_2}{R} \frac{gh}{2} + gb \right) + \left( Q_1 \frac{gh^2}{2} + hP \right) U \right] \\ & = hP \frac{\partial U}{\partial x} + ghR \frac{\partial b}{\partial t} \end{aligned} \qquad (102) $$
with

$$Q_2 = \frac{2}{M^4} \left[ 1 + (M^2 - 1) e^{M^2} \right] = 1 + \frac{2}{3} M^2 + O(M^4) \quad (103)$$

Multiplying Equation (89) by $P/a^2$ and (88) by $4W/3$ and adding them together leads to

$$\frac{\partial}{\partial t} \left( \frac{2}{3} h R W^2 + h R \frac{P^2}{2a^2} \right) + \frac{\partial}{\partial x} \left( \frac{2}{3} h R U W^2 + h R U \frac{P^2}{2a^2} \right) = -h P \frac{\partial U}{\partial x} \quad (104)$$

Adding this equation with (102) gives

$$\frac{\partial h R e_L}{\partial t} + \frac{\partial}{\partial x} (h R U e_L + \Pi U) = g h R \frac{\partial b}{\partial t} \quad (105)$$

with

$$e_L = \frac{U^2}{2} + \frac{2}{3} W^2 + \frac{Q_2}{R} \frac{gh}{2} + gb + \frac{P^2}{2a^2} \quad (106)$$

and

$$\Pi = Q_1 \frac{gh^2}{2} + hP \quad (107)$$

Multiplying Equation (101) by $\rho_s$ and Equation (105) by $\rho_0$ and adding them together gives, using $\rho_\ell = \rho_0 R$,

$$\begin{aligned} & \frac{\partial}{\partial t} \left( \rho_0 h R e_L + \frac{\rho_s}{H} q_2^2 + \frac{S_{12}^2}{\mu H} + \frac{\rho_s c_p^2}{2H} b^2 - \rho_\ell g h_0 b + \rho_\ell g \frac{b^2}{2} \right) \\ & + \frac{\partial}{\partial x} \left( \rho_0 h R U e_L + \rho_0 \Pi U - \frac{2}{H} q_2 S_{12} \right) = -2 \frac{\rho_s \nu'_e}{H^3} q_2^2 \\ & + \frac{\partial}{\partial x} \left( 2 \frac{\rho_s \nu'_e}{H} q_2 \frac{\partial q_2}{\partial x} \right) - 2 \frac{\rho_s \nu'_e}{H} \left( \frac{\partial q_2}{\partial x} \right)^2 \end{aligned} \quad (108)$$

We can add to the left-hand side of this equation the term

$$\frac{\partial}{\partial t} \left[ \frac{1}{2} \frac{H}{\rho_s c_p^2} (\rho_\ell g h_0)^2 \right] \quad (109)$$

which is equal to zero. Using (39), we can define a potential energy $E_{ps}$ and write

$$E_{ps} = \frac{\rho_s c_p^2}{2H} b^2 - \rho_\ell g h_0 b + \frac{1}{2} \frac{H}{\rho_s c_p^2} (\rho_\ell g h_0)^2 = \frac{1}{2} \frac{\rho_s c_p^2}{H} (b + u_{\text{zeq}}(0))^2 \quad (110)$$

where $u_{\text{zeq}}(0)$ is the equilibrium displacement at the solid / liquid interface. This expression is similar to the potential energy of a spring (the well-known
expression $(1/2)k(x - x_0)^2)$ with a spring constant $\rho_s c_p^2/H = (2\mu + \lambda)/H$ and where the position of the interface is $b$ with ocean loading and $-u_{zeq}(0)$ without loading (i.e. without water). Denoting by $E_{ps}$ this potential energy, the final equation of energy balance is

$$ \frac{\partial}{\partial t} \left( \rho_0 h R e_L + \frac{\rho_s}{H} q_2^2 + \frac{S_{12}^2}{\mu H} + E_{ps} + \rho_\ell g \frac{b^2}{2} \right) + \frac{\partial}{\partial x} \left( \rho_0 h R U e_L + \rho_0 \Pi U - \frac{2}{H} q_2 S_{12} \right) = -2 \frac{\rho_s \nu_e'}{H^3} q_2^2 + \frac{\partial}{\partial x} \left( 2 \frac{\rho_s \nu'}{H} q_2 \frac{\partial q_2}{\partial x} \right) - 2 \frac{\rho_s \nu'}{H} \left( \frac{\partial q_2}{\partial x} \right)^2 \quad (111) $$

The terms $\rho_s q_2^2/H$ and $S_{12}^2/(\mu H)$ in the energy are due respectively to the kinetic energy and to the strain energy of the elastic transverse waves. The corresponding energy flux for these elastic waves is $-2q_2 S_{12}/H$. The term $\rho_\ell g b^2/2$ in the energy is due to the fact that we neglected $\rho_s g b/\rho_\ell$ in (68). This term is very small but we are looking for an exact energy conservation equation for the system (86)–(92), so we cannot neglect it here. The terms corresponding to the elastic energy do not appear in the flux because we are in the framework of linear elasticity. The viscosity entails the presence of dissipative terms and of a diffusive term. If $\nu' = 0$, there is no diffusive term and only one dissipative term due to $\nu'_e$.

Thus, the system (86)–(92) admits an exact equation of energy balance. The energy is conserved in the absence of viscosity, which, as expected, introduces dissipation and diffusion in the system.

### 3.3 Dispersive properties

The dispersive properties are studied on the linearized system. The equations are linearized around the equilibrium state defined by $h = h_0$, $\eta = b = 0$, $U = W = 0$, $P = 0$, $q_2 = 0$ and $S_{12} = 0$.

Let's denote by $(\eta', U', W', P', q'_2, S'_{12}, b')^T$ the small perturbations around this equilibrium state. They are written $A \exp[i(kx - \omega t)]$. Inserting these expressions into the system of equations (86)–(89), (93) and (91)–(92) yields
<figure><img src="image_3.png" /><figcaption>Figure 2: Solutions of the dispersion relation (112): Gravity branch (blue curve), acoustic branch (black curve) and elastic branch (red curve) with the values $g = 9.8 \text{ m} \cdot \text{s}^{-2}$, $\rho_s = 3375 \text{ kg} \cdot \text{m}^{-3}$, $\rho_\ell = 1000 \text{ kg} \cdot \text{m}^{-3}$, $\mu = 6.7 \cdot 10^{10} \text{ Pa}$, $\lambda = 8.2 \cdot 10^{10} \text{ Pa}$, $a = 1500 \text{ m} \cdot \text{s}^{-1}$, $h_0 = 4000 \text{ m}$ and $H = 220 \text{ km}$. (a) Frequency as a function of the wavenumber in the case $\nu_e = 0$ (k is real). (b) Frequency as a function of the imaginary part of the wavenumber in the case $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$.</figcaption></figure>

the dispersion relation

$$ \left( \omega^2 - c_s^2 k^2 - \frac{2c_p^2}{H^2} + i\omega \frac{\nu_e'}{H^2} \right) \left\{ \frac{h_0^2 R_0^2}{3a^2} \omega^4 \right. \\ \left. - \omega^2 \left[ 1 + \frac{k^2 h_0^2}{3} \left( 1 + \frac{gh_0}{a^2} R_0^3 e^{-M_0^2} \right) \right] + gk^2 h_0 R_0 e^{-M_0^2} \right\} \\ - \frac{2\rho_\ell g}{\rho_s H} \left( \frac{h_0^2 R_0^2}{3a^2} \omega^2 - 1 - \frac{k^2 h_0^2}{3} \right) \omega^2 = 0 \quad (112) $$

where $M_0^2 = gh_0/a^2$. This equation admits three solution branches (see figure 2(a) where $\nu_e = 0$). The first branch corresponds to the gravity waves (in blue on figure 2(a)). If the seafloor is supposed to be perfectly rigid and if the seawater is supposed to be incompressible ($a \to \infty$), only this branch exists. In this case, the system reduces to the Serre-Green-Naghdi equations (Serre 1953 [22], Green & Naghdi 1976 [9]) and the dispersion relation (112) reduces to

$$ \omega^2 = gh_0 \frac{k^2}{1 + k^2 h_0^2 / 3} \quad (113) $$

The second branch (in black on figure 2(a)) is the acoustic branch. If the seafloor is perfectly rigid and if the water is weakly compressible, only the
gravity branch and the acoustic branch exist. The dispersion relation (112) reduces to

$$ \frac{h_0^2 R_0^2}{3a^2} \omega^4 - \omega^2 \left[ 1 + \frac{k^2 h_0^2}{3} \left( 1 + \frac{g h_0}{a^2} R_0^3 e^{-M_0^2} \right) \right] + g k^2 h_0 R_0 e^{-M_0^2} = 0 \quad (114) $$

which was studied in Richard (2021) [21]. This relation is consistent to the dispersion relation of the linear theory of compressible fluids at the long wave limit. Note that the gravity branch is modified by compressible effects even at the long wave limit.

The third branch (in red on figure 2(a)) is due to the elasticity of the sea floor and can be called the elastic branch. The gravity branch is also modified by elastic effects, even at the long wave limit. The curves on figure 2(a) were obtained with the following parameters: $g = 9.8 \text{ m} \cdot \text{s}^{-2}$, $\rho_s = 3375 \text{ kg} \cdot \text{m}^{-3}$, $\rho_\ell = 1000 \text{ kg} \cdot \text{m}^{-3}$, $\mu = 6.7 \cdot 10^{10} \text{ Pa}$, $\lambda = 8.2 \cdot 10^{10} \text{ Pa}$, $a = 1500 \text{ m} \cdot \text{s}^{-1}$, $h_0 = 4000 \text{ m}$, $H = 220 \text{ km}$ and $\nu_e = 0$. The wavenumber and the angular frequency are represented in dimensionless form. They are scaled as

$$ \tilde{k} = k h_0 ; \quad \tilde{\omega} = \omega \sqrt{\frac{h_0}{g}} \qquad (115) $$

The viscosity $\nu_e$ induces a damping, which is characterized by an imaginary part of the wavenumber $k$. In the case of figure 2(a), $\nu_e = 0$ and $k$ is real. The imaginary part is shown on figure 2(b) where the same parameters as for figure 2(a) where used except that the viscosity is $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$. The imaginary parts of the gravity branch (in blue) and the acoustic branch (in black) are always equal to zero, but the imaginary part of the elastic branch is not zero. This means that, as it can be expected, the viscosity in the solid layer entails a damping effect only for the elastic waves.

The viscosity $\nu_e$ has no effect on the gravity and the acoustic branches. The effect of the viscosity on the elastic branch is shown on figure 3. The same parameters as above are used except for the viscosity which is $\nu_e = 0$ (black curve), $\nu_e = 5 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$ (blue curve), $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$ (green curve), $\nu_e = 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (red curve) and $\nu_e = 5 \cdot 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (orange curve). The real part of $\tilde{k}$ is shown on figure 3(a) and the imaginary part is shown on figure 3(b). There is a cutoff frequency if $\nu_e = 0$ ($\omega$ has a non-zero minimum value for $k = 0$) but there is no cutoff frequency if the viscosity is not zero. For small values of the viscosity, the elastic branch differs from the case $\nu_e = 0$ only for very small values of $\tilde{k}$. The effect of the viscosity on the elastic branch is more and more pronounced when the viscosity is higher. As expected, the imaginary part, and thus the damping, increases with the viscosity.
<figure><img src="image_4.png" /><figcaption>Figure 3: Effect of the viscosity $\nu_e$ on the elastic branch. Same parameters as for figure 2 with $\nu_e = 0$ (black curve), $\nu_e = 5 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$ (blue curve), $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$ (green curve), $\nu_e = 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (red curve) and $\nu_e = 5 \cdot 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (orange curve). (a) Frequency as a function of the real part of the wavenumber. (b) Frequency as a function of the imaginary part of the wavenumber.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Dimensionless phase velocity as a function of the real part of the dimensionless wavenumber. (a) Gravity branch (blue curve), acoustic branch (black curve) and elastic branch (red curve) with the same parameters as for figure 2 and $\nu_e = 0$ (k is real). (b) Elastic branch with the same parameters and $\nu_e = 0$ (black curve), $\nu_e = 5 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$ (blue curve), $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$ (green curve), $\nu_e = 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (red curve) and $\nu_e = 5 \cdot 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (orange curve).</figcaption></figure>
The dimensionless phase velocity $\tilde{v}_{\varphi}$ as a function of the dimensionless wavenumber is shown on figure 4. The phase velocity is scaled as

$$ \tilde{v}_{\varphi} = \frac{v_{\varphi}}{\sqrt{gh_0}} \qquad (116) $$

The phase velocity for the three branches is presented on figure 4(a) with the gravity branch (blue), the acoustic branch (black) and the elastic branch (red). The parameters are the same as above with $\nu_e = 0$. The phase velocity for the acoustic and elastic branches goes to infinity when $k \to 0$ but of course the group velocity goes to zero at the same limit. This is a consequence of the existence of a cutoff frequency.

The viscosity $\nu_e$ has no effect on the phase velocity of the gravity and acoustic branches but has a strong effect on the phase velocity of the elastic branch, which is presented on figure 4(b). The parameters are again the same as above except for the viscosity which is $\nu_e = 0$ (black curve), $\nu_e = 5 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$ (blue curve), $\nu_e = 5 \cdot 10^{10} \text{ m}^2 \cdot \text{s}^{-1}$ (green curve), $\nu_e = 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (red curve) and $\nu_e = 5 \cdot 10^{11} \text{ m}^2 \cdot \text{s}^{-1}$ (orange curve). With a nonzero viscosity the phase velocity reaches a finite value for $k = 0$ because there is no more cutoff frequency. The phase velocity of the elastic branch decreases if $k$ increases if $\nu_e$ is smaller than a critical value, which is here close to $5 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$, and increases with $k$ if $\nu_e$ is larger than this critical value. The limit value of the phase velocity of the elastic branch when $k \to \infty$ is the celerity of the transverse waves $\sqrt{\mu/\rho_s}$ whatever the viscosity $\nu_e$.

For the propagation of tsunamis, the branch of main interest is the gravity branch. The dimensionless phase velocity of the gravity branch is presented on figure 5 as a function of the dimensionless wavenumber. On figure 5(a) and (b) the black curve is the phase velocity given by Airy linear theory of incompressible waves on a perfectly rigid seafloor. This theory leads to

$$ \tilde{v}_{\varphi} = \sqrt{\frac{\tanh \tilde{k}}{\tilde{k}}} \qquad (117) $$

For $\tilde{k} < 0.3$, the shallow water approximation is accurately satisfied and the dispersion relation (113) of the Serre-Green-Naghdi equations gives the same phase velocity with an excellent accuracy. Its expression is

$$ \tilde{v}_{\varphi} = \sqrt{\frac{1}{1 + \frac{\tilde{k}^2}{3}}} \qquad (118) $$

The blue curve on figure 5(a) and (b) is the phase velocity with compressible effects but a perfectly rigid seafloor. It corresponds to the dispersion relation
<figure><img src="image_6.png" /><figcaption>Figure 5: Phase velocity of the gravity branch as a function of the wavenumber. Black curve: incompressible case with a perfectly rigid seafloor. Blue curve: compressible case with a perfectly rigid seafloor. Red, orange and green curves: compressible case with an elastic seafloor. (a) Same parameters as for figure 2 except $H = 410$ km (red curve), $H = 220$ km (orange curve) and $H = 100$ km (green curve). (b) Same parameters as for figure 2 except $\mu = 6.7 \cdot 10^{10}$Pa and $\lambda = 8.2 \cdot 10^{10}$Pa (red curve), $\mu = 4.4 \cdot 10^{10}$Pa and $\lambda = 4.6 \cdot 10^{10}$Pa (orange curve), $\mu = 2.7 \cdot 10^{10}$Pa and $\lambda = 3.4 \cdot 10^{10}$Pa (green curve).</figcaption></figure>

(114). The reduction of the phase velocity from the black curve to the blue curve is due to compressibility alone, since there is no elastic effect. This blue curve is consistent with the linear theory of compressible fluids with an excellent accuracy (Richard 2021 [21]) since $\bar{k}$ is small (smaller than 0.3). The other curves show the effect of elasticity. Since compressibility is included, the reduction of the phase velocity due to elastic effects is the difference between the blue curve and the red, orange or green curves.

On figure 5(a), the parameters are the same as above with $\nu_e = 0$ except that $H = 410$ km (red curve), $H = 220$ km (orange curve) and $H = 100$ km (green curve). The effect of the thickness $H$ of the solid layer, the coefficients of Lamé being the same, is complex. For relatively large values of the wavenumber, the reduction of the phase velocity due to elasticity is larger if $H$ is smaller but for the smallest values of the wavenumber, the reduction of the phase velocity is larger if $H$ is larger. Moreover, if $H$ is large enough, the phase velocity decreases if the wavenumber decreases beyond some value. This is the reverse dispersion effect observed for tsunamis (Watada et al. 2014 [28]) i.e. the tsunami speed becomes slower for tsunamis of long periods (beyond 1000 s). In the examples of figure 5(a), the period beyond which there is a reverse dispersion is between 1200 and 1600 s. The reverse dispersion
sion effect is thus obtained qualitatively with the right order of magnitude. To obtain a quantitatively accurate result, a model with more than one solid layer is probably necessary (see below).

The effect of the rigidity of the seafloor is studied on figure 5(b). The parameters are the same as above with $ν_e = 0$ and $H = 220$ km but $µ = 6.7 \cdot 10^{10}$Pa and $λ = 8.2 \cdot 10^{10}$Pa (red curve), $µ = 4.4 \cdot 10^{10}$Pa and $λ = 4.6 \cdot 10^{10}$Pa (orange curve), $µ = 2.7 \cdot 10^{10}$Pa and $λ = 3.4 \cdot 10^{10}$Pa (green curve). The reduction of the phase velocity due to elasticity is larger for all wave numbers if the rigidity of the solid layer is smaller i.e. if the coefficients of Lamé are smaller.

The model is able to capture the reverse dispersion effect. However, the results above show that tsunamis of relatively small wavelength are influenced mainly by a superficial Earth layer, particularly the crust, with a relatively small thickness and smaller rigidity, whereas the tsunamis of larger wavelength are influenced by a larger thickness of Earth, including part of the upper mantle, in spite of its larger rigidity. For practical applications, a model with at least two solid layers, for example one for the oceanic crust and another one for a part of the upper mantle, could be necessary to obtain accurate results.

## 3.4 Long wave limit

At the long wave limit $k \to 0$ and for a weakly compressible fluid, the phase velocity of the gravity branch is

$$v_{\varphi} \approx \sqrt{\frac{gh_0}{1 + \frac{\rho_{\ell}gH}{\lambda + 2\mu}}} \left(1 - \frac{M_0^2}{4}\right) \quad (119)$$

to within terms of $O(M_0^4)$. The case of an incompressible seawater corresponds to $M_0 = 0$ and a perfectly rigid seafloor is found at the limit $µ \to \infty$ and $λ \to \infty$. This gives the incompressible and perfectly rigid case, $v_{\varphi} = \sqrt{gh_0}$, the case of a weakly compressible seawater and a perfectly rigid seafloor $v_{\varphi} = \sqrt{gh_0}(1 - M_0^2/4)$ (see Richard 2021 [21]) and the case of an incompressible seawater and an elastic seafloor

$$v_{\varphi} \approx \sqrt{\frac{gh_0}{1 + \frac{\rho_{\ell}gH}{\lambda + 2\mu}}} \quad (120)$$
This long wave limit can be also found when the right-hand side of Equation (93) is completely relaxed (the viscosity has no influence here) i.e. if

$$b = - \frac{\rho_{\ell} g H}{(\lambda + 2\mu)} \eta \qquad (121)$$

At this limit, this relation replaces Equation (93) and Equation (91) becomes useless. Using Equation (92), Equation (86) can be rewritten

$$\left(1 + \frac{\rho_{\ell} g H}{\lambda + 2\mu}\right) \frac{\partial \eta}{\partial t} + \frac{\partial h U}{\partial x} = \frac{M^2}{2} Q_0 h \frac{\partial U}{\partial x} \qquad (122)$$

At the long wave limit $k \to 0$, the non-hydrostatic pressure, and thus the dispersion, can be neglected. This implies that $P = 0$ and that Equations (88)–(89) become useless. Equation (87) reduces to

$$\frac{\partial h R U}{\partial t} + \frac{\partial h R U^2}{\partial x} = -g h R \frac{\partial \eta}{\partial x} \qquad (123)$$

The system (122)–(123) is non-dispersive and its phase velocity is

$$v_{\varphi} = \sqrt{\frac{gh_0 R_0 e^{-M_0^2}}{1 + \frac{\rho_{\ell} g H}{\lambda + 2\mu}}} \qquad (124)$$

Since $\sqrt{R_0 \exp(-M_0^2)} = 1 - M_0^2/4 + O(M_0^4)$, this phase velocity reduces to (119).

If, in the system (122)–(123), we further neglect the seawater compressibility ($M_0 = 0$) and if the nonlinearity is small enough that the nonlinear convective term in (123) becomes negligible and that we can write $h \simeq h_0$ in the right-hand side of (123), we find the system

$$\left(1 + \frac{\rho_{\ell} g H}{\lambda + 2\mu}\right) \frac{\partial \eta}{\partial t} + \frac{\partial h U}{\partial x} = 0 \qquad (125)$$

$$\frac{\partial h U}{\partial t} = -g h_0 \frac{\partial \eta}{\partial x} \qquad (126)$$

which is the system of Inazu & Saito (2013) [12] (note the difference of notations: their $\eta$ is for us $\eta - b$ and their $\eta_0$ is for us $-b$). Inazu & Saito (2013) [12] introduced a constant parameter $\beta$ such that, translated into our notations,

$$b = -\beta (\eta - b) \qquad (127)$$
The comparison of this relation with (121) gives the expression of $\beta$

$$ \beta = \frac{\rho_{\ell} g H}{\lambda + 2\mu - \rho_{\ell} g H} \qquad (128) $$

With the various numerical values used above, we calculate values of $\beta$ in the range 0.01–0.02 which are of the same order of magnitude as the values of $\beta$ found by Inazu & Saito (2013) [12] which are in the range 0.005–0.030 and most often 0.015 or 0.020.

The model of Inazu & Saito (2013) [12] is thus the long wave limit of our model when the dispersion, the seawater compressibility and all non-linearities are neglected and when the elastic part is completely relaxed i.e. Equation (93) reduces to (121).

# 4 Numerical simulations

## 4.1 Numerical scheme

The equations (86)–(89) are treated as in Richard (2021) [21] with a splitting method. The viscoelastic part of the system (90)–(92) is solved explicitly. Only the acoustic part is treated implicitly. In the explicit step, the following system is solved:

$$ \frac{\partial \eta}{\partial t} + \frac{\partial hU}{\partial x} = \frac{M^2}{2} Q_0 h \frac{\partial U}{\partial x} + \frac{2q_2}{H} \qquad (129) $$

$$ \frac{\partial hRU}{\partial t} + \frac{\partial hRU^2}{\partial x} = -\frac{\partial hP}{\partial x} - ghR \frac{\partial \eta}{\partial x} \qquad (130) $$

$$ \frac{\partial hRW}{\partial t} + \frac{\partial hRUW}{\partial x} = 0 \qquad (131) $$

$$ \frac{\partial hRP}{\partial t} + \frac{\partial hRUP}{\partial x} = 0 \qquad (132) $$

$$ \frac{\partial q_2}{\partial t} - \frac{1}{\rho_s} \frac{\partial S_{12}}{\partial x} = - \frac{\rho_{\ell}}{\rho_s} g \eta - c_p^2 \frac{b}{H} - \nu_e' \frac{q_2}{H^2} + \nu' \frac{\partial^2 q_2}{\partial x^2} \qquad (133) $$

$$ \frac{\partial S_{12}}{\partial t} - \mu \frac{\partial q_2}{\partial x} = 0 \qquad (134) $$

$$ \frac{\partial b}{\partial t} = \frac{2q_2}{H} \qquad (135) $$

In the implicit step, the system is

$$ \frac{\partial \eta}{\partial t} = 0 \qquad (136) $$
$$ \frac{\partial U}{\partial t} = 0 \qquad (137) $$

$$ \frac{\partial hRW}{\partial t} = \frac{3}{2}P \qquad (138) $$

$$ \frac{\partial hRP}{\partial t} = -a^2 \left( 2W + h \frac{\partial U}{\partial x} \right) \qquad (139) $$

$$ \frac{\partial q_2}{\partial t} = 0 \qquad (140) $$

$$ \frac{\partial S_{12}}{\partial t} = 0 \qquad (141) $$

$$ \frac{\partial b}{\partial t} = 0 \qquad (142) $$

which also implies that $\partial h/\partial t = 0$ and $\partial R/\partial t = 0$ in the implicit step.

A second-order scheme was implemented. The second order in space was obtained with a monotone upstream-centred scheme for conservation laws (MUSCL). There is no need for a limiter since there is no shock. The second order in time is obtained with the diagonally-implicit Runge-Kutta (DIRK) Implicit-Explicit (IMEX) ARS2(2,2,2) scheme of Ascher et al. (1997) [5] (using the notation of Pareschi & Russo 2005 [19]). The ARS2(2,2,2) has two explicit stages and two implicit stages.

The implicit procedure used for the fast subsystem is a backward Euler method. The derivative $\partial U/\partial x$ is calculated by a central finite difference method. There is thus no global system to solve. The expressions of W and P are obtained explicitly, which means that the implicit steps are computationally very cheap.

A finite volume method (Godunov-type) with a Rusanov Riemann solver is used for the explicit steps, except that the term $\partial(hP)/\partial x$ is treated as a source term by a central finite difference. For Equations (129)–(132), as in Richard (2021) [21], the interface speed of the Rusanov solver is calculated as the maximum value between the left and right cells of the hydraulic characteristics $U \pm \sqrt{gh}$. The sound velocity and the velocity of the transverse elastic waves are not included in the Riemann solver for these equations, which reduces considerably the numerical diffusion of the scheme. For Equations (133)–(134), the interface speed is the velocity $c_s$ of the transverse elastic waves since these two equations govern the propagation of these waves.

In the second explicit stage of the IMEX scheme, both systems, (129)–(135) and (136)–(142), have to be solved (see Ascher et al. 1997 [5] and Richard 2021 [21]). In this case, the system (136)–(142) is treated explicitly by a forward Euler method.
The time step is calculated by a standard Courant-Friedrichs-Levy (CFL) condition $\Delta t = C\Delta x/c_s$ where C is the Courant number, $\Delta t$ the time step and $\Delta x$ the cell size. The largest characteristic velocity is the velocity of the transverse elastic waves $c_s$. A Courant number equal to 0.8 was used in the computations.

The diffusive term in (133) is discretized by a central finite difference method. This term reduces considerably the time step since the diffusive CFL condition is in $(\Delta x)^2$. This term has an effect on the propagation of transverse elastic waves. However the effect of this term on the propagation of the tsunami is indiscernible. In the following, the value $\nu' = 0$ is taken – i.e. the diffusive term is removed – and the hyperbolic CFL condition holds.

## 4.2 Arrival delay and leading negative phase

The model is used to simulate an academic 1D-tsunami in order to check whether the main physical features are predicted with a reasonable order of magnitude. In an ocean of constant depth $h_0 = 4000$ m, a tsunami is initiated with an initial condition on the surface elevation $\eta(x, t)$ following the sinusoidal law

$$ \eta(x, 0) = \frac{A}{2} \left[ 1 + \cos \left( \pi \frac{x - x_0}{\ell} \right) \right] \quad (143) $$

if $x_0 - \ell \le x \le x_0 + \ell$ and $\eta(x, 0) = 0$ if not. The initial elevation is a sinusoid centred on $x = x_0$ with a width $\ell$ and a maximum amplitude $A$.

The length of the computational domain is $L = 8000$ km. At both ends of the domain additional cells are used to dampen the waves as in sponge layers. The tsunami is initiated at $x_0 = 1000$ km. The various quantities are registered at $x = 7500$ km after a propagation over a distance of 6500 km. In all simulations, the sound velocity is $a = 1500 \text{ m} \cdot \text{s}^{-1}$. Different values of the Lamé coefficients are used to assess the effect of the Earth elasticity on the tsunami propagation. The three following sets of coefficients are considered, from the most rigid to the most elastic:

$$ \lambda = 8.2 \cdot 10^{10} \text{ Pa}; \mu = 6.7 \cdot 10^{10} \text{ Pa} \quad (144) $$

$$ \lambda = 4.6 \cdot 10^{10} \text{ Pa}; \mu = 4.4 \cdot 10^{10} \text{ Pa} \quad (145) $$

$$ \lambda = 3.4 \cdot 10^{10} \text{ Pa}; \mu = 2.7 \cdot 10^{10} \text{ Pa} \quad (146) $$

Three different values of the initial length $\ell$ are tested: $\ell = 400$ km; $\ell = 200$ km; $\ell = 100$ km. The value of the initial amplitude is $A = 10$ m. The other parameters are $g = 9.8 \text{ m} \cdot \text{s}^{-2}$, $H = 220$ km, $\rho_\ell = 1000 \text{ kg} \cdot \text{m}^{-3}$, $\rho_s = 3375 \text{ kg} \cdot \text{m}^{-3}$, $\nu'_e = 5.0 \cdot 10^9 \text{ m}^2 \cdot \text{s}^{-1}$ and $\nu' = 0$.
<figure><img src="image_7.png" /><figcaption>Figure 6: Water elevation as a function of time: IR model (black curve); CR model (blue curve); IE model (green curve); CE model (red curve).</figcaption></figure>

For comparison, the model with an incompressible seawater and rigid Earth (IR model), the model with an incompressible seawater and elastic Earth (IE model) and the model with a compressible seawater and rigid Earth (CR model) are also solved. If the Earth is supposed to be rigid, $q_2 = 0$ and the equations (90)–(92) are not solved. If the seawater is supposed to be incompressible, $R = Q_0 = 1$ and $M = 0$.

Figure 6 shows the variations with time of the water elevation $\eta$ after a propagation of 6500 km for the IR model (black curve), the CR model (blue curve), the IE model (green curve) and the compressible and elastic model (red curve) i.e. the full model (CE model). When elasticity was included, the Lamé parameters were taken as in (144). The length of the initial perturbation is $l = 200$ km.

With the IR model the wave maximum arrives at $t = 32\,824.9$ s which gives an average velocity of $198\,\text{m·s}^{-1}$, which is almost equal to $\sqrt{gh_0}$. In the IE case, the wave maximum arrives with a delay of 73.6 s. In the CR case, the wave maximum arrives with a delay of 143.3 s. With the full CE model, the wave maximum arrives with a delay of 217.3 s, which is very close to the sum 216.9 s of the delay of the IE model and of the delay of the CR model. The delay in the full model represents 0.66% of the total propagation time.
Given that the most rigid case was considered, the order of magnitude of the arrival time delay is correct.

Figure 6 shows that there is a small negative water elevation before the arrival of the main hump for the models including elasticity (green and red curves). Without Earth elasticity, there is no negative elevation before the main wave. This negative elevation corresponds to the leading negative phase observed before the arrival of the first tsunami positive peak. The minimum elevation of this negative phase is almost the same (−14.7 cm) for the IE model and for the CE model. This negative phase is thus not due to the sea-water compressibility but only to the Earth elasticity. At the beginning of the propagation, there is no such leading negative phase. The initial perturbation has no negative elevation. The leading negative phase appears during the propagation and becomes important at long distances from the initiation. It is not related to the source of the tsunami but only to propagation effects related to the elasticity of the solid layer.

Moreover the dispersive tail of the tsunami (the train of small waves after the main hump) has a higher elevation if elasticity is included in the model (green and red curves) than for a rigid Earth (black and blue curves). This is in a way the opposite effect to the leading negative phase. Before the main wave, elasticity decreases the water elevation and after the main wave, elasticity increases the water elevation. This last effect is less obvious because of the oscillations of the dispersive tail.

The maximum amplitude is the same for both rigid cases, incompressible or compressible, (4.988 m), but it is 6.9 cm smaller if elasticity is included, incompressible or compressible. Figure 6 shows that compressibility entails only a delay of the wave without distorsion (the blue curve is only delayed compared to the black curve and the red curve is only delayed compared to the green curve) whereas Earth elasticity entails both a delay and a distorsion of the wave form.

All these results are in agreement with previous works (in particular Watada *et al.* 2014 [28], Allgeyer & Cummins 2014 [4]). The model is able to capture, qualitatively and with a correct order of magnitude, the main effects of Earth elasticity, the arrival time delay and the leading negative phase.

Figure 7 presents the variations, in space at a time $t = 32.965.5$ s, of the water elevation $\eta$ (black curve), the position of the solid/liquid interface $b$ (red curve) and the depth-averaged quantities $q_2$ (blue curve) and $S_{12}$ (green curve) for the full compressible and elastic model with the previous parameters. The position of the solid/liquid interface is always negative and reaches a minimum of −3.2 cm, which seems to be a correct order of magnitude. The sign of $q_2$ is negative in the front part of the wave and positive in the rear part. Given that $2q_2/H$ is a source term in the evolution equation for the
<figure><img src="image_8.png" /><figcaption>Figure 7: Variations of the water elevation $\eta$ (in m, black curve), the solid/liquid interface $b$ (in cm, red curve), the depth-averaged quantities $q_2$ (blue curve) and $S_{12}$ (green curve) as a function of the abscissa $x$ at time $t = 32 965.5$ s.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Water elevation as a function of time. (a) $l = 200$ km. The Lamé parameters are given by (144) (black curve), (145) (blue curve) and (146) (red curve). (b) $l = 400$ km (blue curve); $l = 200$ km (black curve); $l = 100$ km (red curve). The Lamé parameters are given by (144).</figcaption></figure>

water elevation $\eta$ (see Equation (86)), a negative value of $q_2$ tends to decrease the water elevation and a positive value of $q_2$ tends to increase the water elevation. Since $q_2$ becomes negative well before the arrival of the main hump (in time), this explains the leading negative phase. Likewise, $q_2$ is positive well after the main hump (in time), in the dispersive tail, explaining that the dispersive tail has a higher elevation if elasticity is included in the model.

The effect of the choice of the Lamé parameters are presented on Figure 8(a). The water elevation as a function of time is plotted in the same conditions as for Figure 6 except that the model is the full model with compressibility and elasticity and that the Lamé parameters are given by (144) (black curve, the most rigid), by (145) (blue curve) and by (146) (red curve, the least rigid). The results show that the less rigid the solid layer, the deeper the leading negative phase and the higher the dispersive tail is. The arrival time delay is also larger if the solid layer is less rigid. In the least rigid case (146), the arrival delay is 315.5 s, which is 0.95 % of the total propagation time. The maximum diminution of depth in the negative leading phase is –39.9 cm instead of –14.7 cm in the most rigid case and the amplitude of the main hump decreases by 26.3 cm instead of 6.9 cm.

The effect of the length $l$ of the initial perturbation is shown on Figure 8(b). The conditions are the same as for Figure 6 except that the model is the full model with compressibility and elasticity and that $l = 400$ km (blue curve), $l = 200$ km (black curve) and $l = 100$ km (red curve).

The waves of the dispersive tail are larger if the initial perturbation is shorter. If an initial perturbation more complicated than a sinusoid is used, the dispersive tail can be even larger. This is due to the fact that more
frequencies are included in the initial perturbation. The depth of the leading
negative phase depends also on the initial perturbation. In the present case,
it is smaller for the longest perturbation $\ell = 400$ km.

# 5 Conclusion

In order to derive a model for tsunamis propagation including the effect of seafloor elasticity, the ocean was modelled as a liquid layer above a solid layer of constant thickness. The system is 2D leading to a 1D-model. The solid is modelled with the constitutive law of linear viscoelasticity with constant Lamé parameters and density. Below this solid layer, the Earth is assumed to be perfectly rigid. The liquid layer is modelled as in Richard (2021) [21] and includes seawater compressibility.

The resulting equations are integrated over the thickness of the solid layer for the solid part and are averaged over the depth of the liquid layer for the liquid part. Since in practice the diffusive terms can be neglected, the obtained model is hyperbolic. Its characteristics include the sound velocity and the velocity of the transverse elastic waves. The system admits an exact equation of energy conservation, the viscosity of the solid layer giving dissipative and diffusive terms.

The system is dispersive. The solutions of the dispersion relation correspond to three branches: an acoustic branch as in Richard (2021) [21], an elastic branch and a gravity branch, which is the most important for tsunamis propagation. The water compressibility and the seafloor elasticity have an effect on the phase velocity of the gravity branch even at the long wave limit. The elasticity of the solid layer entails a reverse dispersion effect i.e. the phase velocity reaches a maximum for some value of the wave number and decreases, not only for large wave numbers, but also for smaller wave numbers.

At the long wave limit, when dispersion, compressibility and all nonlinearities are neglected, and if the elastic part is supposed to be completely relaxed, our model reduces to the model of Inazu & Saito (2013) [12].

The equations are solved by an explicit finite-volume method except the acoustic part, which is solved implicitly as in Richard (2021) [21]. However there is no global system to solve and the implicit step has thus a very cheap computational cost. There is no need to include neither the sound velocity nor the velocity of the transverse elastic waves in the interface speed of the Riemann solver for the equations governing the liquid layer, which means that the numerical diffusion of the scheme is very small for the propagation of tsunamis. The velocity of the transverse elastic waves is taken into account
in the Riemann solver for the equations related to the solid layer.

The model was solved to simulate a 1D-tsunami initiated with a sinusoidal perturbation of the water elevation. Seawater compressibility produces an arrival time delay as in Richard (2021) [21]. The elasticity of the solid layer gives an additional delay and produces also a leading negative phase, i.e. a negative water elevation before the main hump, which is absent if the solid layer is assumed to be rigid. The elevation of the dispersive tail of the tsunami is higher if the seafloor elasticity is taken into account. The seawater compressibility entails only a delay of the wave without distorsion whereas the elasticity of the solid layer produces both a delay and a distorsion of the wave form.

All these effects are in agreement with previous works (in particular Watada et al. 2014 [28], Allgeyer & Cummins 2014 [4]). The model derived in this article is thus a promising alternative to the existing approaches. It is able to capture all physical features of tsunamis propagation with a local and hyperbolic system of equations, which has a relatively simple mathematical structure and which is easy to solve numerically.

For practical applications, the model should be extended to the 2D-case with bathymetry, variations of gravity and Coriolis forces, using spherical coordinates. This work suggests also that, in order to reach a sufficient accuracy, two solid layers could be considered, one for the oceanic crust, which is relatively softer, and one for the more rigid upper mantle.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding

This work was supported by the French National Research Agency project NABUCO [grant ANR-17-CE40-0025].

References

[1] Ali Abdolali, Usama Kadri, and James T Kirby. Effect of water compressibility, sea-floor elasticity, and field gravitational potential on tsunami phase speed. *Sci. Rep.*, 9(1):1–8, 2019.

[2] Ali Abdolali and James T Kirby. Role of compressibility on tsunami propagation. *J. Geophys. Res.-Oceans*, 122(12):9780–9794, 2017.
[3] Y Accad and Chaim L Pekeris. Solution of the tidal equations for the m2 and s2 tides in the world oceans from a knowledge of the tidal potential alone. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 290(1368):235-266, 1978.

[4] Sébastien Allgeyer and Phil Cummins. Numerical tsunami simulation including elastic loading and seawater density stratification. *Geophys. Res. Lett.*, 41(7):2368-2375, 2014.

[5] Uri M Ascher, Steven J Ruuth, and Raymond J Spiteri. Implicit-explicit runge-kutta methods for time-dependent partial differential equations. *Appl. Numer. Math.*, 25(2-3):151-167, 1997.

[6] Toshitaka Baba, Sebastien Allgeyer, Jakir Hossen, Phil R Cummins, Hiroaki Tsushima, Kentaro Imai, Kei Yamashita, and Toshihiro Kato. Accurate numerical simulation of the far-field tsunami caused by the 2011 tohoku earthquake, including the effects of boussinesq dispersion, seawater density stratification, elastic loading, and gravitational potential change. *Ocean Model.*, 111:46-54, 2017.

[7] Toshitaka Baba, Yodai Gon, Kentaro Imai, Kei Yamashita, Tetsuo Matsuno, Mitsuru Hayashi, and Hiroshi Ichihara. Modeling of a dispersive tsunami caused by a submarine landslide based on detailed bathymetry of the continental slope in the nankai trough, southwest japan. *Tectono-physics*, 768:228182, 2019.

[8] Adam M Dziewonski and Don L Anderson. Preliminary reference earth model. *Physics of the earth and planetary interiors*, 25(4):297-356, 1981.

[9] AE Green and PM Naghdi. A derivation of equations for wave propagation in water of variable depth. *J. Fluid Mech.*, 78(2):237-246, 1976.

[10] Stephan T Grilli, Jeffrey C Harris, Tayebeh S Tajalli Bakhsh, Timothy L Masterlark, Christodoulos Kyriakopoulos, James T Kirby, and Fengyan Shi. Numerical simulation of the 2011 tohoku tsunami based on a new transient fem co-seismic source: Comparison to far-and near-field observations. *Pure Appl. Geophys.*, 170(6):1333-1359, 2013.

[11] Tung-Cheng Ho, Kenji Satake, and Shingo Watada. Improved phase corrections for transoceanic tsunami data in spatial and temporal source estimation: Application to the 2011 tohoku earthquake. *Journal of Geophysical Research: Solid Earth*, 122(12):10-155, 2017.
[12] Daisuke Inazu and Tatsuhiko Saito. Simulation of distant tsunami propagation with a radial loading deformation effect. *Earth, Planets and Space*, 65(8):835-842, 2013.

[13] Teruyuki Kato, Yukihiro Terada, Hitoyoshi Nishimura, Toshihiko Nagai, and Shun'ichi Koshimura. Tsunami records due to the 2010 chile earthquake observed by gps buoys established along the pacific coast of japan. *Earth, Planets and Space*, 63(6):e5-e8, 2011.

[14] James T Kirby, Fengyan Shi, Babak Tehranirad, Jeffrey C Harris, and Stephan T Grilli. Dispersive tsunami waves in the ocean: Model equations and sensitivity to dispersion and coriolis effects. *Ocean Model.*, 62:39-55, 2013.

[15] Randall J LeVeque et al. *Finite volume methods for hyperbolic problems*, volume 31. Cambridge university press, 2002.

[16] IM Longman. A green's function for determining the deformation of the earth under surface mass loads: 1. theory. *Journal of Geophysical Research*, 67(2):845-850, 1962.

[17] Emile A Okal. Mode-wave equivalence and other asymptotic problems in tsunami theory. *Physics of the Earth and Planetary interiors*, 30(1):1-11, 1982.

[18] Spiros D Pagiatakis. The response of a realistic earth to ocean tide loading. *Geophysical Journal International*, 103(2):541-560, 1990.

[19] Lorenzo Pareschi and Giovanni Russo. Implicit-explicit runge-kutta schemes and applications to hyperbolic systems with relaxation. *J. Sci. Comput.*, 25(1):129-155, 2005.

[20] Alexander B Rabinovich, Philip L Woodworth, and Vasily V Titov. Deep-sea observations and modeling of the 2004 sumatra tsunami in drake passage. *Geophysical Research Letters*, 38(16), 2011.

[21] Gaël L Richard. An extension of the boussinesq-type models to weakly compressible flows. *European Journal of Mechanics-B/Fluids*, 89:217-240, 2021.

[22] François Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. *La Houille Blanche*, (6):830-872, 1953.
[23] Liujuan Tang, Vasily V Titov, Eddie N Bernard, Yong Wei, Christopher D Chamberlin, Jean C Newman, Harold O Mofjeld, Diego Arcas, Marie C Eble, Christopher Moore, et al. Direct energy estimation of the 2011 japan tsunami using deep-ocean pressure measurements. *Journal of Geophysical Research: Oceans*, 117(C8), 2012.

[24] Victor C Tsai, Jean-Paul Ampuero, Hiroo Kanamori, and David J Stevenson. Estimating the effect of earth elasticity and variable water density on tsunami speeds. *Geophys. Res. Lett.*, 40(3):492-496, 2013.

[25] Jean Virieux. P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method. *Geophysics*, 51(4):889-901, 1986.

[26] Kelin Wang, Yan Hu, and Jiangheng He. Deformation cycles of subduction earthquakes in a viscoelastic earth. *Nature*, 484(7394):327-332, 2012.

[27] Steven N Ward. Relationships of tsunami generation and an earthquake source. *Journal of Physics of the Earth*, 28(5):441-474, 1980.

[28] Shingo Watada, Satoshi Kusumoto, and Kenji Satake. Traveltime delay and initial phase reversal of distant tsunamis coupled with the self-gravitating elastic earth. *J. Geophys. Res. Sol. Ea.*, 119(5):4287-4310, 2014.

[29] Masahiro Yoshimoto, Shingo Watada, Yushiro Fujii, and Kenji Satake. Source estimate and tsunami forecast from far-field deep-ocean tsunami waveforms—the 27 february 2010 mw 8.8 maule earthquake. *Geophysical Research Letters*, 43(2):659-665, 2016.

[30] David A Yuen and WR Peltier. Normal modes of the viscoelastic earth. *Geophysical Journal International*, 69(2):495-526, 1982.