# Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible Euler equation

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, François Bouchut, Jacques Sainte-Marie

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

## To cite this version:

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, François Bouchut, Jacques Sainte-Marie. Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible Euler equation. Journal of Fluid Mechanics, 2023, 970, pp.A28. 10.1017/jfm.2023.595. hal-03880423v4

HAL Id: hal-03880423

https://hal.science/hal-03880423v4

Submitted on 2 Aug 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible Euler equation

Juliette Dubois<fnref n="1" />*, Sébastien Imperiale<fnref n="2" />, Anne Mangeney<fnref n="3" />,<fnref n="1" />, François Bouchut<fnref n="4" />, Jacques Sainte-Marie<fnref n="1" />

August 2, 2023

## Abstract

In this paper we construct an accurate linear model describing the propagation of both acoustic and gravity waves in water. This original model is obtained by the linearization of the compressible Euler equations, written in Lagrangian coordinates. The system is studied in the isentropic case, with a free surface, an arbitrary bathymetry, and vertical variations of the background temperature and density. We show that our model is an extension of some models from the literature to the case of a non-barotropic fluid with a variable sound speed. Other models from the literature are recovered from our model through two asymptotic analysis, one for the incompressible regime and one for the acoustic regime. We also propose a method to write the model in Eulerian coordinates. Our model includes many physical properties, such as the existence of internal gravity waves or the variation of the sound speed with depth.

## 1 Introduction

Several authors have proposed to use the propagation of acoustic waves in the ocean to detect tsunamis, as the sound travels in water at approximately 1500 m s⁻<fnref n="1" /> and the velocity of tsunami wave is around 300 m s⁻<fnref n="1" /> (Constantin, 2009). The existence of hydro-acoustic signals generated by tsunami sources such as earthquakes or landslides was shown in Tolstoy (1950). This motivates

<footnote n="1">¹Inria, Team ANGE, Inria Paris, 2 rue Simone Iff, 75012 Paris and Laboratoire Jacques-Louis Lions, Sorbonne Université, 75005 Paris</footnote>

<footnote>*Corresponding author: juliette.dubois@inria.fr</footnote>

<footnote n="2">²Inria, M3DISIM Team, Inria Saclay, 1 rue Honoré d'Estienne d'Orves, 91120 Palaiseau</footnote>

<footnote n="3">³Seismology Group, Institut de Physique du Globe de Paris, Université Paris Diderot, Sorbonne Paris Cité, 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote n="4">⁴Laboratoire d'Analyse et de Mathématiques Appliquées, CNRS, Université Gustave Eifel, Marne-la-Vallée, France</footnote>
the mathematical modelling of the propagation of both surface waves – the tsunami – and underwater acoustic waves, also called hydroacoustic waves, in a compressible formulation.

The idea of using acoustic-gravity waves for tsunami early-warning systems dates back from 1950 (Ewing *et al.*, 1950). A more recent study (Stiassnie, 2010) indicates that the pressure variations induced by the tsunami are significant enough to be used for the improvement of the tsunami early-warning systems.

For the description of the propagation of sound in water, the most common model is a linear wave equation for the fluid potential, i.e. for an irrotational fluid (Jensen *et al.*, 2011). When both surface and acoustic waves are considered, different types of models are available. In his work, Stiassnie (2010) studies the acoustic equation for the fluid potential coupled with a free-boundary condition. The three-dimensional acoustic equation is analyzed by Nosov & Kolesov (2007) and a depth-integrated version is proposed by Sammarco *et al.* (2013) in order to reduce the computational costs. This approach was further developed in a series of papers (Cecioni *et al.*, 2014; Abdolali *et al.*, 2015; Gomez & Kadri, 2021).

Another approach was proposed by Longuet-Higgins (1950) where the equation is still on the fluid potential, but includes a gravity term. This equation, including second order terms, made it possible for the first time to explain the seismic noise generated worldwide by waves interactions in the ocean (Stutzmann *et al.*, 2012). This model was also the starting point of an extensive work to describe the nonlinear interactions between acoustic and gravity waves (Kadri & Stiassnie, 2013). In other works, such as Smith (2015) and Auclair *et al.* (2021), the flow is not assumed irrotational, so that the equations are written for the fluid velocity. They include gravity terms and a vertical stratification for the background density, temperature and salinity. This generalisation allows to study the internal waves caused by the stratification of the fluid, and dispersion relations for the three types of waves (acoustic, internal, surface) are obtained.

The above cited works share one or several of the following assumptions: irrotational flow, homogeneous background density or barotropic fluid, and a constant speed of sound. These modelling choices have a strong influence on the structure of the equations, resulting in a variety of tools for their analysis and their numerical approximation. For example, the irrotationality assumption allows to reduce the number of unknowns, but the validity of this assumption in the compressible case is not clear. Furthermore, in the models that do not assume irrotational flow, the bed is assumed to be flat, even though bed variations are a key element impacting tsunami and acoustic wave propagation (Caplan-Auerbach *et al.*, 2014). In the ocean, the choice of a constant sound speed may be not appropriate since the variation of the sound speed creates the SOFAR channel, a horizontal strip in which the acoustic waves propagate with very little energy loss. Quantifying the impact of these approximations requires the use of simulations based on a more complete model. Finally, the free-surface equation induces a strong nonlinearity in the system. Indeed, the domain on which the equation are written depends on the solution to the equations. The common
approach for the linearization of the system consists in writing the linear equa-
tions on the unperturbed domain. However, by doing so an approximation on
the domain is made, in addition to the approximation made on the solution.
It is not clear how to quantify the magnitude of the error made by the two
combined approximations.

The aim of this work is to adress these modelling choices by deriving an
accurate linear model as rigorously as possible with only very few assumptions
for hydro-acoustic, internal and surface waves propagating in a fluid over an
arbitrary bathymetry. Salinity, thermal dissipation and viscosity are neglected,
and in order to linearize the equations we assume that the ocean is at equilibrium
and at rest before the earthquake or landslide occurs, and that the tsunami
source induces a small displacement of the water. In this model, the speed
of sound results from the imposed background temperature profile, so that the
effects of the SOFAR channel on the propagation of the hydroacoustic waves are
naturally present. The obtained model is comparable to the model of Auclair
et al. (2021), however our model includes a bathymetry and a variable sound
speed. Moreover our approach differs on several aspects:

* The problem is formulated as a second-order equation, which allows the use of numerical solver dedicated to wave propagation problem such as Specfem (Komatitsch & Tromp, 1999). Specfem uses spectral finite elements to compute acoustic and/or elastic wave propagations, and is widely used in the seismology community, as for example to simulate seismic waves generated by landslides (Kuehnert et al., 2020). In addition to the acoustic waves already modelled in Specfem, the model proposed in the present paper includes the linear water waves.

* The method used to write the linearization of a free-surface flow is generic, and can be applied to extend the model. A possible extension would include second-order terms (a similar work was done in Longuet-Higgins (1950) in the barotropic case). Another possibility is to take into account the interaction with the Earth. In particular, one can consider the elastic deformations of the ocean bottom, that are shown to impact the travel time of tsunami waves (Abdolali et al., 2019).

Another advantage of having a model with few assumptions is that a cas-
cade of simplified systems can be obtained from it. We indeed show that with
some simplifying assumptions, our model reduces to the models proposed in the
litterature. The analysis of these simplifications helps to understand the math-
ematical and physical choices made in these models. For example, the most
common model for the propagation of hydro-acoustic waves (Stiassnie, 2010;
Nosov & Kolesov, 2007; Sammarco et al., 2013) is recovered from the proposed
model by assuming a barotropic fluid and a constant background density.

We also show that our model and the simplified models are energy-preserving.
Our model is a linear version of the Euler equations, and the equation account-
ing for the energy conservation may be modified by the linearization. To ensure
that the obtained model is physically relevant, we check that an equation for
the energy conservation holds in the linear case. Beyond this aspect, the energy preservation allows to write stable numerical schemes (Allaire, 2015). Indeed, the properties of a numerical scheme are often related to the preservation of a discrete energy. For these reasons, the energy preservation is a key feature, both in the continuous and in the discrete level.

The paper is organized as follows. In Sec. 2 the compressible Euler equations for a free-surface flow are written, then the system is transformed in Lagrangian coordinates in order to keep an exact description of the free surface. After linearization, a wave-like equation for the fluid velocity is obtained and we show that the energy of the system is preserved. In Sec. 3, we show that with additional assumptions, the model reduces to other linear models from the literature. The barotropic case is studied, then the incompressible limit and the acoustic limit of the wave equation are written. In Sec. 4, we present a method allowing to write the model in Eulerian coordinates. The obtained system can be linearized at the cost of an additional approximation, namely that the equations have to be restricted to a fixed domain, and we show how to obtain a linear free-surface condition. Finally, in Sec 5 we obtain a dispersion relation which includes all of the physical effects mentioned above. In particular, it is a generalization of the dispersion relation studied in the work of Auclair et al. (2021) to the case of a varying sound speed.

## 2 Linearization of compressible Euler equations in Lagrangian coordinates

We derive here a linear model around a state at rest for the isentropic compressible Euler equation with a free surface and an arbitrary bathymetry, valid for a generic equation of state and a generic vertical temperature profile. We aim at deriving a model which is physically relevant in the sense that it preserves or dissipate energy. For this reason, we will analyse the energy equation associated to this system and show that preservation or dissipation of energy requires a condition on the fluid stratification that is related to the internal waves.

We consider a portion of the ocean away from the coast and at equilibrium: there is no mean current and the temperature varies only vertically. In this work, we do not take the presence of salinity into account, hence the ocean is assimilated to pure water. The bottom and the surface of the domain are assumed to be parametrized as graphs, respectively the topography $z_b(x, y) \ge 0$ and the free-surface elevation $\eta(x, y, t)$. The reference level $z = 0$ is situated inside the earth at an arbitrary level. The ground displacement induced by an earthquake or landslide source is assumed to take place away from the coast, so that the domain is considered infinite in the $(x, y)$ plane, see Figure 1. The domain is assumed to have the following description, for all time $t$,

$$ \Omega(t) = \{(x, y, z) \in \mathbb{R}^3 \mid z_b(x, y) < z < \eta(x, y, t)\}. \quad (1) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: The domain $\Omega(t)$, (1a) for time $t = 0$, (1b) for time $t > 0$. In (1a) typical profiles for the temperature and the density at rest are drawn.</figcaption></figure>

The boundary of the domain are then defined by

$$ \Gamma_s(t) = \{(x, y, z) \in \mathbb{R}^3 | z = \eta(x, y, t)\}, \quad (2) $$

and

$$ \Gamma_b(t) = \{(x, y, z) \in \mathbb{R}^3 | z = z_b(x, y) - b(x, y, t)\}. \quad (3) $$

The function $b$ is the source term, namely the normal displacement of the seabed. It can represent for example an earthquake or a landslide. It is assumed that this displacement starts at a time $t_0 > 0$, so that $b(x, y, 0) = 0$.

## 2.1 Euler equation in Eulerian coordinate

### 2.1.1 Equations in the volume

The unknowns are the fluid velocity $U$, its density $\rho$, its pressure $p$, its temperature $T$, its internal energy $e$ and its entropy $s$.

For future reference the equations are written for a viscous fluid with thermal dissipation. The stress tensor of a Newtonian fluid $\mathcal{T}$ has the form

$$ \mathcal{T} = (-p + \lambda \nabla \cdot U)\mathbb{I} + 2\mu \mathbb{D}(U), \quad (4) $$

Where $\mathbb{D}(U)$ is defined by $\mathbb{D}(U) = (\frac{1}{2}(\partial_i U^j + \partial_j U^i))_{i,j=x,y,z}$ and $\mathbb{I}$ is the identity matrix in $\mathbb{R}^3$. The heat flux is denoted by $\mathbf{q}$ and is a function of $\rho$ and $T$.

The conservation of mass, momentum and energy of a Newtonian fluid read,
in the domain $\Omega(t)$,

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{U}) = 0, \qquad (5) $$

$$ \frac{\partial}{\partial t}(\rho \mathbf{U}) + \nabla \cdot (\rho \mathbf{U} \otimes \mathbf{U}) + \nabla p = \rho \mathbf{g} + \nabla(\lambda \nabla \cdot \mathbf{U}) + \nabla \cdot (2\mu \mathbb{D}(\mathbf{U})), \qquad (6) $$

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \rho \frac{|\mathbf{U}|^2}{2} + \rho e \right) + \nabla \cdot \left( \left( \rho \frac{|\mathbf{U}|^2}{2} + \rho e + p \right) \mathbf{U} \right) \\ & = \rho \mathbf{g} \cdot \mathbf{U} + \nabla \cdot (\lambda \mathbf{U} \nabla \cdot \mathbf{U}) + \nabla \cdot (2\mu \mathbb{D}(\mathbf{U}) \cdot \mathbf{U}) - \nabla \cdot \mathbf{q}. \end{aligned} \qquad (7) $$

The acceleration of gravity is $\mathbf{g} = -g \mathbf{e}_3$ with $g > 0$ and $\mathbf{e}_3$ is the unit vector in the vertical direction, oriented upwards.

In order to describe the acoustic waves, we derive an equation for the pressure. Among $(\rho, e, T, p, s)$ only two variables are independent because of the Gibbs law and of the equation of state (Gill, 1982). When considering $\rho$ and $s$ as independent, it is natural to introduce the scalar functions $f_e$, $f_p$ and $f_T$ satisfying

$$ e = f_e(\rho, s), \quad p = f_p(\rho, s), \quad T = f_T(\rho, s). \qquad (8) $$

With the Gibbs law $(\partial f_e / \partial \rho = f_p / \rho^2$ and $\partial f_e / \partial s = f_T)$, one has

$$ \frac{\partial e}{\partial t} + \mathbf{U} \cdot \nabla e = \frac{f_p}{\rho^2} \left( \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho \right) + f_T \left( \frac{\partial s}{\partial t} + \mathbf{U} \cdot \nabla s \right). \qquad (9) $$

Using (7) – $\mathbf{U} \cdot (6)$ and (9), one obtains as an intermediate step the evolution equation of the entropy,

$$ \rho T \left( \frac{\partial s}{\partial t} + \mathbf{U} \cdot \nabla s \right) = \lambda (\nabla \cdot \mathbf{U})^2 + 2\mu \mathbb{D}(\mathbf{U}) : \mathbb{D}(\mathbf{U}) - \nabla \cdot \mathbf{q}. \qquad (10) $$

Now, since $p = f_p(\rho, s)$ we have

$$ \frac{\partial p}{\partial t} + \mathbf{U} \cdot \nabla p = \frac{\partial f_p}{\partial \rho} \left( \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho \right) + \frac{1}{T \rho} \frac{\partial f_p}{\partial s} \left( \frac{\partial s}{\partial t} + \mathbf{U} \cdot \nabla s \right), \qquad (11) $$

hence using (5) and (10) one obtains

$$ \frac{\partial p}{\partial t} + \mathbf{U} \cdot \nabla p = - \frac{\partial f_p}{\partial \rho} (\rho \nabla \cdot \mathbf{U}) + \frac{1}{T \rho} \frac{\partial f_p}{\partial s} (\lambda (\nabla \cdot \mathbf{U})^2 + 2\mu \mathbb{D}(\mathbf{U}) : \mathbb{D}(\mathbf{U}) - \nabla \cdot \mathbf{q}). \qquad (12) $$

