2.1.1. Equations in the volume

The unknowns are the fluid velocity ${\boldsymbol 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 $\boldsymbol{\mathsf{T}}$ has the form

(2.4)\begin{equation} \boldsymbol{\mathsf{T}} = ({-}p + \lambda \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U} ) \boldsymbol{\mathsf{I}} + 2 \mu \boldsymbol{\mathsf{D}}({\boldsymbol U}), \end{equation}

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

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

(2.5)$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho {\boldsymbol U}) = 0, \end{gather}$$

(2.6) $$\begin{gather} \frac{\partial}{\partial t} ( \rho {\boldsymbol U} ) + \boldsymbol{\nabla}\boldsymbol{\cdot}( \rho {\boldsymbol U}\otimes {\boldsymbol U}) + \boldsymbol{\nabla} p = \rho {\boldsymbol g} + \boldsymbol{\nabla} (\lambda \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (2 \mu \boldsymbol{\mathsf{D}}({\boldsymbol U})), \end{gather}$$

(2.7)$$\begin{gather}\hspace{-60pt}\frac{\partial }{\partial t} \left( \rho \frac{|{\boldsymbol U}|^2}{2} + \rho e \right) + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \left(\rho \frac{|{\boldsymbol U}|^2}{2} + \rho e + p\right) {\boldsymbol U} \right)\nonumber\\= \rho {\boldsymbol g}\boldsymbol{\cdot} {\boldsymbol U} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\lambda {\boldsymbol U} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (2 \mu \boldsymbol{\mathsf{D}}({\boldsymbol U}) \boldsymbol{\cdot} {\boldsymbol U}) - \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol q} . \end{gather}$$

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

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 Reference Gill1982). When considering $\rho$ and $s$ as independent, it is natural to introduce the scalar functions $f_e$, $f_p$ and $f_T$ satisfying

(2.8a–c)\begin{equation} e = f_e(\rho,s),\quad p = f_p(\rho,s),\quad T = f_T(\rho,s). \end{equation}

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

(2.9)\begin{equation} \frac{\partial e}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} e = \frac{f_p}{\rho^2} \left( \frac{\partial \rho}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho \right) +f_T \left( \frac{\partial s}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} s \right). \end{equation}

Using (2.7) $-{\boldsymbol U} \, \boldsymbol {\cdot }$ (2.6) and (2.9), one obtains as an intermediate step the evolution equation of the entropy,

(2.10)\begin{equation} \rho T \left(\frac{\partial s}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} s\right) = \lambda (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U})^2 + 2 \mu \boldsymbol{\mathsf{D}}({\boldsymbol U}) : \boldsymbol{\mathsf{D}}({\boldsymbol U}) - \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol q}. \end{equation}

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

(2.11)\begin{equation} \frac{\partial p}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} p = \frac{\partial f_p}{\partial \rho} \left( \frac{\partial \rho}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho \right) + \frac{1}{T \rho} \frac{\partial f_p}{\partial s} \left( \frac{\partial s}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} s \right), \end{equation}

hence using (2.5) and (2.10), one obtains

(2.12)\begin{align} \frac{\partial p}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} p ={-} \frac{\partial f_p}{\partial \rho} ( \rho \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U} ) + \frac{1}{T \rho} \frac{\partial f_p}{\partial s} ( \lambda (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U})^2 + 2 \mu \boldsymbol{\mathsf{D}}({\boldsymbol U}) : \boldsymbol{\mathsf{D}}({\boldsymbol U}) - \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol q}). \end{align}

At this point, we use the common assumption that the viscous term and the thermal dissipation can be neglected compared with the advection term (see Lannes Reference Lannes2013, Chap. 1). With (2.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) \geq 0$, hence we can introduce the speed of sound $c$ defined by

(2.13)\begin{equation} c^2 = \frac{\partial f_p}{\partial \rho} (\rho,s). \end{equation}

Equation (2.12) then reads

(2.14)\begin{equation} \frac{\partial p}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \boldsymbol{\nabla} p + \rho c^2 \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol U} = 0. \end{equation}

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

(2.15)\begin{equation} c^2(p,T) = \frac{\partial f_p}{\partial \rho} ( f_{\rho}(p,T), f_s(p,T)). \end{equation}

In practice, we choose to work directly with the expression $c=c(p,T)$ tabulated by International Association for the Properties of Water and Steam (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:

(2.16)$$\begin{gather} {\boldsymbol U} \boldsymbol{\cdot} {\boldsymbol n}_b = u_b = \partial_t b \quad \text{on } \varGamma_b, \end{gather}$$

(2.17)$$\begin{gather}p = p^a \quad \text{on } \varGamma_s. \end{gather}$$

The bottom boundary condition (2.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 ${\boldsymbol n}_b$ the unit vector normal to the bottom and oriented outwards. The second condition (2.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 a solution of the following kinematic equation:

(2.18)\begin{equation} \dfrac{\partial \eta}{\partial t} + {\boldsymbol U} \boldsymbol{\cdot} \begin{pmatrix} \partial_x \eta \\ \partial_y \eta \\ -1 \end{pmatrix}= 0\quad \text{on} \ \varGamma_s(t). \end{equation}

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,

(2.19)\begin{equation} \varOmega(0) = \{ (x,y,z) \in \mathbb{R}^3 \, | \,z_b(x,y) < z < H \}. \end{equation}

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

(2.20)$$\begin{gather} {\boldsymbol U}(x,y,z,0) = 0, \end{gather}$$

(2.21a–c)$$\begin{gather}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), \end{gather}$$

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 ${\boldsymbol U} \equiv 0$ if the functions $T_0, \rho _0, p_0$ satisfy

(2.22a–c)\begin{equation} \boldsymbol{\nabla} p_0 = \rho_0 {\boldsymbol g}, \quad \rho_0 = f_{\rho}(p_0, T_0), \quad p_0(H) = p^a. \end{equation}

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

(2.23)$$\begin{gather} \frac{{\rm d} p_0}{{\rm d} z} ={-}g f_{\rho}(p_0, T_0), \quad z \in (0, H), \end{gather}$$

(2.24)$$\begin{gather}p_0 = p^a, \quad z = H \end{gather}$$

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 (2.5), (2.6), (2.14) with boundary conditions (2.16), (2.17) and initial conditions (2.20), (2.21) will be linearized around the previously defined equilibrium state.

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

The system for the correction terms reads in $\hat \varOmega$,

(2.50)$$\begin{gather} \hat \rho_0 \frac{\partial \hat{\boldsymbol{U}}_1}{\partial t} + \nabla_{\xi} \hat p_1 - (\nabla_{\xi} {\boldsymbol d}_1)^{\rm T} \nabla_{\xi} \hat p_0 = \hat \rho_1 {\boldsymbol g}, \end{gather}$$

(2.51)$$\begin{gather}\frac{\partial \hat \rho_1}{\partial t} + \hat \rho_0 \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 = 0, \end{gather}$$

(2.52)$$\begin{gather}\frac{\partial \hat p_1}{\partial t} + \hat \rho_0 \hat c_0^2 \ \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 = 0, \end{gather}$$

with the boundary conditions

(2.53)$$\begin{gather} \hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b = \hat u_{b,1} \quad \text{on} \hat \varGamma_b, \end{gather}$$

(2.54)$$\begin{gather}\hat p_1 = 0 \quad \text{on } \hat \varGamma_s. \end{gather}$$

In this system, the speed of sound is evaluated at the limit – or background – pressure and temperature, $\hat c_0 = c(\hat p_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 (2.50) thanks to the other equations: differentiating in time (2.50) and replacing $\hat \rho _1$ and $\hat p_1$ with (2.51), (2.52), we obtain a second-order equation for $\hat {\boldsymbol {U}}_1$,

(2.55)\begin{equation} \hat \rho_0 \frac{\partial^2 \hat{\boldsymbol{U}}_1}{\partial t^2} - \nabla_{\xi} (\hat \rho_0 \hat c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1) - (\nabla_{\xi} \hat{\boldsymbol{U}}_1)^{\rm T} \hat \rho_0 {\boldsymbol g} + \hat \rho_0 \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 \ {\boldsymbol g} = 0 \ \text{ in } \hat{\varOmega} . \end{equation}

Using (2.52), the surface boundary condition (2.54) is formulated for $\hat {\boldsymbol {U}}_1$, hence the two boundary conditions for the wave-like equation (2.55) are

(2.56)$$\begin{gather} \hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b = \hat u_{b,1} \quad \text{on } \hat \varGamma_b, \end{gather}$$

(2.57)$$\begin{gather}\nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 = 0 \quad \text{on }\hat \varGamma_s. \end{gather}$$

The wave-like equation (2.55) is completed with vanishing initial condition for $\hat {\boldsymbol {U}}_1 (0)$ and $\partial _t \hat {\boldsymbol {U}}_1 (0)$. The system (2.55) 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 helioseismology and in aeroacoustics (Legendre Reference Legendre2003; Maeder, Gabard & Marburg Reference Maeder, Gabard and Marburg2020; Hägg & Berggren Reference Hägg and Berggren2021). However, to our knowledge, this equation has never been used to describe the propagation of hydro-acoustic waves.

We show now that the system (2.55), (2.56), (2.57) is energy preserving. A model describing a physical system should either preserve or dissipate 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 (2.55) with $\partial _t \hat {\boldsymbol {U}}_1$ and integrating over the domain. After some computations (see Appendix A), we have

(2.58)\begin{equation} \frac{{\rm d}}{{\rm d}t} \mathcal{E} = \int_{\hat \varGamma_b} \hat \rho_0 ( c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 - \hat \rho_0 g \hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3) \frac{\partial \hat u_{b,1}}{\partial t} \,\text{d} \sigma, \end{equation}

with the energy being the quadratic functional given by

(2.59)\begin{align} \mathcal{E} &= \int_{\hat \Omega} \rho_0 \frac{1}{2} \left| \frac{\partial \hat{\boldsymbol{U}}_1}{\partial t} \right|^2 \text{d} \boldsymbol \xi + \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 \left( \hat c_0 \nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_1 - \frac{g}{\hat c_0} \hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 \right)^2 \text{d} \boldsymbol \xi\nonumber\\ &\quad + \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 N_b (\hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3)^2 \,\text{d} \boldsymbol \xi + \frac{1}{2} \int_{\hat{\varGamma}_s} \hat \rho_0 g (\hat{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3)^2 \,\text{d} \sigma. \end{align}

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

(2.60)\begin{equation} 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). \end{equation}

The Brunt–Väisälä frequency, or buoyancy frequency, is closely related to the internal waves that appear in a stratified medium (Gill Reference Gill1982, Chap. 6). In the ocean, the usual values of $N_b$ are approximately $10^{-8} \ \text {rad}^2\ \text {s}^{-2}$ (King et al. Reference King, Stone, Zhang, Gerkema, Marder, Scott and Swinney2012).

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

(2.61)\begin{equation} \hat \rho_0 \frac{\partial^2 {\boldsymbol d}_1}{\partial t^2} - \nabla_{\xi} ( \hat \rho_0 \hat c_0^2 \nabla_{\xi} \boldsymbol{\cdot} {\boldsymbol d}_1 ) - (\nabla_{\xi} {\boldsymbol d}_1)^{\rm T} \hat \rho_0 {\boldsymbol g} + \hat \rho_0 \nabla_{\xi} \boldsymbol{\cdot} {\boldsymbol d}_1 {\boldsymbol g} = 0 \quad \text{in } \hat{\varOmega} . \end{equation}

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

(2.62)\begin{equation} \frac{{\rm d}}{{\rm d}t} \mathcal{E}_d = \int_{\hat \varGamma_b} \hat \rho_0 ( c_0^2 \nabla_{\xi} \boldsymbol{\cdot} {\boldsymbol d}_1 - \hat \rho_0 g {\boldsymbol d}_1 \boldsymbol{\cdot} {\boldsymbol e}_3) \hat u_{b,1} \,\text{d} \sigma, \end{equation}

with the energy

(2.63)\begin{align} \mathcal{E}_d &= \int_{\hat \Omega} \hat \rho_0 \frac{1}{2} | \hat{\boldsymbol{U}}_1 |^2 \,\text{d} \boldsymbol \xi + \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 \left( \hat c_0 \nabla_{\xi} \boldsymbol{\cdot} {\boldsymbol d}_1 - \frac{g}{\hat c_0} {\boldsymbol d}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 \right)^2 \text{d} \boldsymbol \xi\nonumber\\ &\quad + \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 N_b ({\boldsymbol d}_1 \boldsymbol{\cdot} {\boldsymbol e}_3)^2 \,\text{d} \boldsymbol \xi + \frac{1}{2} \int_{\hat{\varGamma}_s} \hat \rho_0 g ({\boldsymbol d}_1\boldsymbol{\cdot} {\boldsymbol e}_3)^2 \,\text{d} \sigma. \end{align}

The first term in (2.63) is the kinetic energy. We show that the second term of (2.63) corresponds to the acoustic energy. First, using (2.52) with the vanishing initial conditions yields $\hat \rho _0 \hat c_0^2 \boldsymbol {\cdot } \nabla _{\xi } {\boldsymbol d}_1 = - \hat p_1$. We define then the acoustic pressure

(2.64)\begin{equation} p_a = \hat p_1 - \boldsymbol{\nabla} \hat p_0 \boldsymbol{\cdot} {\boldsymbol d}_1. \end{equation}

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 (2.24), it holds that for the second term of (2.63),

(2.65)\begin{equation} \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 \left( \hat c_0 \nabla_{\xi} \boldsymbol{\cdot} {\boldsymbol d}_1 - \frac{g}{\hat c_0} {\boldsymbol d}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 \right)^2 \text{d} \boldsymbol \xi =\frac{1}{2} \int_{\hat \Omega} \frac{p_a^2}{\hat \rho_0 \hat c_0^2} \,\text{d} \boldsymbol \xi, \end{equation}

which is the usual expression for the acoustic energy (Lighthill Reference Lighthill1978). The last term of (2.63) is the potential energy associated with the surface waves. Finally, the third term of (2.63) is the potential energy associated with the internal gravity waves (Lighthill Reference Lighthill1978), under the condition

(2.66)\begin{equation} N_b > 0. \end{equation}

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'={\rm d}\hat \rho _0/{\rm d}z$ has to be negative and satisfy

(2.67)\begin{equation} \frac{|\hat \rho_0'|}{\hat \rho_0} > \frac{g}{\hat c_0^2}. \end{equation}

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 Reference Gill1982, 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.

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 Reference Constantin2009). 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

(3.18a,b)\begin{equation} \text{Fr} = \frac{U}{\sqrt{gH}}, \quad \text{Ma} = \frac{U}{C}. \end{equation}

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\ {\rm m\ s}^{-1}$, $U \sim 1 \ {\rm m\ s}^{-1}$ and $\sqrt {gH} \sim 100 \text { m s}^{-1}$. The dimensionless numbers are then

(3.19a,b)\begin{equation} \text{Fr} = 0.01, \quad \text{Ma} = 6.10^{{-}4}. \end{equation}

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 \varOmega$ and its surface and bottom boundary are respectively $\tilde \varGamma _s$ and $\tilde \varGamma _b$. The non-dimensional system reads, after simplification by the factor $\bar \rho U$,

(3.20)\begin{equation} \frac{\tilde \rho_0 }{\tau^2} \frac{\partial^2 \tilde{\boldsymbol{U}}_1}{\partial t^2} - \frac{C^2}{L^2} \nabla_{\xi} (\tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1) + \frac{g}{L} \tilde \rho_0 (\nabla_{\xi} ( \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1 {\boldsymbol e}_3 ) = 0, \end{equation}

with the boundary conditions

(3.21)$$\begin{gather} \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} \widetilde{{\boldsymbol n}_b} = \tilde u_{b,1} \quad \text{on } \tilde \varGamma_b, \end{gather}$$

(3.22)$$\begin{gather}\nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1 = 0 \quad \text{on } \tilde \varGamma_s, \end{gather}$$

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.

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}$. Equation (3.20) becomes

(3.23)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_1}{\partial t^2} - \frac{\text{Fr}}{\text{Ma}} \nabla_{\xi} ( \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1) + \tilde \rho_0 ( \nabla_{\xi} ( \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1 {\boldsymbol e}_3 )= 0. \end{equation}

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

(3.24)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_1}{\partial t^2} - \frac{1}{\delta^2} \nabla_{\xi} ( \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1) + \tilde \rho_0 ( \nabla_{\xi} ( \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1 {\boldsymbol e}_3 ) = 0, \end{equation}

and the goal is now to calculate the limit of (3.24) when $\delta$ goes to zero. We make the following ansatz for $\tilde {\boldsymbol {U}}_1$:

(3.25)\begin{equation} \tilde{\boldsymbol{U}}_1 =\tilde{\boldsymbol{U}}_{1, 0} +\delta^2 \tilde{\boldsymbol{U}}_{1, 2}+ {O}(\delta^3), \end{equation}

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

The equation obtained with the terms in $\delta ^{-2}$ reads

(3.26)\begin{equation} \nabla_{\xi} (\tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0}) = 0, \end{equation}

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

(3.27)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \nabla_{\xi} ( \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 2}) + \tilde \rho_0 ( \nabla_{\xi} ( \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0}\, {\boldsymbol e}_3 ) = 0 .\end{equation}

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