At this point we use the common assumption that the viscous term and the thermal dissipation can be neglected compared to the advection term (see Lannes, 2013, Chap. 1). With (10) we see that this is equivalent to assuming that the flow is isentropic. Moreover, from physical considerations the function $f_p$ must satisfy $\partial_\rho f_p(\rho, s) \ge 0$, hence we can introduce the speed of sound $c$ defined by

$$ c^2 = \frac{\partial f_p}{\partial \rho}(\rho, s). \qquad (13) $$
The equation (12) then reads

$$ \frac{\partial p}{\partial t} + \mathbf{U} \cdot \nabla p + \rho c^2 \nabla \cdot \mathbf{U} = 0. \qquad (14) $$

Eq. (14) is used for the study of a compressible fluid in the isentropic case, see Gill (1982), Chap. 4. Note that the speed of sound c can also be viewed as a function of p and T, in that case we have

$$ c^2(p, T) = \frac{\partial f_p}{\partial \rho} (f_\rho(p, T), f_s(p, T)). \qquad (15) $$

In practice we choose to work directly with the expression $c = c(p, T)$ tabulated in IAPWS-SR7 (2009). Note that here the temperature intervenes as a side variable, because it is necessary to compute the speed of sound. However, we will see later that only the temperature profile of the state at rest is needed to close the system.

### 2.1.2 Boundary conditions

The following boundary conditions hold,