(3.28)$$\begin{gather} \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0} = 0 \quad \text{on } \tilde \varGamma_s, \end{gather}$$

(3.29)$$\begin{gather}\tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} \widetilde{{\boldsymbol n}}_b = \tilde u_{b,1} \quad \text{on } \tilde \varGamma_b. \end{gather}$$

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

(3.30)$$\begin{gather} \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 2} = 0 \quad \text{on } \tilde \varGamma_s, \end{gather}$$

(3.31)$$\begin{gather}\tilde{\boldsymbol{U}}_{1, 2} \boldsymbol{\cdot} \widetilde{{\boldsymbol n}}_b = 0 \quad \text{on } \tilde \varGamma_b. \end{gather}$$

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

(3.32)\begin{equation} \tilde{\boldsymbol{U}}_{1, 0} = \nabla_{\xi} \varphi_{1,0} + \nabla_{\xi} \times \boldsymbol \psi_{1,0}, \end{equation}

where $\varphi _{1,0}$ vanishes on $\tilde \varGamma _s$ and $\tilde \varGamma _b$. Injecting the decomposition of $\tilde {\boldsymbol {U}}_{1, 0}$ in (3.26) yields

(3.33)\begin{equation} \nabla_{\xi}( \tilde \rho_0 \tilde c_0^2 \ \Delta_{\xi} \varphi_{1,0} ) = 0, \end{equation}

hence the term inside the gradient is constant in space. Since the velocity $\tilde {\boldsymbol {U}}_{1, 0}$ is equal to zero at infinity, we obtain that $\Delta _{\xi } \varphi _{1,0} = 0$ in $\tilde \varOmega$ (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 \varOmega$. Then, taking the divergence of $\tilde {\boldsymbol {U}}_{1, 0}$ yields

(3.34)\begin{equation} \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0} = \nabla_{\xi} \boldsymbol{\cdot} (\nabla_{\xi} \times \boldsymbol \psi_{1,0}) = 0, \end{equation}

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

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

(3.35)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \nabla_{\xi} (\tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 2}) + \nabla_{\xi} ( \tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \tilde \rho_0 ' ( \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3 = 0. \end{equation}

Taking the curl of this equation yields

(3.36)\begin{equation} \nabla_{\xi} \times \left( \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \tilde \rho_0 ' ( \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3\right) = 0 . \end{equation}

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$,

(3.37)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \tilde \rho_0 ' ( \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3 ={-} \nabla_{\xi} \tilde \varphi_0. \end{equation}

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 (3.37) in (3.35), we have

(3.38)\begin{equation} \nabla_{\xi} ( - \tilde \varphi_0 - \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 2} + \tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) = 0, \end{equation}

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

(3.39)\begin{equation} \tilde \varphi_0 ={-} \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 2} + \tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3. \end{equation}

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