$$ \left\{ \begin{array}{ll} \mathbf{U} \cdot \mathbf{n}_b = u_b = \partial_t b & \text{on } \Gamma_b, \\ p = p^a & \text{on } \Gamma_s. \end{array} \right. \qquad (16) $$

$$ p = p^a \quad \text{on } \Gamma_s. \qquad (17) $$

The bottom boundary condition (16) is a non-penetration condition with a source term. It models the tsunami source as a displacement of the ocean bottom with velocity $u_b$. We denote by $\mathbf{n}_b$ the unit vector normal to the bottom and oriented outwards. The second condition (17) is a dynamic condition, where we assume that the surface pressure is at equilibrium with a constant atmospheric pressure $p^a$. Note that the elevation $\eta$ is solution of the following kinematic equation

$$ \frac{\partial \eta}{\partial t} + \mathbf{U} \cdot \begin{pmatrix} \partial_x \eta \\ \partial_y \eta \\ -1 \end{pmatrix} = 0 \quad \text{on } \Gamma_s(t). \qquad (18) $$

### 2.1.3 Initial conditions and equilibrium state

It is assumed that the initial state corresponds to the rest state, meaning that $\eta(x, y, 0) = H$ with the elevation at rest $H$ being independent of space and $H > z_b(x, y)$, therefore

$$ \Omega(0) = \{(x, y, z) \in \mathbb{R}^3 \mid z_b(x, y) < z < H\}. \qquad (19) $$

We choose the following initial conditions for the velocity, the temperature, the density and the pressure,

$$ \mathbf{U}(x, y, z, 0) = 0, \qquad (20) $$

$$ T(x, y, z, 0) = T_0(z), \quad \rho(x, y, z, 0) = \rho_0(z), \quad p(x, y, z, 0) = p_0(z). \qquad (21) $$
Where $T_0, \rho_0, p_0$ are functions defined on $(0, H)$ but because of the topography $z_b(x, y)$, the functions $T, \rho, p$ need not to be defined from $z=0$ for all $(x, y)$.

When the source term $u_b$ vanishes, we have an equilibrium state around $U \equiv 0$ if the functions $T_0, \rho_0, p_0$ satisfy

$$ \nabla p_0 = \rho_0 g, \quad \rho_0 = f_\rho(p_0, T_0), \quad p_0(H) = p^a. \qquad (22) $$

Hence, if $T_0(z)$ is given, the system

$$ \left\{ \begin{array}{l} \frac{dp_0}{dz} = -g f_\rho(p_0, T_0), \quad z \in (0, H), \\ p_0 = p^a, \quad z = H. \end{array} \right. \qquad (23) $$

$$ p_0 = p^a, \quad z = H. \qquad (24) $$

can be solved to compute $p_0$, and then $\rho_0$ is computed with $\rho_0 = f_\rho(p_0, T_0)$. Note that, in the forthcoming sections the system (5), (6), (14) with boundary conditions (16), (17) and initial conditions (20), (21) will be linearized around the previously defined equilibrium state.

## 2.2 Lagrangian description

Although most of the works on free-surface flows are done in Eulerian coordinates, the Lagrangian formalism is sometimes preferred, see for example the paper from Nouguier et al. (2015) and the references therein, or the work of Godlewski et al. (1999) for a precise derivation of linear models. Here we choose the Lagrangian description in order to avoid any approximation on the shape of the domain when we linearize the equations. The usual approximation made on the surface for the linear models in Eulerian coordinates consists in evaluating the surface condition on pressure at a fixed height, rather than at the actual, time-dependent free surface. The kinematic boundary condition is also replaced by its linear approximation. For the derivation and justification of the approximation, see Lighthill (1978, Chap. 3).

Let $\hat{\Omega}$ be the domain of the ocean at a reference time, with its surface boundary $\hat{\Gamma}_s$ and bottom boundary $\hat{\Gamma}_b$. The reference time is chosen before the tsunami generation, so that the surface of the domain is horizontal. In fact the following natural choice is made

$$ \hat{\Omega} = \Omega(0), \quad \hat{\Gamma}_s = \Gamma_s(0), \quad \hat{\Gamma}_b = \Gamma_b(0). \qquad (25) $$

The position at the reference time of a fluid particle is denoted

$$ \boldsymbol{\xi} = (\xi^1, \xi^2, \xi^3) \in \hat{\Omega}. \qquad (26) $$

At time $t$, the fluid has moved, the domain is $\Omega(t)$ and the new position of a fluid particle is $\mathbf{x} = (x(\xi, t), y(\xi, t), z(\xi, t)) \in \Omega(t)$. We denote by $\phi$ the transformation from $\hat{\Omega}$ to $\Omega(t)$ that maps each particle from its reference position $\xi$ to its position $\mathbf{x}$ at time $t$ (see figure 2).

$$ \phi: \left\{ \begin{array}{l} \hat{\Omega} \rightarrow \Omega(t) \\ \boldsymbol{\xi} \mapsto \mathbf{x}(\boldsymbol{\xi}, t) \end{array} \right. \qquad (27) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: The mapping $\phi_t$ between the reference domain $\hat{\Omega} = \Omega(0)$ and the domain $\Omega(t)$.</figcaption></figure>

Hence one has $\mathbf{x} = \phi(\xi, t)$. The transformation is assumed invertible, in particular we do not consider the case of wave breaking. We also define the displacement of the fluid. For each fluid particle with initial position $\xi$, its displacement is defined by

$$ \mathbf{d}(\xi, t) = \phi(\xi, t) - \xi. \qquad (28) $$

The gradient of $\phi$ with respect to $\xi$ is denoted $F$,

$$ F = \nabla_{\xi}\phi, \qquad (29) $$

and its determinant is denoted $J$. Both $F$ and $J$ can be expressed as functions of the displacement,

$$ F = \mathbb{I} + \nabla_{\xi}\mathbf{d}, \quad J = \det F, \qquad (30) $$

where $\nabla_{\xi}$ is the gradient with respect to the coordinate system $\xi$. For a function $X(\mathbf{x}, t)$ defined on the domain $\Omega(t)$, we introduce $\hat{X}(\xi, t)$ defined on $\hat{\Omega}$ by

$$ \hat{X}(\xi, t) = X(\phi(\xi, t), t). \qquad (31) $$

Finally, note that the velocity $\hat{U}(\xi, t) = U(\phi(\xi, t), t)$ is the time derivative of the displacement $\mathbf{d}$.

$$ \hat{U} = \frac{\partial \mathbf{d}}{\partial t}. \qquad (32) $$

With this change of coordinates, the system (5), (6), (14) is now defined in the time-independent reference domain $\hat{\Omega}$ and it reads

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \hat{\rho}}{\partial t} + \frac{\hat{\rho}}{|J|} \nabla_{\xi} \cdot (|J| F^{-1} \hat{U}) = 0, \\ \\ \displaystyle \hat{\rho} \frac{\partial \hat{U}}{\partial t} + F^{-T} \nabla_{\xi} \hat{p} = \hat{\rho} g, \\ \\ \displaystyle \frac{\partial \hat{p}}{\partial t} + \frac{\hat{\rho} \hat{c}^2}{|J|} \nabla_{\xi} \cdot (|J| F^{-1} \hat{U}) = 0. \end{array} \right. \qquad (33) $$

$$ \hat{\rho} \frac{\partial \hat{U}}{\partial t} + F^{-T} \nabla_{\xi} \hat{p} = \hat{\rho} g, \qquad (34) $$

$$ \frac{\partial \hat{p}}{\partial t} + \frac{\hat{\rho} \hat{c}^2}{|J|} \nabla_{\xi} \cdot (|J| F^{-1} \hat{U}) = 0. \qquad (35) $$
The boundary conditions become

$$ \begin{cases} \hat{U} \cdot \hat{n}_b = \hat{u}_b & \text{on } \hat{\Gamma}_b, \\ \hat{p} = p^a & \text{on } \hat{\Gamma}_s, \end{cases} $$

where $\hat{n}_b$ is a unit vector normal to $\hat{\Gamma}_b$ and pointing toward the exterior of the domain. The variables $\hat{\rho}, \hat{p}, \hat{T}$ satisfy the same equation of state

$$ \hat{p} = f_p(\hat{\rho}, \hat{T}), \qquad (36) $$

and the speed of sound is a function of the new variables, $\hat{c} = c(\hat{\rho}, \hat{s})$.

## 2.3 Linearization and wave equation

We assume that the source of the tsunami is a displacement of magnitude $a$ at the seafloor occuring in an ocean at rest as described in Sec. 2.1.3. In particular, for this rest state there is no mean current and the temperature, pressure and density vary only vertically. The magnitude of the displacement is assumed small compared to the water height $H$. The ratio of the bottom displacement amplitude to the water height is denoted $\epsilon = a/H \ll 1$, and the source term can be expressed as

$$ \hat{u}_b = \epsilon \hat{u}_{b,1} + \mathcal{O}(\epsilon^2). \qquad (37) $$

The linearization of the equations (33)-(35) around the rest state corresponds to the following asymptotic expansion,

$$ d(\xi, t) = \epsilon d_1(\xi, t) + \mathcal{O}(\epsilon^2), \qquad (38) $$

$$ \hat{\rho}(\xi, t) = \hat{\rho}_0(\xi) + \epsilon \hat{\rho}_1(\xi, t) + \mathcal{O}(\epsilon^2), \qquad (39) $$

$$ \hat{p}(\xi, t) = \hat{p}_0(\xi) + \epsilon \hat{p}_1(\xi, t) + \mathcal{O}(\epsilon^2). \qquad (40) $$

Note that the displacement has no zero order term, because the reference configuration used to define the Lagrangian description is the state given by the initial conditions. It holds then $d_0 = 0$, $\hat{U}_0 = 0$ and $\hat{\Omega} = \Omega(0)$.

*Remark:* In comparison with the linearization done in Auclair et al. (2021), where the expansion of the density and the pressure is justified with the decomposition into hydrostatic and non-hydrostatic components, the asymptotic expansion (39) - (40) is obtained in a more straightforward way. Indeed it only requires the assumption of a small perturbation.

From the expansion one deduces the following Taylor expansions for the other functions,

$$ \hat{U} = \epsilon \hat{U}_1 + \mathcal{O}(\epsilon^2), \qquad (41) $$

$$ F = \mathbb{I} + \epsilon \nabla_{\xi} d_1 + \mathcal{O}(\epsilon^2), \qquad (42) $$

$$ (F)^{-1} = \mathbb{I} - \epsilon \nabla_{\xi} d_1 + \mathcal{O}(\epsilon^2), \qquad (43) $$

$$ J = 1 + \epsilon \nabla_{\xi} \cdot d_1 + \mathcal{O}(\epsilon^2). \qquad (44) $$
Injecting these expressions in Eq.(33)-(35) yields the system

$$ \frac{\partial}{\partial t}(\hat{\rho}_0 + \varepsilon\hat{\rho}_1) + \varepsilon\hat{\rho}_0\nabla_{\xi}\cdot\hat{U}_1 = \mathcal{O}(\varepsilon^2), \quad (45) $$

$$ \varepsilon\hat{\rho}_0 \frac{\partial \hat{U}_1}{\partial t} + (I - \varepsilon\nabla_{\xi}d_1)^T \nabla_{\xi}\hat{p}_0 + \varepsilon\nabla_{\xi}\hat{p}_1 = (\hat{\rho}_0 + \varepsilon\hat{\rho}_1)g + \mathcal{O}(\varepsilon^2), \quad (46) $$

$$ \frac{\partial}{\partial t}(\hat{\rho}_0 + \varepsilon\hat{p}_1) + \varepsilon\hat{\rho}_0 c^2(\hat{p}_0, \hat{T}_0)\nabla_{\xi}\cdot\hat{U}_1 = \mathcal{O}(\varepsilon^2). \quad (47) $$

By separating the powers of $\varepsilon$ we obtain two systems: a limit system when $\varepsilon \to 0$ and a system for the first order corrections. Since the limit system correspond to the initial conditions described in Sec. 2.1.3, the model reduces to the first order system.

## First-order system: a wave-like equation for the velocity

The system for the correction terms reads in $\hat{\Omega}$,

$$ \hat{\rho}_0 \frac{\partial \hat{U}_1}{\partial t} + \nabla_{\xi} \hat{p}_1 - (\nabla_{\xi} d_1)^T \nabla_{\xi} \hat{p}_0 = \hat{\rho}_1 g, \quad (48) $$

$$ \frac{\partial \hat{\rho}_1}{\partial t} + \hat{\rho}_0 \nabla_{\xi} \cdot \hat{U}_1 = 0, \quad (49) $$

$$ \frac{\partial \hat{p}_1}{\partial t} + \hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1 = 0, \quad (50) $$

with the boundary conditions

$$ \hat{U}_1 \cdot \hat{n}_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \quad (51) $$

$$ \hat{p}_1 = 0 \quad \text{on } \hat{\Gamma}_s. \quad (52) $$

In this system the speed of sound is evaluated at the limit – or background – pressure and temperature, $\hat{c}_0 = c(\hat{\rho}_0, \hat{T}_0)$. In particular, $\hat{c}_0$ can be written as a function of depth. With an adapted temperature profile it is then possible to recover the typical speed of sound profile creating the SOFAR channel.

The pressure $\hat{p}_1$ and density $\hat{\rho}_1$ can be eliminated in (48) thanks to the other equations: differentiating in time (48) and replacing $\hat{\rho}_1$ and $\hat{p}_1$ with (49), (50) we obtain a second order equation for $\hat{U}_1$,

$$ \hat{\rho}_0 \frac{\partial^2 \hat{U}_1}{\partial t^2} - \nabla_{\xi} (\hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1) - (\nabla_{\xi} \hat{U}_1)^T \hat{\rho}_0 g + \hat{\rho}_0 \nabla_{\xi} \cdot \hat{U}_1 g = 0 \quad \text{in } \hat{\Omega}. \quad (53) $$

Using (50), the surface boundary condition (52) is formulated for $\hat{U}_1$, hence the two boundary conditions for the wave-like equation (53) are

$$ \hat{U}_1 \cdot \hat{n}_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \quad (54) $$

$$ \nabla_{\xi} \cdot \hat{U}_1 = 0 \quad \text{on } \hat{\Gamma}_s. \quad (55) $$
The wave-like equation (53) is completed with vanishing initial condition for $\hat{U}_1(0)$ and $\partial_t \hat{U}_1(0)$. The system (53) includes both gravity and acoustic terms. This equation, which describes the velocity of a compressible, non-viscous fluid, in Lagrangian description, is called the Galbrun equation. It is used in helio-seismology and in aeroacoustics (Legendre, 2003; Maeder et al., 2020; Hägg & Berggren, 2021). However, to our knowledge, this equation has never been used to describe the propagation of hydro-acoustic waves.

We show now that the system (53), (54), (55) is energy-preserving. A model describing a physical system should either preserve or dissipate an energy, and this property makes it also possible to write a stable numerical scheme. Here the energy equation is obtained by taking the scalar product of (53) with $\partial_t \hat{U}_1$ and integrating over the domain. After some computations (see Appendix, Sec. 6), we have

$$ \frac{d}{dt} \mathcal{E} = \int_{\hat{\Omega}_b} \hat{\rho}_0 \left( c_0^2 \nabla_{\xi} \cdot \hat{U}_1 - \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \right) \frac{\partial \hat{u}_{b,1}}{\partial t} d\sigma, \quad (56) $$

with the energy being the quadratic functional given by

$$ \mathcal{E} = \int_{\hat{\Omega}} \rho_0 \frac{1}{2} \left| \frac{\partial \hat{U}_1}{\partial t} \right|^2 d\xi + \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 \left( \hat{c}_0 \nabla_{\xi} \cdot \hat{U}_1 - \frac{g}{\hat{c}_0} \hat{U}_1 \cdot e_3 \right)^2 d\xi \\ + \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 N_b (\hat{U}_1 \cdot e_3)^2 d\xi + \frac{1}{2} \int_{\Gamma_s} \hat{\rho}_0 g (\hat{U}_1 \cdot e_3)^2 d\sigma. \quad (57) $$

The scalar $N_b$ is the squared Brunt-Väsisälä frequency, defined by

$$ N_b(\xi^3) = - \left( \frac{g^2}{\hat{c}_0(\xi^3)^2} + g \frac{\hat{\rho}'_0(\xi^3)}{\hat{\rho}_0(\xi^3)} \right). \quad (58) $$

The Brunt-Väsisälä frequency, or buoyancy frequency, is closely related to the internal waves that appear in a stratified medium (Gill, 1982, Chap. 6). In the ocean, the usual values of $N_b$ are around $10^{-8} \text{ rad}^2\text{s}^{-2}$ (King et al., 2012).

The physical interpretation for the different terms in the energy is clearer when one writes the wave-like equation (53) in terms of the displacement $d_1$ instead of the velocity $\hat{U}_1$. Using $\partial_t d_1 = \hat{U}_1$ and integrating Eq. (53) once in time with the vanishing initial conditions for the displacement, one obtains

$$ \hat{\rho}_0 \frac{\partial^2 d_1}{\partial t^2} - \nabla_{\xi} (\hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot d_1) - (\nabla_{\xi} d_1)^T \hat{\rho}_0 g + \hat{\rho}_0 \nabla_{\xi} \cdot d_1 g = 0 \quad \text{in } \hat{\Omega}. \quad (59) $$

The exact same steps of Appendix 6, with $d_1$ instead of $\hat{U}_1$, yield the energy equation

$$ \frac{d}{dt} \mathcal{E}_d = \int_{\hat{\Omega}_b} \hat{\rho}_0 (c_0^2 \nabla_{\xi} \cdot d_1 - \hat{\rho}_0 g d_1 \cdot e_3) \hat{u}_{b,1} d\sigma, \quad (60) $$
with the energy

$$
\begin{equation}
\begin{split}
\mathcal{E}_d = \int_{\hat{\Omega}} \hat{\rho}_0 \frac{1}{2} |\hat{U}_1|^2 d\xi &+ \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 \left( \hat{c}_0 \nabla_{\xi} \cdot \mathbf{d}_1 - \frac{g}{\hat{c}_0} \mathbf{d}_1 \cdot \mathbf{e}_3 \right)^2 d\xi \\
&+ \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 N_b (\mathbf{d}_1 \cdot \mathbf{e}_3)^2 d\xi + \frac{1}{2} \int_{\Gamma_s} \hat{\rho}_0 g (\mathbf{d}_1 \cdot \mathbf{e}_3)^2 d\sigma.
\end{split}
\tag{61}
\end{equation}
$$

The first term in Eq. (61) is the kinetic energy. We show that the second term of Eq. (61) corresponds to the acoustic energy. First, using Eq. (50) with the vanishing initial conditions yields $\hat{\rho}_0 \hat{c}_0^2 \cdot \nabla_{\xi} \mathbf{d}_1 = -\hat{p}_1$. We define then the acoustic pressure

$$ p_a = \hat{p}_1 - \nabla \hat{p}_0 \cdot \mathbf{d}_1. \tag{62} $$

Indeed, in Lagrangian coordinates, the pressure perturbation $\hat{p}_1$ has two contributions: the small variations in acoustic pressure, and the background pressure being evaluated at a new position. With the definition of $p_a$ and the equation (24), it holds for the second term of (61)

$$ \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 \left( \hat{c}_0 \nabla_{\xi} \cdot \mathbf{d}_1 - \frac{g}{\hat{c}_0} \mathbf{d}_1 \cdot \mathbf{e}_3 \right)^2 d\xi = \frac{1}{2} \int_{\hat{\Omega}} \frac{p_a^2}{\hat{\rho}_0 \hat{c}_0^2} d\xi, \tag{63} $$

which is the usual expression for the acoustic energy (Lighthill, 1978). The last term of (61) is the potential energy associated to the surface waves. Finally, the third term of (61) is the potential energy associated to the internal gravity waves (Lighthill, 1978), under the condition

$$ N_b > 0. \tag{64} $$

When $N_b$ is positive, it is denoted $N_b = N^2$, where $N$ is the buoyancy frequency. The sign of $N_b$ depends on the choice of the state at equilibrium: $\hat{\rho}'_0 = d\hat{\rho}_0/dz$ has to be negative and satisfy

$$ \frac{|\hat{\rho}'_0|}{\hat{\rho}_0} > \frac{g}{\hat{c}_0^2}. \tag{65} $$

With the term in $g^2/\hat{c}_0^2$, we see that the compressibility tends to take the fluid away from its equilibrium. The stratification of the fluid must be strong enough to counter this effect and keep the system stable (see the discussion in Gill, 1982, Chap. 3). As a consequence, if one wants the model to preserve the energy of the system, the background density should not be assumed homogeneous. In the following, we assume that the fluid has a stable stratification, namely that the function $N_b$ is assumed always positive. We will use the notation $N^2$ in the rest of this paper.

**Remark:** According to the equation of state (when the salinity is neglected) $\rho = f_{\rho}(p, T)$, the background density varies because of the variations in temperature and in pressure. The temperature profile can be chosen homogeneous, but
the effect of gravity – see Eq. (24) – prevents the pressure to be independent of depth. Hence in a model with gravity, the fluid is always stratified with the density increasing with depth.

Remark: One can notice that the condition (65) is not explicit in Eq. (53). We obtain this condition when imposing that the energy $E_d$ is positive.

# 3 Derivation of simplified models

In order to compare with existing models, we present several simplifications of our model. We first show that in the barotropic case the system (53) - (55) is equivalent to the first-order scalar equation of Longuet-Higgins (1950). Our model also reduces to well-known models in the acoustic and incompressible asymptotic regimes, as demonstrated below. Further numerical implementations of our model will make it possible to quantify the impact of assumptions made in more simple models, in particular in the case of acoustic-gravity wave generation by earthquakes or landslides in the ocean.

## 3.1 The barotropic case

We consider the barotropic case, which is a very common assumption for the study of hydro-acoustic waves (see for example Longuet-Higgins (1950), Stissnie (2010)). For a barotropic fluid the pressure is a function of the density only,

$$f_p(\rho, s) = f_p(\rho) = p.$$

Then, using Eq. (22) and the definition of the speed of sound,

$$\hat{p}'_0 = \hat{p}'_0 \frac{df_p}{d\rho}(\hat{p}_0) \Rightarrow -\hat{p}_0 g = \hat{p}'_0 \hat{c}_0^2, \quad (66)$$

meaning that the Brunt-Väisälä frequency vanishes, $N^2 = 0$. This corresponds to the case where the density is stratified because of the variation of pressure only. In order to use this equality, we divide Eq. (53) by $\hat{p}_0$,

$$\frac{\partial^2 \hat{U}_1}{\partial t^2} - \nabla_{\xi} (\hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1) - \left( \frac{\hat{p}'_0}{\rho_0} \hat{c}_0^2 + g \right) \nabla_{\xi} \cdot \hat{U}_1 e_3 + \nabla_{\xi} (\hat{U}_1 \cdot e_3) g = 0, \quad (67)$$

and when Eq. (66) holds, the equation (67) can be simplified and reads

$$\frac{\partial^2 \hat{U}_1}{\partial t^2} - \nabla_{\xi} (\hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1) + \nabla_{\xi} (\hat{U}_1 \cdot e_3) g = 0. \quad (68)$$

Taking the curl of Eq. (68) yields

$$\frac{\partial^2 \nabla_{\xi} \times \hat{U}_1}{\partial t^2} = 0 \text{ in } \hat{\Omega}. \quad (69)$$
With the vanishing initial conditions we obtain that the velocity of a barotropic fluid is irrotational. This is a well-known result, since the fluid is also inviscid and subject to a potential force only (Guyon, 2001, Chap. 7). By the Helmholtz decomposition theorem (Girault & Raviart, 1986) the fluid velocity is written as the gradient of a potential $\psi$ defined up to a constant. The expression $\hat{U}_1 = \nabla_{\xi}\psi$ is used in Eq. (68), to obtain

$$ \nabla_{\xi} \left( \frac{\partial^2 \psi}{\partial t^2} - \hat{c}_0^2 \Delta_{\xi} \psi + g \frac{\partial \psi}{\partial \xi^3} \right) = 0. \qquad (70) $$

The potential $\psi$ being defined up to a constant, it can always be sought as the solution of

$$ \frac{\partial^2 \psi}{\partial t^2} - \hat{c}_0^2 \Delta_{\xi} \psi + g \frac{\partial \psi}{\partial \xi^3} = 0. \qquad (71) $$

The equation (71) is multiplied by $\hat{\rho}_0/\hat{c}_0^2$, and we use $g/\hat{c}_0^2 = -\hat{\rho}'_0/\hat{\rho}_0$,

$$ \frac{\hat{\rho}_0}{\hat{c}_0^2} \frac{\partial^2 \psi}{\partial t^2} - \hat{\rho}_0 \Delta_{\xi} \psi - \hat{\rho}'_0 \frac{\partial \psi}{\partial \xi^3} = 0. \qquad (72) $$

And since $\hat{\rho}_0$ depends only on $\xi^3$, the two last terms can be rewritten,

$$ \frac{\hat{\rho}_0}{\hat{c}_0^2} \frac{\partial^2 \psi}{\partial t^2} - \nabla_{\xi} \cdot (\hat{\rho}_0 \nabla_{\xi} \psi) = 0. \qquad (73) $$

Hence $\psi$ satisfies a wave equation. The boundary conditions are then deduced from Eq. (54) and Eq. (55),

$$ \left\{ \begin{array}{l} \nabla_{\xi} \psi \cdot \hat{n}_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \\ \\ \hat{c}_0^2 \Delta_{\xi} \psi = \frac{\partial^2 \psi}{\partial t^2} + g \frac{\partial \psi}{\partial \xi^3} = 0 \quad \text{on } \hat{\Gamma}_s. \end{array} \right. \qquad (74) $$

$$ \hat{c}_0^2 \Delta_{\xi} \psi = \frac{\partial^2 \psi}{\partial t^2} + g \frac{\partial \psi}{\partial \xi^3} = 0 \quad \text{on } \hat{\Gamma}_s. \qquad (75) $$

The system (71),(74),(75) is the first-order system obtained in Longuet-Higgins (1950). In Longuet-Higgins (1950), the derivation is quite different since the irrotationality assumption is made independantly from the fact that the fluid is barotropic, and the boundary conditions are obtained from a linearized surface condition. The linearization made in Longuet-Higgins (1950) gives exactly the same result as the linearization strategy we have presented.

We show that the system (73),(74),(75) is energy-preserving. The equation (73) is multiplied by $\partial_t \psi$ and integrated by parts,

$$ \int_{\hat{\Omega}} \hat{\rho}_0 \frac{\partial \psi}{\partial t} \frac{\partial^2 \psi}{\partial t^2} d\xi + \int_{\hat{\Omega}} \hat{\rho}_0 \nabla \left( \frac{\partial \psi}{\partial t} \right) \cdot \nabla \psi d\xi \\ - \int_{\hat{\Gamma}_s} \hat{\rho}_0 \frac{\partial \psi}{\partial t} \nabla \psi \cdot e_3 d\sigma + \int_{\hat{\Gamma}_b} \hat{\rho}_0 \frac{\partial \psi}{\partial t} \nabla \psi \cdot n_b d\sigma = 0. \qquad (76) $$

With the boundary conditions (74) - (75) and after simplifications it holds

$$ \frac{d}{dt} \mathcal{E}_{bar} = - \int_{\Gamma_b} \hat{\rho}_0 \frac{\partial \psi}{\partial t} \hat{u}_{b,1} d\sigma, \qquad (77) $$
where the energy $\mathcal{E}_{bar}$ is defined by

$$ \mathcal{E}_{bar} = \frac{1}{2} \int_{\hat{\Omega}} \frac{\hat{\rho}_0}{\hat{c}_0^2} \left( \frac{\partial \psi}{\partial t} \right)^2 d\xi + \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 |\nabla \psi|^2 d\xi + \frac{1}{2} \int_{\Gamma_s} \frac{\hat{\rho}_0}{g} \left( \frac{\partial \psi}{\partial t} \right)^2 d\xi. \quad (78) $$

The first term of Eq. (78) is the acoustic energy. Indeed, with (14) and (66) one can show that the acoustic pressure $p_a$ and the potential $\psi$ satisfy the usual relation $p_a = -\hat{\rho}_0 \partial_t \psi$ (Lighthill, 1978, Chap.3). The second term of (78) is the kinetic energy. Finally, with (75), one sees that the third term of (78) is the potential energy of the surface waves. In order to obtain the energy equation for the barotropic system (73), it is necessary to use the background density $\hat{\rho}_0$ even if it does not appear in Eq. (73). The correct manipulation for writing the energy equation was found by comparison with the general case described by Eq. (53).

Finally, note that when assuming a homogeneous density in the equation (73), the system (73)-(75) reduce to

$$ \frac{\partial^2 \psi}{\partial t^2} - \hat{c}_0^2 \Delta \psi = 0 \quad \text{in } \hat{\Omega} \quad (79) $$

$$ \nabla_{\xi} \psi \cdot \hat{n}_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \quad (80) $$

$$ \frac{\partial^2 \psi}{\partial t^2} + g \frac{\partial \psi}{\partial \xi^3} = 0 \quad \text{on } \hat{\Gamma}_s. \quad (81) $$

and the energy equation (78) is not modified by this assumption. However, assuming a homogeneous density is not compatible with the derivation of the system (73)-(75), which relies on the equality $g/\hat{c}_0^2 = -\hat{\rho}'_0/\hat{\rho}_0$. The model (79)-(81) can be understood as a barotropic model with the additional assumption that both $-\hat{\rho}'_0/\hat{\rho}_0$ and $g/\hat{c}_0^2$ are neglected inside the domain.

## 3.2 Two asymptotic regimes of the system

In this section we write the limit models for two asymptotic regimes of the system (53)-(55). We consider the incompressible regime, where the acoustic waves are neglected, and the acoustic regime, where the effect of gravity is neglected. The wave equation (53) is written in non-dimensional form, and we show that it depends on a small non-dimensional parameter. A simplified model is then obtained by passing formally to the limit when the small parameter vanishes. By making the appropriate choice for the time scale we obtain first an incompressible approximation, then an acoustic approximation.

### 3.2.1 Non-dimensional equation

We introduce the following characteristic scales for the system: a time $\tau$, a horizontal scale $L$, a vertical scale $H$, a density $\bar{\rho}$, and a fluid velocity $U$. Since the speed of sound is not assumed constant, we denote by $C$ its characteristic magnitude. Finally, the surface waves velocity is of the order of $\sqrt{gH}$ (Constantin,
2009). We focus on a non-shallow water formulation, hence we take $L = H$. For a shallow water version of the equation one would choose $H \ll L$.

Two dimensionless numbers are introduced: the Froude number and the Mach number, respectively defined by

$$ Fr = \frac{U}{\sqrt{gH}}, \quad Ma = \frac{U}{C}. \qquad (82) $$

To fix the idea, we choose the following numerical values respectively for the speed of sound, the fluid velocity and the surface waves velocity: $C \sim 1480 \text{ m s}^{-1}$, $U \sim 1 \text{ m s}^{-1}$ and $\sqrt{gH} \sim 100 \text{ m s}^{-1}$. The dimensionless numbers are then

$$ Fr = 0.01, \quad Ma = 6.10^{-4}. \qquad (83) $$

The characteristic scale for time will be fixed later, as it will depend on the regime we want to study. The variables are put in non-dimensional form and the dimensionless variables are denoted with a $\tilde{\cdot}$, except for the space and time variable for the sake of conciseness. The adimensionned domain is denoted by $\tilde{\Omega}$ and its surface and bottom boundary are respectively $\tilde{\Gamma}_s$ and $\tilde{\Gamma}_b$. The non-dimensional system reads, after simplification by the factor $\tilde{\rho}U$,

$$ \frac{\tilde{\rho}_0}{\tau^2} \frac{\partial^2 \tilde{U}_1}{\partial t^2} - \frac{C^2}{L^2} \nabla_{\tilde{\xi}} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\tilde{\xi}} \cdot \tilde{U}_1) + \frac{g}{L} \tilde{\rho}_0 (\nabla_{\tilde{\xi}} (\tilde{U}_1 \cdot e_3) - \nabla_{\tilde{\xi}} \cdot \tilde{U}_1 e_3) = 0, \quad (84) $$

with the boundary conditions

$$ \tilde{U}_1 \cdot \tilde{n}_b = \tilde{u}_{b,1} \quad \text{on } \tilde{\Gamma}_b, \qquad (85) $$

$$ \nabla_{\tilde{\xi}} \cdot \tilde{U}_1 = 0 \quad \text{on } \tilde{\Gamma}_s, \qquad (86) $$

where $\tilde{u}_{b,1}$ is a dimensionless source term.

### 3.2.2 Incompressible limit

We show that in the incompressible regime, our model is an extension of the classical free-surface Poisson equation to the case of a variable background density.

In order to study the incompressible limit, the characteristic time $\tau$ is chosen to follow the surface waves, which are much slower than the acoustic waves. We take $L/\tau = \sqrt{gH}$. The equation (84) becomes

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_1}{\partial t^2} - \frac{\text{Fr}}{\text{Ma}} \nabla_{\tilde{\xi}} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\tilde{\xi}} \cdot \tilde{U}_1) + \tilde{\rho}_0 (\nabla_{\tilde{\xi}} (\tilde{U}_1 \cdot e_3) - \nabla_{\tilde{\xi}} \cdot \tilde{U}_1 e_3) = 0. \quad (87) $$

The small parameter $\delta = Ma/\text{Fr} \sim 6.10^{-2}$ is introduced in the equation,

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_1}{\partial t^2} - \frac{1}{\delta^2} \nabla_{\tilde{\xi}} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\tilde{\xi}} \cdot \tilde{U}_1) + \tilde{\rho}_0 (\nabla_{\tilde{\xi}} (\tilde{U}_1 \cdot e_3) - \nabla_{\tilde{\xi}} \cdot \tilde{U}_1 e_3) = 0, \quad (88) $$
and the goal is now to calculate the limit of equation (88) when $\delta$ goes to zero. We make the following ansatz for $\tilde{U}_1$,

$$ \tilde{U}_1 = \tilde{U}_{1,0} + \delta^2 \tilde{U}_{1,2} + \mathcal{O}(\delta^3), \qquad (89) $$

where $\tilde{U}_{1,0}$, $\tilde{U}_{1,1}$ and $\tilde{U}_{1,2}$ are independent of $\delta$. Since Eq. (88) has only even powers of $\delta$, the term $\tilde{U}_{1,1}$ is equal to zero. Replacing $\tilde{U}_1$ by its ansatz in the wave equation (88) and separating the powers of $\delta$ yields an equation for each term of the asymptotic development of $\tilde{U}_1$. The equation obtained with the terms in $\delta^{-2}$ reads

$$ \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,0}) = 0, \qquad (90) $$

and the equation obtained with the terms $\delta^0$ reads

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,2}) + \tilde{\rho}_0 (\nabla_{\xi} (\tilde{U}_{1,0} \cdot e_3) - \nabla_{\xi} \cdot \tilde{U}_{1,0} e_3) = 0. \quad (91) $$

With the terms in $\delta^0$ of the boundary conditions we have

$$ \left\{ \begin{array}{ll} \nabla_{\xi} \cdot \tilde{U}_{1,0} = 0 & \text{on } \tilde{\Gamma}_s, \\ \tilde{U}_{1,0} \cdot \tilde{n}_b = \tilde{u}_{b,1} & \text{on } \tilde{\Gamma}_b. \end{array} \right. \qquad (92) $$

$$ \tilde{U}_{1,0} \cdot \tilde{n}_b = \tilde{u}_{b,1} \quad \text{on } \tilde{\Gamma}_b. \qquad (93) $$

And the terms in $\delta^2$ of the boundary conditions read

$$ \left\{ \begin{array}{ll} \nabla_{\xi} \cdot \tilde{U}_{1,2} = 0 & \text{on } \tilde{\Gamma}_s, \\ \tilde{U}_{1,2} \cdot \tilde{n}_b = 0 & \text{on } \tilde{\Gamma}_b. \end{array} \right. \qquad (94) $$

$$ \tilde{U}_{1,2} \cdot \tilde{n}_b = 0 \quad \text{on } \tilde{\Gamma}_b. \qquad (95) $$

We show now that the limit model represents an incompressible flow. The Helmholtz decomposition of $\tilde{U}_{1,0}$ reads

$$ \tilde{U}_{1,0} = \nabla_{\xi} \varphi_{1,0} + \nabla_{\xi} \times \psi_{1,0}, \qquad (96) $$

where $\varphi_{1,0}$ vanishes on $\tilde{\Gamma}_s$ and $\tilde{\Gamma}_b$. Injecting the decomposition of $\tilde{U}_{1,0}$ in Eq. (90) yields

$$ \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \Delta_{\xi} \varphi_{1,0}) = 0, \qquad (97) $$

hence the term inside the gradient is constant in space. Since the velocity $\tilde{U}_{1,0}$ is equal to zero at infinity, we obtain that $\Delta_{\xi} \varphi_{1,0} = 0$ in $\tilde{\Omega}$ (the quantity $\tilde{\rho}_0 \tilde{c}_0$ being always strictly positive). With the vanishing boundary conditions for $\varphi_{1,0}$, we obtain that $\varphi_{1,0}$ is equal to zero everywhere in $\tilde{\Omega}$. Then, taking the divergence of $\tilde{U}_{1,0}$ yields

$$ \nabla_{\xi} \cdot \tilde{U}_{1,0} = \nabla_{\xi} \cdot (\nabla_{\xi} \times \psi_{1,0}) = 0, \qquad (98) $$

hence $\tilde{U}_{1,0}$ is divergence-free.

Now, using the property $\nabla \cdot \tilde{U}_{1,0}$ in the equation (91) and rearranging some terms, we obtain

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,2}) + \nabla_{\xi} (\tilde{\rho}_0 \tilde{U}_{1,0} \cdot e_3) - \tilde{\rho}'_0 (\tilde{U}_{1,0} \cdot e_3) e_3 = 0. \quad (99) $$
Taking the curl of this equation yields

$$ \nabla_{\xi} \times \left( \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \tilde{\rho}'_0 (\tilde{U}_{1,0} \cdot e_3) e_3 \right) = 0, \quad (100) $$

This means that these terms can be expressed as the gradient of a potential function defined up to a constant and denoted $-\tilde{\varphi}_0$,

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \tilde{\rho}'_0 (\tilde{U}_{1,0} \cdot e_3) e_3 = -\nabla_{\xi} \tilde{\varphi}_0. \quad (101) $$

The new function $\tilde{\varphi}_0$ can be understood as the Lagrange multiplier for the incompressibility constraint. However, one must be cautious that $\tilde{\varphi}_0$ is not similar to a pressure in this case, and rather plays the role of a velocity potential, as we will see later in the case of homogeneous density. The function $\tilde{\varphi}_0$ can be expressed differently. By using the definition (101) in the equation (99) we have

$$ \nabla_{\xi} \left( -\tilde{\varphi}_0 - \tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,2} + \tilde{\rho}_0 \tilde{U}_{1,0} \cdot e_3 \right) = 0, \quad (102) $$

and since the potential $\tilde{\varphi}_0$ is defined up to a constant, it can be chosen such that, in $\hat{\Omega}$, we have

$$ \tilde{\varphi}_0 = -\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,2} + \tilde{\rho}_0 \tilde{U}_{1,0} \cdot e_3. \quad (103) $$

We deduce from this equality and (94) the boundary condition

$$ \tilde{\varphi}_0 = \tilde{\rho}_0 \tilde{U}_{1,0} \cdot e_3 \quad \text{on } \hat{\Gamma}_s. \quad (104) $$

To recover a dimensional system, the terms are multiplied by their corresponding characteristic scales, and $\hat{\varphi}_0 = \bar{\rho} U \tilde{\varphi}_0$ is defined. The limit solution $\hat{U}_{1,0} = U \tilde{U}_{1,0}$ satisfies

$$ \hat{\rho}_0 \frac{\partial^2 \hat{U}_{1,0}}{\partial t^2} - g \hat{\rho}'_0 (\hat{U}_{1,0} \cdot e_3) e_3 + g \nabla_{\xi} \hat{\varphi}_0 = 0 \quad \text{in } \hat{\Omega}, \quad (105) $$

$$ \nabla_{\xi} \cdot \hat{U}_{1,0} = 0 \quad \text{in } \hat{\Omega}, \quad (106) $$

with the boundary conditions

$$ \hat{U}_{1,0} \cdot \hat{n}_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \quad (107) $$

$$ \nabla_{\xi} \cdot \hat{U}_{1,0} = 0 \quad \text{on } \hat{\Gamma}_s, \quad (108) $$

$$ \hat{\varphi}_0 = \hat{\rho}_0 \hat{U}_{1,0} \cdot e_3 \quad \text{on } \hat{\Gamma}_s, \quad (109) $$

We show that the model (105)-(109) preserves an energy. Taking the scalar product of (105) with $\partial_t \hat{U}_{1,0}$ and integrating over $\hat{\Omega}$ yields

$$ \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 \left| \frac{\partial \hat{U}_{1,0}}{\partial t} \right|^2 d\xi - \int_{\hat{\Omega}} g \hat{\rho}'_0 (\hat{U}_{1,0} \cdot e_3) e_3 \cdot \frac{\partial \hat{U}_{1,0}}{\partial t} d\xi + \int_{\hat{\Omega}} g \frac{\partial \hat{U}_{1,0}}{\partial t} \cdot \nabla_{\xi} \hat{\varphi}_0 d\xi = 0. \quad (110) $$
The last term of Eq. (110) is integrated by parts. With the vanishing divergence of $\hat{U}_{1,0}$ and the bottom condition (107) it holds

$$ \int_{\hat{\Omega}} g \frac{\partial \hat{U}_{1,0}}{\partial t} \cdot \nabla_{\xi} \hat{\varphi}_0 \, d\xi = \int_{\hat{\Gamma}_s} g \hat{\varphi}_0 \frac{\partial \hat{U}_{1,0}}{\partial t} \cdot e_3 \, d\sigma - \int_{\hat{\Gamma}_b} g \hat{\varphi}_0 \frac{\partial \hat{u}_{b,1}}{\partial t} \, d\sigma, \quad (111) $$

then $\hat{\varphi}_0$ is replaced in the surface integral using Eq. (109),

$$ \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 \left| \frac{\partial \hat{U}_{1,0}}{\partial t} \right|^2 d\xi - \int_{\hat{\Omega}} g \hat{\rho}'_0 (\hat{U}_{1,0} \cdot e_3) e_3 \cdot \frac{\partial \hat{U}_{1,0}}{\partial t} d\xi \\ + \int_{\hat{\Gamma}_s} g \hat{\rho}_0 \hat{U}_{1,0} \cdot e_3 \frac{\partial \hat{U}_{1,0}}{\partial t} \cdot e_3 \, d\sigma = \int_{\hat{\Gamma}_b} g \hat{\varphi}_0 \frac{\partial \hat{u}_{b,1}}{\partial t} \, d\sigma. \quad (112) $$

By defining the energy

$$ \mathcal{E}_{incomp} = \frac{1}{2} \int_{\hat{\Omega}} \hat{\rho}_0 \left| \frac{\partial \hat{U}_{1,0}}{\partial t} \right|^2 d\xi - \frac{1}{2} \int_{\hat{\Omega}} g \hat{\rho}'_0 |\hat{U}_{1,0} \cdot e_3|^2 d\xi + \frac{1}{2} \int_{\hat{\Gamma}_s} g \hat{\rho}_0 |\hat{U}_{1,0} \cdot e_3|^2 d\sigma, \quad (113) $$

the equation (110) can be formulated in the following way,

$$ \frac{d}{dt} \mathcal{E}_{incomp} = \int_{\Gamma_b} g \varphi_0 \frac{\partial \hat{u}_{b,1}}{\partial t}. \quad (114) $$