(3.40)\begin{equation} \tilde \varphi_0 = \tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 \quad \text{on } \tilde \varGamma_s. \end{equation}

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 {\boldsymbol {U}}_{1, 0} = U \tilde {\boldsymbol {U}}_{1, 0}$ satisfies

(3.41)$$\begin{gather} \hat \rho_0 \frac{\partial^2 \hat {\bf U}_{1,0}}{\partial t^2} - g \hat \rho_0' (\hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3) {\boldsymbol e}_3 + g \nabla_{\xi} \hat \varphi_0 = 0 \quad \text{in } \hat \Omega, \end{gather}$$

(3.42)$$\begin{gather}\nabla_{\xi} \boldsymbol{\cdot} \hat {\bf U}_{1,0} = 0 \quad \text{ in } \hat \Omega, \end{gather}$$

with the boundary conditions

(3.43)$$\begin{gather} \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b = \hat u_{b,1} \quad \text{on } \hat\varGamma_b, \end{gather}$$

(3.44)$$\begin{gather}\nabla_{\xi} \boldsymbol{\cdot} \hat{\boldsymbol{U}}_{1, 0} = 0 \quad \text{on } \hat\varGamma_s, \end{gather}$$

(3.45)$$\begin{gather}\hat \varphi_0 = \hat \rho_0 \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 \quad \text{on } \hat \varGamma_s. \end{gather}$$

We show that the model (3.41)–(3.45) preserves an energy. Taking the scalar product of (3.41) with $\partial _t \tilde {\boldsymbol {U}}_{1, 0}$ and integrating over $\hat \varOmega$ yields

(3.46)\begin{align} \frac{1}{2} \frac{{\rm d}}{{\rm d}t }\int_{\hat \Omega} \hat \rho_0 \left| \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \right| ^2 \text{d} \boldsymbol \xi -\int_{\hat \Omega} g \hat \rho_0 ' ( \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3 \boldsymbol{\cdot} \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \,\text{d} \boldsymbol \xi + \int_{\hat \Omega} g \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \boldsymbol{\cdot} \nabla_{\xi} \hat \varphi_0 \,\text{d} \boldsymbol \xi = 0. \end{align}

The last term of (3.46) is integrated by parts. With the vanishing divergence of $\hat {\boldsymbol {U}}_{1, 0}$ and the bottom condition (3.43), it holds that

(3.47)\begin{equation} \int_{\hat \Omega} g \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \boldsymbol{\cdot} \nabla_{\xi} \hat \varphi_0 \,\text{d} \boldsymbol \xi = \int_{\hat \varGamma_s} g \hat \varphi_0 \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \boldsymbol{\cdot} {\boldsymbol e}_3 \,\text{d} \sigma - \int_{\hat \varGamma_b} g \hat \varphi_0 \frac{\partial \hat u_{b,1}}{\partial t} \,\text{d} \sigma, \end{equation}

then $\hat \varphi _0$ is replaced in the surface integral using (3.45),

(3.48)\begin{align} &\frac{1}{2} \frac{{\rm d}}{{\rm d}t }\int_{\hat \Omega} \hat \rho_0 \left| \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \right| ^2 \text{d} \boldsymbol \xi - \int_{\hat \Omega} g \hat \rho_0 ' ( \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3 \boldsymbol{\cdot} \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \,\text{d} \boldsymbol \xi\nonumber\\ &\quad + \int_{\hat \varGamma_s} g \hat \rho_0 \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \boldsymbol{\cdot} {\boldsymbol e}_3 \,\text{d} \sigma = \int_{\hat \varGamma_b} g \hat \varphi_0 \frac{\partial \hat u_{b,1}}{\partial t} \,\text{d} \sigma. \end{align}