Each term of $\mathcal{E}_{incomp}$ has the same interpretation as in $\mathcal{E}$. Note that the acoustic term of $\mathcal{E}$ is not present in $\mathcal{E}_{incomp}$. The potential energy associated to the internal waves is also written differently, as in the formal limit $\hat{c}_0 \to \infty$ the buoyancy frequency reads $N^2 = -g\hat{\rho}'_0/\hat{\rho}_0$.

**Remark:** The condition $|\hat{\rho}'_0|/\hat{\rho}_0 > g/\hat{c}_0^2$ is no longer required because the destabilizing effects in the energy equation (56) come from the compressibility, and here it is neglected. This can be seen by formally assuming that the sound speed is infinite, then the squared buoyancy frequency reads $N^2 = -g\hat{\rho}'_0/\hat{\rho}_0$. Density must still decrease with depth, but can be homogeneous.

The system (106)-(105) represents an incompressible fluid. However, this system is different from the classical Poisson equation found in the literature (Lighthill, 1978) because of the assumption of a nonhomogeneous background density. For the sake of comparison with other models, assume now that the ocean at rest has a homogeneous density, $\hat{\rho}'_0 = 0$. Taking the divergence of Eq. (105) yields

$$ \Delta_{\xi} \hat{\varphi}_0 = 0. \quad (115) $$

The boundary conditions are written differently to ease the comparison. The bottom boundary condition is obtained by taking the scalar product of Eq. (105) with $\hat{n}_b$, and replacing the first term with Eq. (107) differentiated twice in time,

$$ - \hat{\rho}_0 \frac{\partial^2 \hat{u}_{b,1}}{\partial t^2} - g \hat{\rho}'_0 (\hat{U}_{1,0} \cdot e_3) e_3 \cdot \hat{n}_b + g \nabla_{\xi} \hat{\varphi}_0 \cdot \hat{n}_b = 0. \quad (116) $$
For the surface condition, the equation (109) is differentiated twice in time and the term in $\partial_{tt}^2 \hat{U}_{1,0}$ is replaced with (105),

$$ \frac{\partial^2 \hat{\varphi}_0}{\partial t^2} - g \rho_0' (\hat{U}_{1,0} \cdot e_3) + g \frac{\partial \hat{\varphi}_0}{\partial \xi^3} = 0 \quad \text{on } \hat{\Gamma}_s. \qquad (117) $$

With the assumption of a homogeneous density, the boundary conditions (116), (117) read then

$$ \left\{ \begin{array}{ll} \nabla_{\xi} \hat{\varphi}_0 \cdot \hat{n}_b = -\hat{\rho}_0 g \hat{u}_{b,1} & \text{on } \hat{\Gamma}_b, \\ \\ \frac{\partial^2 \hat{\varphi}_0}{\partial t^2} + g \frac{\partial \hat{\varphi}_0}{\partial \xi^3} = 0 & \text{on } \hat{\Gamma}_s. \end{array} \right. \qquad (118) $$

$$ \frac{\partial^2 \hat{\varphi}_0}{\partial t^2} + g \frac{\partial \hat{\varphi}_0}{\partial \xi^3} = 0 \quad \text{on } \hat{\Gamma}_s. \qquad (119) $$

The Poisson equation (115) with boundary conditions (118) - (119) is the system satisfied by the velocity flow of an incompressible homogeneous free-surface fluid (Lighthill, 1978, Chap. 3.1). Note that it was required that $\rho_0' \neq 0$ in the system (53) to obtain a priori positive energy. Here this assumption is dropped, however a rather simple expression for the preserved energy can be derived: multiplying Eq. (115) by $\partial_t \hat{\varphi}_0$, integrating by parts and using Eq. (118) - (119), we obtain

$$ \int_{\hat{\Omega}} \Delta_{\xi} \hat{\varphi}_0 \frac{\partial \hat{\varphi}_0}{\partial t} d\xi = - \int_{\hat{\Omega}} \nabla_{\xi} \varphi \cdot \nabla_{\xi} \left( \frac{\partial \hat{\varphi}_0}{\partial t} \right) d\xi - \int_{\Gamma_s} \frac{1}{g} \frac{\partial \hat{\varphi}_0}{\partial t} \frac{\partial^2 \hat{\varphi}_0}{\partial t^2} d\sigma + \int_{\Gamma_b} \hat{\rho}_0 g \frac{\partial \hat{\varphi}_0}{\partial t} \hat{u}_{b,1} d\sigma. \qquad (120) $$

We define the energy

$$ \mathcal{E}_{\text{Poisson}} = \frac{1}{2} \left( \int_{\hat{\Omega}} |\nabla_{\xi} \hat{\varphi}|^2 d\xi + \int_{\hat{\Gamma}_s} \frac{1}{g} \left( \frac{\partial \hat{\varphi}_0}{\partial t} \right)^2 d\sigma \right), \qquad (121) $$

Then it holds

$$ \frac{d}{dt} \mathcal{E}_{\text{Poisson}} = - \int_{\Gamma_b} \hat{\rho}_0 g \frac{\partial \hat{\varphi}_0}{\partial t} \hat{u}_{b,1} d\sigma. \qquad (122) $$

By comparison with the energy of the barotropic system (78), we see that the first term of (121) is the kinetic energy and the second term of (121) is the potential energy associated to the surface waves.

### 3.2.3 Acoustic limit

Another possible simplification of the system (53) - (55) is to keep only the acoustic terms. This choice is justified for short time scale, because the propagation speed of the acoustic waves and the gravity waves have different orders of magnitude (Longuet-Higgins, 1950). Here we show that in the acoustic limit, the model reduces to a classical acoustic equation.
With the timescale $L/\tau = C$, corresponding to the acoustic wave, and with the same small parameter $\delta = Ma/Fr$ as before, the system (84) becomes

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_1}{\partial t^2} - \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_1) + \delta^2 \tilde{\rho}_0 (\nabla_{\xi} (\tilde{U}_1 \cdot e_3) - \nabla_{\xi} \cdot \tilde{U}_1 e_3) = 0 \quad \text{in } \tilde{\Omega}, \quad (123) $$

With the boundary conditions

$$ \left\{ \begin{array}{ll} \tilde{U}_1 \cdot \tilde{n}_b = \tilde{u}_{b,1} & \text{on } \tilde{\Gamma}_b, \\ \tilde{U}_1 \cdot \tilde{n}_b = 0 & \text{on } \tilde{\Gamma}_s. \end{array} \right. \quad (124) $$

$$ \tilde{U}_1 \cdot \tilde{n}_b = 0 \quad \text{on } \tilde{\Gamma}_s. \quad (125) $$

As before, we make the following ansatz for $\tilde{U}_1$,

$$ \tilde{U}_1 = \tilde{U}_{1,0} + \delta^2 \tilde{U}_{1,2} + \mathcal{O}(\delta^3). \quad (126) $$

One can see that the limit term $\delta \to 0$ for the volumic equation (123) is

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \nabla_{\xi} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,0}) = 0. \quad (127) $$

Taking the curl of this equation yields

$$ \frac{\partial^2}{\partial t^2} \left( \nabla_{\xi} \times (\tilde{\rho}_0 \tilde{U}_{1,0}) \right) = 0, \quad (128) $$

hence the curl of $\tilde{\rho}_0 \tilde{U}_{1,0}$ is affine in time. Moreover, it is equal to zero due to the vanishing initial conditions. By the Helmholtz decomposition theorem, the term $\tilde{\rho}_0 \tilde{U}_{1,0}$ can be expressed as the gradient of some function $\tilde{\psi}_0$ defined up to a constant,

$$ \tilde{\rho}_0 \tilde{U}_{1,0} = \nabla_{\xi} \tilde{\psi}_0. \quad (129) $$

By substituting in equation (127) we have

$$ \nabla_{\xi} \left( \frac{\partial^2 \tilde{\psi}_0}{\partial t^2} - \tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot (\tilde{\rho}_0^{-1} \nabla_{\xi} \tilde{\psi}_0) \right) = 0, \quad (130) $$

then it holds

$$ \frac{\partial^2 \tilde{\psi}_0}{\partial t^2} - \tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot (\tilde{\rho}_0^{-1} \nabla_{\xi} \tilde{\psi}_0) = 0, \quad (131) $$

since $\tilde{\psi}_0$ is defined up to a constant. We need the boundary conditions to conclude. Evaluating (127) at the surface yields

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} - \frac{\partial}{\partial \xi^3} (\tilde{\rho}_0 \tilde{c}_0^2 \nabla_{\xi} \cdot \tilde{U}_{1,0}) e_3 = 0 \quad \text{on } \tilde{\Gamma}_s. \quad (132) $$

Using the surface condition (125) in the equation (132) yields

$$ \tilde{\rho}_0 \frac{\partial^2 \tilde{U}_{1,0}}{\partial t^2} = 0 \quad \text{on } \tilde{\Gamma}_s. \quad (133) $$
With the definition of the potential $\tilde{\psi}_0$ it holds

$$ \frac{\partial^2 \nabla \tilde{\psi}_0}{\partial t^2} = 0 \quad \text{on } \tilde{\Gamma}_s, \qquad (134) $$

hence one has

$$ \frac{\partial^2 \tilde{\psi}_0}{\partial t^2} = C(t) \quad \text{on } \tilde{\Gamma}_s, \qquad (135) $$

where $C$ does not depend on space. Moreover, since $\tilde{\psi}_0$ vanishes at infinity the constant $C$ is equal to zero, hence $\partial_{tt}^2 \tilde{\psi}_0 = 0$ on $\tilde{\Gamma}_s$. With the vanishing initial conditions this implies that $\tilde{\psi}_0 = 0$ on $\tilde{\Gamma}_s$. To recover a dimensional system, the terms are multiplied by their corresponding characteristic scales, and $\tilde{\psi}_0 = \bar{\rho}UL \tilde{\psi}_0$ is defined. The system reads then

$$ \frac{\partial^2 \hat{\psi}_0}{\partial t^2} - \hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot (\hat{\rho}_0^{-1} \nabla_{\xi} \hat{\psi}_0) = 0 \quad \text{in } \hat{\Omega}, \qquad (136) $$

with the boundary conditions

$$ \left\{ \begin{array}{ll} \nabla_{\xi} \hat{\psi}_0 \cdot \hat{n}_b = \hat{u}_{b,1} & \text{on } \hat{\Gamma}_b, \\ \hat{\psi}_0 = 0 & \text{on } \hat{\Gamma}_s. \end{array} \right. \qquad (137) $$

$$ \hat{\psi}_0 = 0 \quad \text{on } \hat{\Gamma}_s. \qquad (138) $$

The system (136)-(138) is the classical wave equation for the potential $\hat{\psi}_0$, with a propagation speed $\hat{c}_0^2$ and a non-homogeneous density.

An energy equation can be obtained by multiplying Eq. (136) by $\partial_t \psi / (\rho_0 \hat{c}_0^2)$ and integrating over the domain. The result reads after an integration by parts

$$ \frac{d}{dt} \mathcal{E}_{\text{acoustic}} = - \int_{\hat{\Gamma}_b} \frac{1}{\hat{\rho}_0} \frac{\partial \hat{\psi}_0}{\partial t} \hat{u}_{b,1} \, d\sigma, \qquad (139) $$

where the acoustic energy is

$$ \mathcal{E}_{\text{acoustic}} = \frac{1}{2} \int_{\hat{\Omega}} \frac{1}{\hat{\rho}_0 \hat{c}_0^2} \left( \frac{\partial \hat{\psi}_0}{\partial t} \right)^2 d\xi + \frac{1}{2} \int_{\hat{\Omega}} \frac{1}{\hat{\rho}_0} |\nabla \hat{\psi}_0|^2 d\xi. \qquad (140) $$

With the same analysis as in the previous cases, one can show that the first term of (140) is the acoustic energy, and the second term is the kinetic energy.

**Remark:** In the sections 3.2.2 and 3.2.3, the equations (115) - (119) and (136)-(138) use the Lagrangian description whereas the equations from the literature use the Eulerian description. In the general case, the use of different coordinate systems would cause two problems. First, when doing the change of coordinates, new terms should appear from the space or time differentiation. Second, the description of the domain is different, and this implies that the boundary conditions are not evaluated at the same location. In the next section we will show that the first problem does not exist in our case, due to the lack of
a background velocity. As for the second problem, the linear eulerian models are obtained by evaluating the boundary conditions at a fixed water height. To this regard they use the same boundary as if they were in a Lagrangian description of the domain, so that the comparison remains valid.

The equations with their boundary conditions and the associated energy, for the general model and its different simplifications, are summarized in Table 1.

# 4 The model in Eulerian coordinates

The equations we have been working on are defined on the reference domain $\hat{\Omega}$. However, the linear equations for the acoustic-gravity waves are generally written in Eulerian coordinates. In order to compare our model with those from the literature, the equations must be formulated on the moving domain $\Omega(t)$. In this section we present a method to write the system in Eulerian coordinate.

## 4.1 General method

The aim is to write the equation on a moving domain $\Omega(t)$, hence a transformation $\phi : \hat{\Omega} \to \Omega(t)$ is needed. We start by using a first order approximation of the real transformation $\phi$. The transformation $\phi$ is developed for small displacements,

$$ \phi(\xi, t) = \mathbb{I} + \varepsilon\phi_1(\xi, t) + \mathcal{O}(\varepsilon^2). \quad (141) $$

Let $\phi_{\varepsilon}(\xi, t) = \mathbb{I} + \varepsilon\phi_1(\xi, t)$ be its first order approximation. $\phi_{\varepsilon}$ is used to define the equivalent domain and its boundary,

$$ \Omega_{\varepsilon}(t) = \phi_{\varepsilon}(\hat{\Omega}), \quad \Gamma_{s,eq} = (\phi_{\varepsilon}(\hat{\Gamma}_s)), \quad \Gamma_{b,eq} = (\phi_{\varepsilon}(\hat{\Gamma}_b)). \quad (142) $$

The coordinates on the equivalent domain are written $\mathbf{x} = (x, y, z)$. For any generic function $\hat{X}(\xi, t)$ defined in $\hat{\Omega}$, a function $X(\mathbf{x}, t)$ is defined in $\Omega_{\varepsilon}$ by the following change of variables

$$ X(\mathbf{x}, t) = \hat{X}(\phi_{\varepsilon}^{-1}(\mathbf{x}, t), t), \quad (143) $$

which is equivalent to

$$ \hat{X}(\xi, t) = X(\phi_{\varepsilon}(\xi, t), t), \quad (144) $$

as long as $\phi_{\varepsilon}$ is invertible. Then, if the function $\hat{X}$ has a first-order approximation $\hat{X} = \hat{X}_0 + \varepsilon\hat{X}_1 + \mathcal{O}(\varepsilon^2)$, then the function $X$ also has a first-order approximation $X = X_0 + \varepsilon X_1 + \mathcal{O}(\varepsilon^2)$ and it holds (see Appendix, Sec. 6)

$$ \nabla_{\xi} \hat{X}_0 = \nabla X_0, \quad (145) $$

$$ \frac{\partial \hat{X}_0}{\partial t} = \frac{\partial X_0}{\partial t}, \quad (146) $$

$$ \nabla_{\xi} \hat{X}_1 = (\nabla_{\xi} \mathbf{d}_1)^T \nabla X_0 + \nabla X_1, \quad (147) $$

$$ \frac{\partial \hat{X}_1}{\partial t} = \frac{\partial X_1}{\partial t} + \mathbf{U}_1 \cdot \nabla X_0. \quad (148) $$
In the following, when writing the equations satisfied by the free surface of $\Omega_\epsilon$, we will also use

$$ \frac{\partial \phi_{\epsilon}}{\partial t} = \epsilon \hat{U}_1. \qquad (149) $$

## 4.2 The model in Eulerian coordinates

Using the change of variable (144) in the system (48)-(50) and with the equalities (145)-(148) we obtain the following system for $U_1, p_1, \rho_1$ defined in $\Omega_\epsilon$,

$$ \rho_0 \frac{\partial U_1}{\partial t} + \nabla p_1 = \rho_1 g, \qquad (150) $$

$$ \frac{\partial \rho_1}{\partial t} + \nabla \cdot (\rho_0 U_1) = 0, \qquad (151) $$

$$ \frac{\partial p_1}{\partial t} + \nabla p_0 \cdot U_1 + \rho_0 c_0^2 \nabla \cdot U_1 = 0. \qquad (152) $$

And $p_0, \rho_0$ satisfy the limit equations

$$ \frac{\partial \rho_0}{\partial t} = 0, \qquad (153) $$

$$ \nabla p_0 = \rho_0 g. \qquad (154) $$

In order to close the System (150)-(154), boundary conditions should be prescribed. To get a linear problem one wants to prescribe this condition on the fixed domain $\hat{\Omega}$. To do so we assume in the following that the equations (150)-(154) are defined in $\hat{\Omega}$. It would be true if $\hat{\Omega} \subset \Omega_\epsilon$, but the inclusion is in general not verified. Because of this approximation, errors of order $O(\epsilon)$ may be introduced. For this reason the system in Lagrangian coordinates should be preferred, at least for future extension of this work.

### Boundary conditions and free surface description

Following the approach of Nouguier et al. (2015), we show that a description for the free surface can be obtained. In the following, the components of the fluid velocity are denoted $U_1 = (U_1^1, U_1^2, U_1^3)^T$. The surface is defined by $\Gamma_{s,eq} = \phi_\epsilon(\hat{\Gamma}_s)$, and we assume that at each time $t$ it can be parametrized as the graph $\eta_\epsilon$. The elevation $\eta_\epsilon$ is a function of $x, y$, and $t$ and can be decomposed in the following way,

$$ \eta_{\epsilon}(x, y, t) = H + \epsilon \eta_{1}(x, y, t). \qquad (155) $$

From the correspondence between the free surface and the particle displacement, it holds

$$ \phi_{\epsilon}^{3}(\xi^{1}, \xi^{2}, H, t) = \eta_{\epsilon}(x(\xi^{1}, \xi^{2}, H, t), y(\xi^{1}, \xi^{2}, H, t), t). \qquad (156) $$

Differentiating (156) in time and using the equation (149) yields

$$ \epsilon \hat{U}_1^3(\xi^1, \xi^2, H, t) = \frac{\partial \eta_\epsilon}{\partial t} + \epsilon \hat{U}_1^1(\xi^1, \xi^2, H, t) \frac{\partial \eta_\epsilon}{\partial x} + \epsilon \hat{U}_1^2(\xi^1, \xi^2, H, t) \frac{\partial \eta_\epsilon}{\partial y}. \qquad (157) $$
We use the change of variables $\phi_{\varepsilon}(\xi, t) = \mathbb{I} + \varepsilon\phi_1(\xi, t)$,

$$ \varepsilon U_1^3(\phi_{\varepsilon}(\xi^1, \xi^2, H, t), t) = \frac{\partial \eta_{\varepsilon}}{\partial t} + \varepsilon U_1^1(\phi_{\varepsilon}(\xi^1, \xi^2, H, t), t) \frac{\partial \eta_{\varepsilon}}{\partial x} \\ + \varepsilon U_1^2(\phi_{\varepsilon}(\xi^1, \xi^2, H, t), t) \frac{\partial \eta_{\varepsilon}}{\partial y}. \quad (158) $$

After a Taylor development and keeping only the terms in $\varepsilon$, it holds

$$ U_1^3(x, y, H, t) = \frac{\partial \eta_1}{\partial t}, \quad (159) $$

this is the linearized equation for the free surface. Then the dynamic boundary conditions are linearized. With the change of variables, the boundary conditions (24), (51) and (52) become

$$ U_1 \cdot n_b = u_{b,1} \quad \text{on } \Gamma_{b,eq}, \quad (160) $$

$$ p_0 = p^a \quad \text{on } \Gamma_{s,eq}, \quad (161) $$

$$ p_1 = 0 \quad \text{on } \Gamma_{s,eq}. \quad (162) $$

If we linearize (162) only we would miss the first-order term coming from (161). From (161) and (162) we deduce the boundary condition for the pressure

$$ p_0 + \varepsilon p_1 = p^a \quad \text{on } \Gamma_{s,eq}. \quad (163) $$

A Taylor development of $p_0$ and $p_1$ around $z = H$ on $\Gamma_{s,eq}$ yields

$$ p_0(H) + \varepsilon(p_1(x, y, H, t) + p_0'(H)\eta_1) + \mathcal{O}(\varepsilon^2) = p^a. \quad (164) $$

After an identification of the powers of $\varepsilon$ it holds

$$ p_0(H) = p^a, \quad p_1(x, y, H, t) = \rho_0(x, y, H, t)g \eta_1(x, y, t). \quad (165) $$

In a similar way, the linearization of Eq. (160) reads

$$ U_1^3(x, y, z_b) - U_1^1(x, y, z_b) \partial_x z_b - U_1^2(x, y, z_b) \partial_y z_b = u_{b,1}(x, y, t). \quad (166) $$

Hence the equations for $U_1$, $\rho_1$, $p_1$ can be fully defined on the domain $\hat{\Omega}$, with an error in $\mathcal{O}(\varepsilon^2)$. Finally, note that the system (150)-(152) with the boundary conditions (165), (166) and the kinematic condition (159), is shown to be energy-preserving, locally as well as over a whole water column (Lighthill, 1978; Lotto & Dunham, 2015).

In this section we have derived the linear equation in Eulerian coordinates, even though an approximation on the domain in which the equations are defined was necessary. The computations of Sec. 4.1 also justify that in the absence of mean flow and with the evaluation of the boundary conditions at a fixed height, the linear system in Eulerian coordinates is similar to the one in Lagrangian coordinate, up to terms in $\mathcal{O}(\varepsilon^2)$. At the same time, the linearization in the Lagrangian coordinates is better defined. For this reason the system in Lagrangian coordinates is preferred for the rest of this work. We conclude this paper by the study of the dispersion relation obtained from Eq. (53).
## 5 Dispersion relation

A key aspect of wave models is the related dispersion relation, that we derive here from Eq. (53) and solve numerically. First note that if one defines the equivalent pressure $p_ε$, the equivalent density $ρ_ε$ and the equivalent velocity $U_ε$ by

$$p_{\varepsilon} = p_{0} + \varepsilon p_{1}, \quad \rho_{\varepsilon} = \rho_{0} + \varepsilon \rho_{1}, \quad U_{\varepsilon} = \varepsilon U_{1}, \qquad (167)$$

then a combination of the equations (150)-(154) yields the following system for $p_ε$, $ρ_ε$ and $U_ε$,

$$\rho_0 \frac{\partial U_\varepsilon}{\partial t} + \nabla p_\varepsilon = \rho_\varepsilon g + \mathcal{O}(\varepsilon^2), \qquad (168)$$

$$\frac{\partial \rho_{\varepsilon}}{\partial t} + \nabla \cdot (\rho_0 U_{\varepsilon}) = \mathcal{O}(\varepsilon^2), \qquad (169)$$

$$\frac{\partial p_{\varepsilon}}{\partial t} + \nabla p_0 \cdot U_{\varepsilon} + \rho_0 c_0^2 \nabla \cdot U_{\varepsilon} = \mathcal{O}(\varepsilon^2). \qquad (170)$$

This system is comparable – up to the terms in $\mathcal{O}(\varepsilon^2)$ – to the system studied in the paper from Auclair et al. (2021). In Auclair et al. (2021), the author thoroughly analyze the dispersion relation for the model of a stratified compressible fluid with a constant sound speed.

In order to make the computations clearer, the problem is restricted to a 2-dimensional configuration in $\xi^1$ and $\xi^3$. We also assume that the bottom is flat. Following the approach of Auclair et al. (2021), the wave angular frequency $\omega$ and the horizontal wave-number frequency $k_x$ are defined, and we seek a solution of the form

$$\hat{\rho}_0 \hat{U}_1(\xi^1, \xi^3, t) = \begin{pmatrix} \tilde{U}^1(\xi^3) \\ \tilde{U}^3(\xi^3) \end{pmatrix} e^{i(k_x \xi^1 - \omega t)}. \qquad (171)$$

First Eq. (53) is written differently to make the unknown $\hat{\rho}_0 \hat{U}_1$ appear.

$$\frac{\partial^2 \hat{\rho}_0 \hat{U}_1}{\partial t^2} - \nabla_{\xi} (\hat{c}_0^2 \nabla_{\xi} \cdot (\hat{\rho}_0 \hat{U}_1)) - \nabla_{\xi} \left( \frac{\hat{c}_0^2 N_0^2}{g} \rho_0 \hat{U}_1 \cdot e_3 \right) - g \nabla_{\xi} \cdot (\hat{\rho}_0 \hat{U}_1) e_3 = 0 \quad (172)$$

Injecting the ansatz (171) in Eq. (172) yields

$$\omega^2 \tilde{U}^1 + ik_x \left( \hat{c}_0^2 (ik_x \tilde{U}^1 + (\tilde{U}^3)') + \frac{\hat{c}_0^2 N^2}{g} \tilde{U}^3 \right) = 0, \qquad (173)$$

$$\omega^2 \tilde{U}^3 + \partial_3 \left( \hat{c}_0^2 (ik_x \tilde{U}^1 + (\tilde{U}^3)') + \partial_3 \left( \frac{\hat{c}_0^2 N^2}{g} \tilde{U}^3 \right) + g(ik_x \tilde{U}^1 + (\tilde{U}^3)') \right) = 0. \qquad (174)$$

Using the equation (173), the horizontal component $\tilde{U}^1$ is expressed as a function of the vertical component,

$$\tilde{U}^1 = -ik_x \frac{\hat{c}_0^2 D(\tilde{U}^3)' + (\hat{c}_0^2 - gD)\tilde{U}^3}{D(\omega^2 - \hat{c}_0^2 k_x^2)}, \qquad (175)$$
where D is a depth scale, defined by

$$ \frac{1}{D} = \frac{N^2}{g} + \frac{g}{\hat{c}_0^2} = \frac{\hat{\rho}'_0}{\hat{\rho}_0}. \qquad (176) $$

We also define the quantity

$$ S = 2 \frac{\hat{c}'_0}{\hat{c}_0}. \qquad (177) $$

Replacing $\tilde{U}^1$ in the equation (174) yields, after some computations,

$$ (\tilde{U}^3)'' + \left(\frac{1}{D} + \omega^2 S^2\right) (\tilde{U}^3)' + \left(\frac{\omega^2}{\hat{c}_0^2} + k_x^2 \frac{N^2 - \omega^2}{\omega^2} - \frac{D'}{D^2} + S \left(\frac{g}{\hat{c}_0^2} + \frac{N^2}{g} \frac{\omega^2}{\omega^2 - \hat{c}_0^2 k_x^2}\right)\right) \tilde{U}^3 = 0. \qquad (178) $$

In order to write an harmonic equation the following change of variable is made,

$$ \tilde{U}^3(z) = \tilde{U}^3(H)F(z) \exp \left( \int_z^H \frac{\alpha}{2} dz' \right), \quad \alpha = \frac{1}{D} + \omega^2 S. \qquad (179) $$

Then $F(0) = 0$, $F(H) = 1$ and $F$ satisfies the equation

$$ F'' + k_z^2 F = 0 \qquad (180) $$

where the vertical wave-number $k_z$ is defined by

$$ k_z^2 + k_x^2 \frac{N^2 - \omega^2}{\omega^2} + \frac{\omega^2}{\hat{c}_0^2} - \frac{1 + 2D'}{4D^2} - \frac{1}{2} \omega^2 S' + S \left( \frac{g}{\hat{c}_0^2} + \frac{N^2}{g} \frac{\omega^2}{\omega^2 - \hat{c}_0^2 k_x^2} - \frac{\omega^2}{2D} - \frac{1}{4} \omega^4 S \right) = 0, \qquad (181) $$

The equation (181) is the dispersion relation for the two wave-numbers $k_x$, $k_z$ and the frequency $\omega$. It is a generalization of the inner dispersion relation in Auclair et al. (2021) to the case of a non-constant sound speed. Indeed, with a constant sound speed one has $S = 0$ and (181) is exactly the inner dispersion relation in Auclair et al. (2021).

**Remark:** In the most general case, the scalars $N, D$ and $S$ depend on the depth $z$, hence $k_z$ also depends on $z$. It is then not clear whether the solution to Eq. (180) and the profile $\tilde{U}^3$, can be written explicitely. When $k_z$ does not depend on $z$, as in the study from Auclair et al. (2021), the expression of the profile $\tilde{U}^3$ is used with the boundary conditions to obtain a boundary dispersion relation. In our case $k_z$ is not a constant, and the boundary dispersion relation is not easily deduced.

## Numerical approximation of the dispersion relation

An evaluation of the equation (181) is possible once the limit state for the pressure and the density is computed. The differential equation for the pressure (22) is numerically solved for the temperature profile shown in figure 3a.
<figure><img src="image_5.png" /><figcaption>Figure 3: Temperature, density and sound speed profiles used for the computation of the dispersion relation where $\xi^3 = 0$ is the seafloor and $\xi^3 = 4000 \text{ m}$ is the ocean surface.</figcaption></figure>
Then the density and the speed of sound are computed from the tabulations given in IAPWS-SR7 (2009). Figure 3b, 3c show the obtained density and speed of sound. With these profiles the dispersion relation (181) is computed. Figure 4 shows the contours of the vertical wave-number as a function of the horizontal wave-number and the angular frequency, at different depths. For the sake of comparison, the plotted variables are the adimensionned variables $\delta_x = k_x H$, $\delta_z = k_z H$ and $\log_{10}(\delta_\omega)$, where $\delta_\omega = \omega\sqrt{H/g}$.

Although Fig 4 is close to the one in the paper from Auclair *et al.* (2021), one can notice the influence of the ocean depth on the contours. This first result suggests that the variation of the parameters $\hat{c}_0, N, D$ with depth plays a non-negligible role in the waves dispersion.

# 6 Conclusion and future work

In this work we have presented a system describing the propagation of acoustic-gravity waves in a free-surface fluid over a varying bed (bathymetry) and with a variable sound speed, applicable to describe in particular hydro-acoustic and tsunami waves generated by earthquakes or landslides in the ocean. Through a rigorous linearization of the compressible Euler equation we have obtained a model able to represent many physical phenomena, such as the SOFAR channel or the propagation of internal waves. The variety of these phenomena is well represented in the dispersion relation.

In the derivation only few assumptions are made and some common simplifying hypotheses were avoided. In particular, the fluid is not necessarily assumed barotropic and it is not assumed irrotational. Thanks to this approach, many terms representing different physical phenomena are kept in the wave-like equation. With a numerical approximation one could then compute their respective magnitude, and justify which terms can be neglected. Note also that in the present work the source term is a displacement of the seabed, but this is not restrictive and other source terms could be used (a change in the surface pressure for example).

With additional assumptions compatible with the derivation of the system, such as considering a barotropic fluid, or restricting the model to the incompressible regime or to the acoustic regime, we are able to recover simpler models widely studied in the literature. The mathematical study of the more complete model can help gain insight on the other ones. For example, we could clearly identify the assumptions made in the hydro-acoustic waves model used by Stiassnie (2010); Sammarco *et al.* (2013) and others. Namely, in those models the fluid is assumed barotropic, and the effects of stratification and gravity are neglected inside the domain. The study of the more complete model also helped to write the conservation of energy in each simplified case. The linear model in Lagrangian coordinates can also be used to recover the linearized Euler equations in Eulerian coordinates. This brings a clear understanding of the usual – nevertheless non satisfactory – assumption that is used to derive the aforementioned
<figure><img src="image_6.png" /><figcaption>Figure 4: Contour of the vertical wave-number as a function of the horizon-
tal wave-number $\delta_x$ and the angular frequency $\delta_\omega$, at different depths $\xi^3$. 4a:
$\xi^3 = 2000$ m, 4b: $\xi^3 = 3600$ m, 4c: $\xi^3 = 4000$ m.</figcaption></figure>
models in Eulerian coordinates.

The wave-like formulation of the model makes it a good candidate for a numerical approximation by the finite elements method. The fact that it preserves an energy suggests that the problem is well-posed, which motivates a more thorough study of the mathematical problem. Numerical implementation of this model will make it possible to simulate acoustic-gravity waves generated by earthquakes and landslide sources accounting for the complex bathymetry, thus contributing to improve early-warning systems. It will also help quantifying the errors made in more simple models, such as the hypothesis of an irrotational flow. These two aspects will be investigated in a future work.

**Funding** This work was supported by grants from Région Ile-de-France and the project ERC-CG-2013-PE10-617472 SLIDEQUAKES and the European project DT-GEO (Digital Twin for GEOphysical Extremes.

**Declaration of interests** The authors report no conflict of interest.

**Acknowledgements** The authors thank the two reviewers for their thorough reading of the paper and for their useful suggestions.

## References

ABDOLALI, ALI, KADRI, USAMA & KIRBY, JAMES T. 2019 Effect of Water Compressibility, Sea-floor Elasticity, and Field Gravitational Potential on Tsunami Phase Speed. *Scientific Reports* **9** (1), 16874.

ABDOLALI, ALI, KIRBY, JAMES T. & BELLOTTI, GIORGIO 2015 Depth-integrated equation for hydro-acoustic waves with bottom damping. *Journal of Fluid Mechanics* **766**, R1.

ALLAIRE, GRÉGOIRE 2015 *Analyse numérique et optimisation, deuxième édition*. Les Éditions de l'École Polytechnique.

AUCLAIR, FRANCIS, DEBREU, LAURENT, DUVAL, EMILIE, MARCHESIELLO, PATRICK, BLAYO, ERIC & HILT, MARGOT 2021 Theory and analysis of acoustic-gravity waves in a free-surface compressible and stratified ocean. *Ocean Modelling* **168**, 31.

CAPLAN-AUERBACH, J., DZIAK, R. P., BOHNENSTIEHL, D. R., CHADWICK, W. W. & LAU, T.-K. 2014 Hydroacoustic investigation of submarine land-slides at West Mata volcano, Lau Basin: LANDSLIDES AT WEST MATA VOLCANO. *Geophysical Research Letters* **41** (16), 5927–5934.

CECIONI, C., BELLOTTI, G., ROMANO, A., ABDOLALI, A., SAMMARCO, P. & FRANCO, L. 2014 Tsunami Early Warning System based on Real-time Measurements of Hydro-acoustic Waves. *Procedia Engineering* **70**, 311–320.
<table><caption>Table 1: Summary of the equations, boundary conditions and energy for the general model and its simplifications. For conciseness the operator $\nabla_{\xi}$ and $\Delta_{\xi}$ are denoted respectively $\nabla$ and $\Delta$.</caption><thead><tr><th></th><th>Equation in &Omega;</th><th>Boundary conditions</th><th>Energy</th></tr></thead><tbody><tr><td>General model<br>(Sec. 2.3)</td><td>&rho;&#x0302;<sub>0</sub> &part;<sup>2</sup>&Ucirc;<sub>1</sub>/&part;t<sup>2</sup> &minus; &nabla; (&rho;<sub>0</sub>c&#x0302;<sub>0</sub><sup>2</sup>&nabla; &Ucirc;<sub>1</sub>) &minus; &rho;<sub>0</sub>g(&nabla;(&Ucirc;<sub>1</sub>&middot;e<sub>3</sub>) + &nabla; &Ucirc;<sub>1</sub>e<sub>3</sub>) = 0</td><td>&Ucirc;<sub>1</sub> &middot; n&#x0302;<sub>b</sub> = u&#x0302;<sub>b,1</sub>, &Ucirc;<sub>b</sub><br>&nabla; &Ucirc;<sub>1</sub> = 0, &Ucirc;<sub>s</sub></td><td>&epsilon; = <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Omega;</sub> &rho;<sub>0</sub> (<sup>|&nabla;U&#x0302;<sub>1</sub>|<sup>2</sup></sup>&frasl;<sub>&part;t</sub> + (&c&#x0302;<sub>0</sub>&nabla; &Ucirc;<sub>1</sub> &minus; <sup>g</sup>&frasl;<sub>c&#x0302;<sub>0</sub></sub> &Ucirc;<sub>1</sub>&middot;e<sub>3</sub>)<sup>2</sup>) + <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Omega;</sub> &rho;<sub>0</sub>N<sup>2</sup>(&Ucirc;<sub>1</sub>&middot;e<sub>3</sub>)<sup>2</sup> d&xi; + <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Gamma;<sub>s</sub></sub> &rho;<sub>0</sub>g(&Ucirc;<sub>1</sub>&middot;e<sub>3</sub>)<sup>2</sup> d&sigma;</td></tr><tr><td>Barotropic<br>case<br>(Sec. 3.1)</td><td>&rho;&#x0302;<sub>0</sub> &part;<sup>2</sup>&psi;/&part;t<sup>2</sup> &minus; &nabla; (&rho;<sub>0</sub>&nabla;&psi;) = 0, = 0,<br>&Ucirc;<sub>1</sub> = &nabla;&psi;</td><td>&nabla;&psi; &middot; n&#x0302;<sub>b</sub> = u&#x0302;<sub>b,1</sub>, &Ucirc;<sub>b</sub><br>&part;<sup>2</sup>&psi;/&part;t<sup>2</sup> + g&part;&psi;/&part;&xi;<sub>3</sub> = 0, &Ucirc;<sub>s</sub></td><td>&epsilon;<sub>bar</sub> = <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Omega;</sub> &rho;<sub>0</sub> (<sup>1</sup>&frasl;<sub>c&#x0302;<sub>0</sub><sup>2</sup></sub> (<sup>&part;&psi;</sup>&frasl;<sub>&part;t</sub>)<sup>2</sup> + |&nabla;&psi;|<sup>2</sup>) d&xi;<br>+ <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Gamma;<sub>s</sub></sub> &rho;<sub>0</sub> (<sup>&part;&psi;</sup>&frasl;<sub>&part;t</sub>)<sup>2</sup> d&sigma;</td></tr><tr><td>Incompressible<br>limit<br>(Sec. 3.2.2)</td><td>&rho;&#x0302;<sub>0</sub> &part;<sup>2</sup>&Ucirc;<sub>1,0</sub>/&part;t<sup>2</sup> &minus; g&rho;&#x0302;<sub>0</sub>(&Ucirc;<sub>1,0</sub>&middot;e<sub>3</sub>)e<sub>3</sub> + g&nabla;&phi;&#x0302;<sub>0</sub> = 0,<br>&nabla; &Ucirc;<sub>1,0</sub> = 0</td><td>&Ucirc;<sub>1,0</sub> &middot; n&#x0302;<sub>b</sub> = u&#x0302;<sub>b,1</sub>, &Ucirc;<sub>b</sub><br>&nabla; &Ucirc;<sub>1,0</sub> = 0, &Ucirc;<sub>s</sub><br>&phi;&#x0302;<sub>0</sub> = g&rho;<sub>0</sub>&Ucirc;<sub>1,0</sub>&middot;e<sub>3</sub>, &Ucirc;<sub>s</sub></td><td>&epsilon;<sub>incomp</sub> = <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Omega;</sub> (&rho;<sub>0</sub> |&part;&Ucirc;<sub>1,0</sub>/&part;t|<sup>2</sup> &minus; g&rho;&#x0302;<sub>0</sub>|&Ucirc;<sub>1,0</sub>&middot;e<sub>3</sub>|<sup>2</sup>) d&xi;<br>+ <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Gamma;<sub>s</sub></sub> g&rho;<sub>0</sub>|&Ucirc;<sub>1,0</sub>&middot;e<sub>3</sub>|<sup>2</sup> d&sigma;</td></tr><tr><td>Incompressible<br>limit,<br>homogeneous<br>density<br>(Sec.<br>3.2.2)</td><td>&Delta;&phi;&#x0302;<sub>0</sub> = 0,<br>&part;<sup>2</sup>&Ucirc;<sub>1,0</sub>/&part;t<sup>2</sup> = &minus;<sup>g</sup>&frasl;<sub>&rho;<sub>0</sub></sub> &nabla;&phi;&#x0302;<sub>0</sub></td><td>&nabla;&phi;&#x0302;<sub>0</sub> &middot; n&#x0302;<sub>b</sub> = &minus;&rho;<sub>0</sub>g u&#x0302;<sub>b,1</sub>, &Ucirc;<sub>b</sub><br>&part;<sup>2</sup>&phi;&#x0302;<sub>0</sub>/&part;t<sup>2</sup> + g&part;&phi;&#x0302;<sub>0</sub>/&part;&xi;<sub>3</sub> = 0, &Ucirc;<sub>s</sub></td><td>&epsilon;<sub>Poisson</sub> = <sup>1</sup>&frasl;<sub>2</sub> (<span style="display: inline-block; vertical-align: middle; text-align: center;"><span style="display: block; border-bottom: 1px solid black; padding: 0 0.2em;">(&int;<sub>&Omega;</sub> |&nabla;&phi;&#x0302;|<sup>2</sup> d&xi; + &int;<sub>&Gamma;<sub>s</sub></sub> <sup>1</sup>&frasl;<sub>g</sub> (<sup>&part;&rho;&#x0302;<sub>0</sub></sup>&frasl;<sub>&part;t</sub>)<sup>2</sup> d&sigma;</span><span style="display: block; padding: 0 0.2em;">&int;<sub>&Omega;</sub> |&nabla;&phi;&#x0302;|<sup>2</sup> d&xi; + |&nabla;&psi;&#x0302;<sub>0</sub>|<sup>2</sup> d&xi;</span></span>)</td></tr><tr><td>Acoustic limit<br>(Sec. 3.2.3)</td><td>&part;<sup>2</sup>&psi;&#x0302;<sub>0</sub>/&part;t<sup>2</sup> &minus; &rho;<sub>0</sub>c&#x0302;<sub>0</sub><sup>2</sup>&nabla; (&rho;&#x0302;<sub>0</sub><sup>&minus;1</sup>&nabla;&psi;&#x0302;<sub>0</sub>) = 0,<br>&Ucirc;<sub>1</sub> = &nabla;&psi;&#x0302;<sub>0</sub></td><td>&nabla;&psi;&#x0302;<sub>0</sub> &middot; n&#x0302;<sub>b</sub> = u&#x0302;<sub>b,1</sub>, &Ucirc;<sub>b</sub><br>&psi;&#x0302;<sub>0</sub> = 0, &Ucirc;<sub>s</sub></td><td>&epsilon;<sub>acoustic</sub> = <sup>1</sup>&frasl;<sub>2</sub> &int;<sub>&Omega;</sub> <sup>1</sup>&frasl;<sub>c&#x0302;<sub>0</sub><sup>2</sup></sub> (<sup>1</sup>&frasl;<sub>c&#x0302;<sub>0</sub><sup>2</sup></sub> (<sup>&part;&psi;&#x0302;<sub>0</sub></sup>&frasl;<sub>&part;t</sub>)<sup>2</sup> + |&nabla;&psi;&#x0302;<sub>0</sub>|<sup>2</sup>) d&xi;</td></tr></tbody></table>
CONSTANTIN, ADRIAN 2009 On the propagation of tsunami waves, with emphasis on the tsunami of 2004. *Discrete & Continuous Dynamical Systems - B* **12** (3), 525–537.

EWING, MAURICE, TOLSTOY, IVAN & PRESS, FRANK 1950 Proposed use of the T phase in tsunami warning systems*. *Bulletin of the Seismological Society of America* **40** (1), 53–58.

GILL, ADRIAN E. 1982 *Atmosphere-ocean dynamics*. *International geophysics series*. New York: Academic Press.

GIRAULT, VIVETTE & RAVIART, PIERRE-ARNAUD 1986 *Finite element methods for Navier-Stokes equations: theory and algorithms*, 1st edn. Springer series in computational mathematics 5. Berlin Heidelberg: Springer.

GODLEWSKI, EDWIGE, OLAZABAL, MARINA & RAVIART, PIERRE-ARNAUD 1999 On the linearization of systems of conservation laws for fluids at a material contact discontinuity. *Journal de Mathématiques Pures et Appliquées* **78** (10), 1013–1042.

GOMEZ, BERNABE & KADRI, USAMA 2021 Near real-time calculation of sub-marine fault properties using an inverse model of acoustic signals. *Applied Ocean Research* **109**, 102557.

GUYON, ETIENNE, ed. 2001 *Physical hydrodynamics*. Oxford ; New York: Oxford University Press, oCLC: ocm47048231.

HÄGG, LINUS & BERGGREN, MARTIN 2021 On the well-posedness of Galbrun’s equation. *Journal de Mathématiques Pures et Appliquées* **150**, 112–133.

INTERNATIONAL ASSOCIATION FOR THE PROPERTIES OF WATER AND STEAM 2009 IAPWS SR7-09, Supplementary Release on a Computationally Efficient Thermodynamic Formulation for Liquid Water for Oceanographic Use.

JENSEN, FINN B., KUPERMAN, WILLIAM A., PORTER, MICHAEL B. & SCHMIDT, HENRIK 2011 *Computational Ocean Acoustics*. New York, NY: Springer New York.

KADRI, USAMA & STIASSNIE, MICHAEL 2013 Generation of an acoustic-gravity wave by two gravity waves, and their subsequent mutual interaction. *Journal of Fluid Mechanics* **735**, R6.

KING, BENJAMIN, STONE, MARK, ZHANG, H. P., GERKEMA, THEO, MARDER, M., SCOTT, ROBERT B. & SWINNEY, HARRY L. 2012 Buoyancy frequency profiles and internal semidiurnal tide turning depths in the oceans. *Journal of Geophysical Research: Oceans* **117**.

KOMATITSCH, DIMITRI & TROMP, JEROEN 1999 Introduction to the spectral element method for three-dimensional seismic wave propagation. *Geophysical Journal International* **139** (3), 806–822.
KUEHNERT, J., MANGENEY, A., CAPDEVILLE, Y., MÉTAXIAN, J. P., BONILLA, L. F., STUTZMANN, E., CHALJUB, E., BOISSIER, P., BRUNET, C., KOWALSKI, P., LAURET, F. & HIBERT, C. 2020 Simulation of Topography Effects on Rockfall-Generated Seismic Signals: Application to Piton de la Fournaise Volcano. *Journal of Geophysical Research: Solid Earth* **125** (10).

LANNES, DAVID 2013 *The water waves problem: mathematical analysis and asymptotics. Mathematical surveys and monographs* volume 188. Providence, Rhode Island: American Mathematical Society.

LEGENDRE, GUILLAUME 2003 Rayonnement acoustique dans un fluide en écoulement: analyse mathématique et numérique de l'équation de Galbrun. PhD thesis, Paris VI.

LIGHTHILL, JAMES 1978 *Waves in fluids*. Cambridge University Press.

LONGUET-HIGGINS, M. S. 1950 A Theory of the Origin of Microseisms. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences* **243** (857), 1–35.

LOTTO, GABRIEL C. & DUNHAM, ERIC M. 2015 High-order finite difference modeling of tsunami generation in a compressible ocean from offshore earthquakes. *Computational Geosciences* **19** (2), 327–340.

MAEDER, MARCUS, GABARD, GWÉNAËL & MARBURG, STEFFEN 2020 90 Years of Galbrun’s Equation: An Unusual Formulation for Aeroacoustics and Hydroacoustics in Terms of the Lagrangian Displacement. *Journal of Theoretical and Computational Acoustics* **28** (04), 2050017.

NOSOV, M. A. & KOLESOV, S. V. 2007 Elastic oscillations of water column in the 2003 Tokachi-oki tsunami source: in-situ measurements and 3-D numerical modelling. *Natural Hazards and Earth System Sciences* **7** (2), 243–249.

NOUGUIER, FRÉDÉRIC, CHAPRON, BERTRAND & GUÉRIN, CHARLES-ANTOINE 2015 Second-order Lagrangian description of tri-dimensional gravity wave interactions. *Journal of Fluid Mechanics* **772**, 165–196.

SAMMARCO, P., CECIONI, C., BELLOTTI, G. & ABDOLALI, A. 2013 Depth-integrated equation for large-scale modelling of low-frequency hydroacoustic waves. *Journal of Fluid Mechanics* **722**, R6.

SMITH, JEROME A. 2015 Revisiting Oceanic Acoustic Gravity Surface Waves. *Journal of Physical Oceanography* **45** (12), 2953–2958.

STIASSNIE, MICHAEL 2010 Tsunamis and acoustic-gravity waves from underwater earthquakes. *Journal of Engineering Mathematics* **67** (1-2), 23–32.

STUTZMANN, E., ARDHUIN, F., SCHIMMEL, M., MANGENEY, A. & PATAU, G. 2012 Modelling long-term seismic noise in various environments: *Modelling seismic noise*. *Geophysical Journal International* **191** (2), 707–722.
TOLSTOY, IVAN 1950 The T Phase of Shallow-Focus Earthquakes. *Bulletin of the Seismological Society of America 1950* **40** (1), 25–51.

Appendix A

In this section an energy equation for the system (53) is obtained. Recall that the system (53) reads in $\hat{\Omega}$,

$$ \hat{\rho}_0 \frac{\partial^2 \hat{U}_1}{\partial t^2} - \nabla_{\xi} (\hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1) - (\nabla_{\xi} \hat{U}_1)^T \hat{\rho}_0 g + \hat{\rho}_0 \nabla_{\xi} \cdot \hat{U}_1 g = 0, \quad (182) $$

with the boundary conditions

$$ \left\{ \hat{U}_1 \cdot n_b = \hat{u}_{b,1} \quad \text{on } \hat{\Gamma}_b, \right. \quad (183) $$

$$ \left. \nabla_{\xi} \cdot \hat{U}_1 = 0 \quad \text{on } \hat{\Gamma}_s. \right. \quad (184) $$

By taking the scalar product of (53) with $\partial_t \hat{U}_1$ and integrating over the domain we have

$$ \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \hat{\rho}_0 \frac{\partial^2 \hat{U}_1}{\partial t^2} \right) d\xi - \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \nabla_{\xi} \left( \hat{\rho}_0 \hat{c}_0^2 \nabla_{\xi} \cdot \hat{U}_1 \right) \right) d\xi \\ + \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \nabla_{\xi} (\hat{U}_1 \cdot e_3) \hat{\rho}_0 g \right) d\xi - \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \hat{\rho}_0 \nabla_{\xi} \cdot \hat{U}_1 g e_3 \right) d\xi = 0. \quad (185) $$

For the first integral of (185) it holds

$$ \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \hat{\rho}_0 \frac{\partial^2 \hat{U}_1}{\partial t^2} \right) d\xi = \frac{d}{dt} \int_{\hat{\Omega}} \rho_0 \frac{1}{2} \left| \frac{\partial \hat{U}_1}{\partial t} \right|^2 d\xi. \quad (186) $$

The second term of (185) is integrated by parts, using $\nabla_{\xi} \cdot \hat{U}_1 = 0$ on the surface and $\hat{U}_1 \cdot \hat{n}_b = \hat{u}_{b,1}$ at the bottom (hence $\partial_t (\hat{U}_1 \cdot \hat{n}_b) = \partial_t \hat{u}_{b,1}$),

$$ - \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \nabla_{\xi} \left( \rho_0 c_0^2 \nabla_{\xi} \cdot \hat{U}_1 \right) d\xi = \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 \hat{c}_0^2 |\nabla_{\xi} \cdot \hat{U}_1|^2 d\xi \\ - \int_{\hat{\Gamma}_b} \rho_0 c_0^2 \nabla_{\xi} \cdot \hat{U}_1 \frac{\partial \hat{u}_{b,1}}{\partial t} d\sigma. \quad (187) $$

For the computation of the two last integral of (185), we define

$$ (I) = \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \nabla_{\xi} (\hat{U}_1 \cdot e_3) \hat{\rho}_0 g \right) d\xi - \int_{\hat{\Omega}} \frac{\partial \hat{U}_1}{\partial t} \cdot \left( \hat{\rho}_0 \nabla_{\xi} \cdot \hat{U}_1 g e_3 \right) d\xi, \quad (188) $$
an we denote by $\hat{n}_b$ be the vector normal to the boundary $\partial\Omega$. (I) is integrated by parts and reads

$$ (I) = \int_{\partial\Omega} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial \hat{U}_1}{\partial t} \cdot \hat{n}_b \, d\sigma - \int_{\hat{\Omega}} g \hat{U}_1 \cdot e_3 \nabla_{\xi} \cdot (\hat{\rho}_0 \frac{\partial \hat{U}_1}{\partial t}) \, d\xi \\ - \int_{\hat{\Omega}} \hat{\rho}_0 g \frac{\partial \hat{U}_1}{\partial t} \cdot e_3 \nabla_{\xi} \cdot \hat{U}_1 \, d\xi. \quad (189) $$

The boundary term is simplified using $\partial_t(\hat{U}_1 \cdot \hat{n}_b) = \partial_t \hat{u}_{b,1}$ at the bottom. On the boundary $\hat{\Gamma}_s$, the surface is horizontal hence the normal vector is the unit vector $e_3$, so it holds

$$ (I) = \int_{\hat{\Gamma}_b} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial \hat{u}_{b,1}}{\partial t} \, d\sigma + \int_{\hat{\Gamma}_s} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial \hat{U}_1}{\partial t} \cdot e_3 \, d\sigma \\ - \int_{\hat{\Omega}} g \hat{U}_1 \cdot e_3 \nabla_{\xi} \cdot (\hat{\rho}_0 \frac{\partial \hat{U}_1}{\partial t}) \, d\xi - \int_{\hat{\Omega}} \hat{\rho}_0 g \frac{\partial \hat{U}_1}{\partial t} \cdot e_3 \nabla_{\xi} \cdot \hat{U}_1 \, d\xi. \quad (190) $$