By defining the energy

(3.49)\begin{align} \mathcal{E}_{incomp} &= \frac{1}{2} \int_{\hat \Omega} \hat \rho_0 \left| \frac{\partial \hat{\boldsymbol{U}}_{1, 0} }{\partial t} \right| ^2 \text{d} \boldsymbol \xi - \frac{1}{2} \int_{\hat \Omega} g \hat \rho_0 ' | \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 |^2 \,\text{d} \boldsymbol \xi\nonumber\\ &\quad + \frac{1}{2} \int_{\hat \varGamma_s} g \hat \rho_0 |\hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3|^2 \,\text{d} \sigma , \end{align}

(3.46) can be formulated in the following way:

(3.50)\begin{equation} \frac{{\rm d} }{{\rm d}t} \mathcal{E}_{incomp} = \int_{\hat{\varGamma}_b} g \varphi_0 \frac{\partial \hat u_{b,1} }{\partial t} . \end{equation}

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 with 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$.

The system (3.42)–(3.41) represents an incompressible fluid. However, this system is different from the classical Poisson equation found in the literature (Lighthill Reference Lighthill1978) because of the assumption of a non-homogeneous 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 (3.41) yields

(3.51)\begin{equation} \Delta_{\xi} \hat \varphi_{0} = 0. \end{equation}

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

(3.52)\begin{equation} - \hat \rho_0 \frac{\partial^2 \hat u_{b,1}}{\partial t^2} - g \hat \rho_0 ' ( \hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3 ) {\boldsymbol e}_3 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b + g \nabla_{\xi} \hat \varphi_0 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b = 0. \end{equation}

For the surface condition, (3.45) is differentiated twice in time and the term in $\partial ^2_{tt} \hat {\boldsymbol {U}}_{1, 0}$ is replaced with (3.41),

(3.53)\begin{equation} \frac{\partial^2 \hat \varphi_0 }{\partial t^2} - g \hat \rho_0' (\hat{\boldsymbol{U}}_{1, 0} \boldsymbol{\cdot} {\boldsymbol e}_3) + g \frac{\partial \hat \varphi_0}{\partial \xi^3} = 0 \quad \text{ on } \tilde \varGamma_s. \end{equation}

With the assumption of a homogeneous density, the boundary conditions (3.52), (3.53) read then

(3.54)$$\begin{gather} \nabla_{\xi} \hat \varphi_0 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b ={-} \hat \rho_0 g \hat u_{b,1} \quad \text{on } \hat \varGamma_b, \end{gather}$$

(3.55)$$\begin{gather}\frac{\partial^2 \hat \varphi_0}{\partial t^2} + g \frac{\partial \hat \varphi_{0}}{\partial \xi^3} = 0 \quad \text{on } \hat \varGamma_s. \end{gather}$$

The Poisson equation (3.51) with boundary conditions (3.54)–(3.55) is the system satisfied by the velocity flow of an incompressible homogeneous free-surface fluid (Lighthill Reference Lighthill1978, Chap. 3.1). Note that it was required that $\tilde \rho _0' \neq 0$ in the system (2.55) to obtain an a priori positive energy. Here this assumption is dropped, however, a rather simple expression for the preserved energy can be derived: multiplying (3.51) by $\partial _t \hat \varphi _0$, integrating by parts and using (3.54)–(3.55), we obtain

(3.56)\begin{align} \int_{\hat \Omega} \Delta_{\xi} \hat \varphi_{0} \frac{\partial \hat \varphi_0}{\partial t} \,\text{d} \boldsymbol \xi &={-} \int_{\hat \Omega} \nabla_{\xi} \varphi \boldsymbol{\cdot} \nabla_{\xi} \left( \frac{\partial \hat \varphi_0}{\partial t}\right) \text{d} \boldsymbol \xi\nonumber\\ &\quad - \int_{\hat{\varGamma}_s} \frac{1}{g} \frac{\partial \hat \varphi_0}{\partial t} \frac{\partial^2 \hat \varphi_0}{\partial t^2} \,\text{d} \sigma + \int_{\hat{\varGamma}_b} \hat \rho_0 g \frac{\partial \hat \varphi_0}{\partial t} \hat u_{b,1} \,\text{d} \sigma. \end{align}

We define the energy

(3.57)\begin{equation} \mathcal{E}_{Poisson} = \frac{1}{2} \left( \int_{\hat \Omega} |\nabla_{\xi} \hat \varphi |^2 \,\text{d} \boldsymbol \xi + \int_{\hat \varGamma_s} \frac{1}{g} \left(\frac{\partial \hat \varphi_0}{\partial t} \right)^2 \text{d} \sigma \right). \end{equation}

Then it holds that

(3.58)\begin{equation} \frac{{\rm d}}{{\rm d}t} \mathcal{E}_{Poisson} ={-} \int_{\hat{\varGamma}_b} \hat \rho_0 g \frac{\partial \hat \varphi_0}{\partial t} \hat u_{b,1} \,\text{d} \sigma. \end{equation}

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

3.2.3. Acoustic limit

Another possible simplification of the system (2.55)–(2.57) 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 Reference Longuet-Higgins1950). Here we show that in the acoustic limit, the model reduces to a classical acoustic equation.

With the time scale $L/\tau = C$, corresponding to the acoustic wave, and with the same small parameter $\delta = \text {Ma}/\text {Fr}$ as before, the system (3.20) becomes

(3.59)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_1}{\partial t^2} - \nabla_{\xi} ( \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1) + \delta^2 \tilde \rho_0 ( \nabla_{\xi} ( \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} {\boldsymbol e}_3 ) - \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1\ {\boldsymbol e}_3) = 0 \quad \text{ in } \tilde \varOmega, \end{equation}

with the boundary conditions

(3.60)$$\begin{gather} \tilde{\boldsymbol{U}}_1 \boldsymbol{\cdot} \widetilde{{\boldsymbol n}}_b = \tilde u_{b,1} \quad \text{on } \tilde \varGamma_b, \end{gather}$$

(3.61)$$\begin{gather}\nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_1 = 0 \quad \text{on } \tilde \varGamma_s. \end{gather}$$

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

(3.62)\begin{equation} \tilde{\boldsymbol{U}}_1 =\tilde{\boldsymbol{U}}_{1, 0} +\delta^2 \tilde{\boldsymbol{U}}_{1, 2}+ {O}(\delta^3). \end{equation}

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

(3.63)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \nabla_{\xi} (\tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0}) = 0. \end{equation}

Taking the curl of this equation yields

(3.64)\begin{equation} \frac{\partial^2 }{\partial t^2} ( \nabla_{\xi} \times (\tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0})) = 0, \end{equation}

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

(3.65)\begin{equation} \tilde \rho_0 \tilde{\boldsymbol{U}}_{1, 0} = \nabla_{\xi} \tilde \psi_0. \end{equation}

By substituting in (3.63), we have

(3.66)\begin{equation} \nabla_{\xi} \left(\frac{\partial^2 \tilde \psi_0}{\partial t^2} - \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} ( \tilde \rho_0^{{-}1} \nabla_{\xi} \tilde \psi_0)\right) = 0, \end{equation}

then it holds that

(3.67)\begin{equation} \frac{\partial^2 \tilde \psi_0}{\partial t^2} - \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} ( \tilde \rho_0^{{-}1} \nabla_{\xi} \tilde \psi_0 ) = 0, \end{equation}

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

(3.68)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} - \frac{\partial}{\partial \xi^3} ( \tilde \rho_0 \tilde c_0^2 \nabla_{\xi} \boldsymbol{\cdot} \tilde{\boldsymbol{U}}_{1, 0} ){\boldsymbol e}_3 = 0 \quad \text{on } \tilde \varGamma_s. \end{equation}

Using the surface condition (3.61) in (3.68) yields

(3.69)\begin{equation} \tilde \rho_0 \frac{\partial^2 \tilde{\boldsymbol{U}}_{1, 0}}{\partial t^2} = 0 \quad \text{on } \tilde \varGamma_s. \end{equation}

With the definition of the potential $\tilde \psi _0$, it holds that

(3.70)\begin{equation} \frac{\partial^2 \boldsymbol{\nabla} \tilde \psi_0}{\partial t^2} = 0 \quad \text{on } \tilde \varGamma_s, \end{equation}

hence one has

(3.71)\begin{equation} \frac{\partial^2 \tilde \psi_0}{\partial t^2} = C(t) \quad \text{on } \tilde \varGamma_s, \end{equation}

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

(3.72)\begin{equation} \frac{\partial^2 \hat \psi_0}{\partial t^2} - \hat \rho_0 \hat c_0^2 \nabla_{\xi} \boldsymbol{\cdot} ( \hat \rho_0^{{-}1} \nabla_{\xi} \hat \psi_0) = 0 \quad \text{in } \hat \varOmega, \end{equation}

with the boundary conditions

(3.73)$$\begin{gather} \nabla_{\xi} \hat \psi_0 \boldsymbol{\cdot} \hat{\boldsymbol{n}}_b = \hat u_{b,1} \quad \text{on } \hat \varGamma_b, \end{gather}$$

(3.74)$$\begin{gather}\hat \psi_0 = 0 \quad \text{on } \hat \varGamma_s. \end{gather}$$

The system (3.72)–(3.74) 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 (3.72) by $\partial _t \psi / (\rho _0 \hat c_0^2)$ and integrating over the domain. The result reads after an integration by parts