Next we develop the gradient in the third integral of Eq. (190). Note that $\hat{\rho}_0$ depends only on the vertical coordinate, then we have

$$ - \int_{\hat{\Omega}} g \hat{U}_1 \cdot e_3 \frac{\partial \hat{U}_1}{\partial t} \cdot \nabla_{\xi} \hat{\rho}_0 = - \int_{\hat{\Omega}} g \hat{U}_1 \cdot e_3 \frac{\partial \hat{U}_1 \cdot e_3}{\partial t} \frac{d\hat{\rho}_0}{d\xi^3} = - \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}'_0 g |\hat{U}_1 \cdot e_3|^2, \quad (191) $$

hence we obtain

$$ (I) = \int_{\hat{\Gamma}_b} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial \hat{u}_{b,1}}{\partial t} \, d\sigma + \frac{1}{2} \frac{d}{dt} \int_{\hat{\Gamma}_s} \hat{\rho}_0 g |\hat{U}_1 \cdot e_3|^2 \, d\sigma \\ - \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} g \hat{\rho}'_0 |\hat{U}_1 \cdot e_3|^2 \, d\xi - \int_{\hat{\Omega}} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial}{\partial t} (\nabla_{\xi} \cdot \hat{U}_1) \, d\xi \\ - \int_{\hat{\Omega}} \hat{\rho}_0 g \frac{\partial \hat{U}_1}{\partial t} \cdot e_3 \nabla_{\xi} \cdot \hat{U}_1 \, d\xi. \quad (192) $$