(3.75)\begin{equation} \frac{{\rm d}}{{\rm d}t} \mathcal{E}_{acoustic} ={-} \int_{\hat \varGamma_b} \frac{1}{\hat \rho_0} \frac{\partial \hat \psi_0}{\partial t} \hat u_{b,1}\, \text{d} \sigma, \end{equation}

where the acoustic energy is

(3.76)\begin{equation} \mathcal{E}_{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\text{d} \boldsymbol \xi + \frac{1}{2} \int_{\hat \Omega} \frac{1}{\hat \rho_0} |\boldsymbol{\nabla} \hat \psi_0|^2 \,\text{d} \boldsymbol \xi. \end{equation}

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

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

Table 1. Summary of the equations, boundary conditions and energy for the general model and its simplifications. For conciseness, the operators $\nabla _{\xi }$ and $\Delta _\xi$ are denoted respectively $\boldsymbol {\nabla }$ and $\Delta$.

4.2.1. Boundary conditions and free surface description

Following the approach of Nouguier et al. (Reference Nouguier, Chapron and Guérin2015), we show that a description for the free surface can be obtained. In the following, the components of the fluid velocity are denoted ${\boldsymbol U}_1 = (U_1^1, U_1^2,U_1^3)^{\rm T}$. The surface is defined by $\varGamma _{s,eq} = \boldsymbol \phi _\varepsilon (\hat \varGamma _s)$, and we assume that at each time $t$, it can be parametrized as the graph $\eta _\varepsilon$. The elevation $\eta _\varepsilon$ is a function of $x,y$ and $t$ and can be decomposed in the following way:

(4.15)\begin{equation} \eta_\varepsilon(x,y,t) = H + \varepsilon \eta_1(x,y,t). \end{equation}

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

(4.16)\begin{equation} \phi_\varepsilon^3(\xi^1, \xi^2, H, t) = \eta_\varepsilon (x(\xi^1, \xi^2, H, t),\ y(\xi^1, \xi^2, H, t),t). \end{equation}

Differentiating (4.16) in time and using the (4.9) yields

(4.17)\begin{equation} \varepsilon \hat U_1^3(\xi^1, \xi^2, H,t) = \frac{\partial \eta_\varepsilon}{\partial t} + \varepsilon \hat U_1^1(\xi^1, \xi^2, H,t) \frac{\partial \eta_\varepsilon}{\partial x} + \varepsilon \hat U_1^2(\xi^1, \xi^2, H,t) \frac{\partial \eta_\varepsilon}{\partial y} . \end{equation}

We use the change of variables $\phi _\varepsilon (\boldsymbol \xi,t) = \boldsymbol{\xi} + \varepsilon \boldsymbol \phi _1(\boldsymbol \xi,t)$,

(4.18)\begin{align} \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}\nonumber\\ &\quad + \varepsilon U_1^2(\phi_\varepsilon(\xi^1, \xi^2, H,t), t) \frac{\partial \eta_\varepsilon}{\partial y} . \end{align}

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

(4.19)\begin{equation} U_1^3(x,y, H,t) = \frac{\partial \eta_1}{\partial t}, \end{equation}

which is the linearized equation for the free surface. Then the dynamic boundary conditions are linearized. With the change of variables, the boundary conditions (2.24), (2.53) and (2.54) become

(4.20)$$\begin{gather} {\boldsymbol U}_1 \boldsymbol{\cdot} {\boldsymbol n}_b =u_{b,1} \quad \text{on } \varGamma_{b,eq}, \end{gather}$$

(4.21)$$\begin{gather}p_0 = p^a \quad \text{on } \varGamma_{s,eq}, \end{gather}$$

(4.22)$$\begin{gather}p_1 = 0 \quad \text{on } \varGamma_{s,eq}. \end{gather}$$

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

(4.23)\begin{equation} p_0 + \varepsilon p_1 = p^a \quad \text{on } \varGamma_{s,eq}. \end{equation}

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

(4.24)\begin{equation} p_0(H) + \varepsilon (p_1(x,y,H,t) + p_0'(H) \eta_1) + {O}(\varepsilon^2) = p^a. \end{equation}

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

(4.25a,b)\begin{equation} p_0(H) = p^a, \quad p_1(x,y,H,t) = \rho_0(x,y,H,t) g \eta_1(x,y,t). \end{equation}

In a similar way, the linearization of (4.20) reads

(4.26)\begin{equation} 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). \end{equation}

Hence the equations for ${\boldsymbol U}_1, \rho _1, p_1$ can be fully defined on the domain $\hat \varOmega$, with an error in ${O}(\varepsilon ^2)$. Finally, note that the system (4.10)–(4.12) with the boundary conditions (4.25), (4.26) and the kinematic condition (4.19) is shown to be energy preserving, locally as well as over a whole water column (Lighthill Reference Lighthill1978; Lotto & Dunham Reference Lotto and Dunham2015).

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 § 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 that in Lagrangian coordinates, up to terms in ${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 with the study of the dispersion relation obtained from (2.55).