The two last terms of (192) are put together,

$$ (I) = \int_{\hat{\Gamma}_b} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \frac{\partial \hat{u}_{b,1}}{\partial t} \, d\sigma + \frac{1}{2} \frac{d}{dt} \int_{\hat{\Gamma}_s} \hat{\rho}_0 g |\hat{U}_1 \cdot e_3|^2 \, d\sigma \\ - \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} g \hat{\rho}'_0 |\hat{U}_1 \cdot e_3|^2 \, d\xi - \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 g \hat{U}_1 \cdot e_3 \nabla_{\xi} \cdot \hat{U}_1 \, d\xi. \quad (193) $$
Summing the terms (186), (187) and (193) yields

$$
\begin{align}
& \frac{d}{dt} \int_{\hat{\Omega}} \rho_0 \frac{1}{2} \left| \frac{\partial \hat{U}_1}{\partial t} \right|^2 d\xi + \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 \left( \hat{c}_0 \nabla_{\xi} \cdot \hat{U}_1 - \frac{g}{\hat{c}_0} \hat{U}_1 \cdot e_3 \right)^2 d\xi \nonumber \\
& - \frac{1}{2} \frac{d}{dt} \int_{\hat{\Omega}} \hat{\rho}_0 (\hat{U}_1 \cdot e_3)^2 \left( \frac{g^2}{\hat{c}_0^2} + \frac{g \hat{\rho}'_0}{\hat{\rho}_0} \right) d\xi + \frac{1}{2} \frac{d}{dt} \int_{\Gamma_s} \hat{\rho}_0 g (\hat{U}_1 \cdot e_3)^2 d\sigma \nonumber \\
& = \int_{\hat{\Gamma}_b} \rho_0 (c_0^2 \nabla_{\xi} \cdot \hat{U}_1 - \hat{\rho}_0 g \hat{U}_1 \cdot e_3) \frac{\partial \hat{u}_{b,1}}{\partial t} d\sigma, \tag{194}
\end{align}
$$

and by defining

$$
N_b = - \left( \frac{g^2}{\hat{c}_0^2} + \frac{g \hat{\rho}'_0}{\hat{\rho}_0} \right), \qquad (195)
$$

we obtain the energy equation (56).

**Appendix B**

In this section we derive the relations between the zero- and first-order approx-
imation in Eulerian and in Lagrangian coordinates, when differentiating with
respect to time or space. First note that $\phi_0$ and $\phi_1$ can be expressed in terms
of the displacement **d**. From the assumption of small displacements it holds
**d** = $\varepsilon$$\mathbf{d}$<sub>1</sub> + $\mathcal{O}(\varepsilon^2)$, then identifying the powers of $\varepsilon$ and summing yields

$$
\phi_{\varepsilon}(\xi, t) = \xi + \varepsilon d_1(\xi, t). \tag{196}
$$

From the change of coordinate we have

$$
\nabla_{\xi} \hat{X} = (\nabla_{\xi} \phi_{\varepsilon})^{T} \nabla X = (Id + \varepsilon \nabla_{\xi} \mathbf{d}_{1})^{T} \nabla X, \quad (197)
$$

and using this identity for $\hat{X} = \hat{X}_0 + \varepsilon \hat{X}_1$ yields

$$
\nabla_{\xi}(\hat{X}_{0} + \varepsilon \hat{X}_{1}) = \nabla X_{0} + \varepsilon ((\nabla_{\xi} \mathbf{d}_{1})^{T} \nabla X_{0} + \nabla X_{1}) + \mathcal{O}(\varepsilon^{2}). \quad (198)
$$

By identifying the powers of $\varepsilon$ it holds

$$
\nabla_{\xi} \hat{X}_0 = \nabla X_0, \quad \nabla_{\xi} \hat{X}_1 = (\nabla_{\xi} \mathbf{d}_1)^T \nabla X_0 + \nabla X_1. \tag{199}
$$

The same method is used for the time derivative. Starting with

$$
\frac{\partial \hat{X}}{\partial t}(\boldsymbol{\xi}, t) = \frac{\partial X}{\partial t}(\phi_{\varepsilon}(\boldsymbol{\xi}), t) + \frac{\partial \phi_{\varepsilon}}{\partial t}(\boldsymbol{\xi}) \cdot \nabla X(\phi_{\varepsilon}(\boldsymbol{\xi}), t), \quad (200)
$$

we obtain after replacing $X$ and $\hat{X}$ by their first order approximation,

$$
\frac{\partial \hat{X}_0}{\partial t} + \varepsilon \frac{\partial \hat{X}_1}{\partial t} = \frac{\partial X_0}{\partial t} + \varepsilon \left( \frac{\partial X_1}{\partial t} + \frac{\partial \mathbf{d}_1}{\partial t} \cdot \nabla X_0 \right) + \mathcal{O}(\varepsilon^2). \quad (201)
$$
With $\partial_t d_1(\xi, t) = \hat{U}_1(\xi, t) = U_1(x, t)$ it holds

$$ \frac{\partial X_0}{\partial t} + \varepsilon \left( \frac{\partial X_1}{\partial t} + U_1 \cdot \nabla X_0 \right) + \mathcal{O}(\varepsilon^2). \qquad (202) $$

We identify the powers of $\varepsilon$,

$$ \frac{\partial \hat{X}_0}{\partial t} = \frac{\partial X_0}{\partial t}, \quad \frac{\partial \hat{X}_1}{\partial t} = \frac{\partial X_1}{\partial t} + U_1 \cdot \nabla X_0. \qquad (203) $$