# Well-posedness and potential-based formulation for the propagation of hydro-acoustic waves and tsunamis

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, Jacques Sainte-Marie

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

To cite this version:

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, Jacques Sainte-Marie. Well-posedness and potential-based formulation for the propagation of hydro-acoustic waves and tsunamis. 2024. hal-04558414

HAL Id: hal-04558414

https://inria.hal.science/hal-04558414

Preprint submitted on 25 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Well-posedness and potential-based formulation for the propagation of hydro-acoustic waves and tsunamis

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, Jacques Sainte-Marie

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

## To cite this version:

Juliette Dubois, Sébastien Imperiale, Anne Mangeney, Jacques Sainte-Marie. Well-posedness and potential-based formulation for the propagation of hydro-acoustic waves and tsunamis. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 10.1051/m2an/2024076. hal-04558414v2

HAL Id: hal-04558414

https://inria.hal.science/hal-04558414v2

Submitted on 18 Oct 2024

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
# Well-posedness and potential-based formulation for the propagation of hydro-acoustic waves and tsunamis

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

April 25, 2024

## Abstract

We study a linear model for the propagation of hydro-acoustic waves and tsunami in a stratified free-surface ocean. A formulation was previously obtained by linearizing the compressible Euler equations. In this paper, we introduce a new formulation written with a generalized potential. The new formulation is obtained by studying the functional spaces and operators associated to the model. The mathematical study of this new formulation is easier and the discretization is also more efficient than for the previous formulation. We prove that both formulations are well posed and show that the solution to the first formulation can be obtained from the solution to the second. Finally, the formulations are discretized using a spectral element method, and we simulate tsunamis generation from submarine earthquakes and landslides.

## Introduction

We present and analyze a model describing the propagation of hydro-acoustic waves and tsunami. Hydro-acoustic waves are acoustic waves propagating in water. They are increasingly used to investigate movements of the seabed, as they can complete information brought by seismic waves [7, 17, 6]. Hydro-acoustic waves could improve early-warning systems of tsunamis generated by submarine earthquakes or landslides. To this aim, it is relevant to study models coupling acoustic waves and tsunamis.

The traditional model for hydro-acoustic waves and tsunamis combines a linear acoustic equation in the domain and the linearized free-surface equation of an incompressible fluid [12]. Another linear model based on an irrotational flow assumption was obtained with an Eulerian-Lagrangian approach [25]. However, the mathematical analysis for those models is lacking. In particular, it is not clear whether the systems preserve an energy. Energy preservation is a key element for ensuring stable numerical schemes. In a previous work [10] we derived a new formulation from the compressible Euler equation in Lagrangian coordinates. The resulting model consists in a wave-like equation and the system preserves an energy under a realistic condition on the background stratification.

The system can be seen as a particular case of the Galbrun equation with a vanishing mean flow and a non-homogeneous boundary condition of Dirichlet type. Most studies on the Galbrun equation focus on the harmonic regime, and the functional framework for its well-posedness is still subject to studies (see the review paper [26]). The literature on the analysis of the evolution problem is scarcer, but we can mention [2, 19]. The Galbrun equation with a uniform mean flow and with homogeneous boundary conditions is studied in [2] using a regularization method. In [19], it is shown that the solution to the Galbrun equation can be deduced from the solution to the linearized Euler equation.

In this work, we show existence and uniqueness for the velocity field solution to the linear compressible Euler equations with gravity, non-homogeneous boundary conditions and no mean flow. Moreover we propose

<footnote n="1">¹RWTH Aachen University, Math 111810, Templergraben 55, 52062 Aachen, Germany</footnote>

<footnote>*Corresponding author: dubois@eddy.rwth-aachen.de</footnote>

<footnote n="2">²Inria, M3DISIM Team, Inria Saclay, 1 rue Honoré d'Estienne d'Orves, 91120 Palaiseau, France</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">⁴Inria, Team ANGE, Inria Paris, 2 rue Simone Iff, 75012 Paris and Laboratoire Jacques-Louis Lions, Sorbonne Université, 75005 Paris, France</footnote>

<footnote n="5">⁵Laboratoire Jacques-Louis Lions, Sorbonne Université, 75005 Paris, France</footnote>
# WELL-POSEDNESS AND POTENTIAL-BASED FORMULATION FOR THE PROPAGATION OF HYDRO-ACOUSTIC WAVES AND TSUNAMIS

JULIETTE DUBOIS<fnref n="1" />, SÉBASTIEN IMPERIALE<fnref n="2" />, ANNE MANGENEY<fnref n="3" />,<fnref n="4" /> AND JACQUES SAINTE-MARIE<fnref n="4" />,<fnref n="5" />

**Abstract.** We study a linear model for the propagation of acoustic and surface gravity waves in a stratified free-surface ocean. A formulation was previously obtained by linearizing the compressible Euler equations. In this paper, we introduce a new formulation written with a generalized potential. The new formulation is obtained by studying the functional spaces and operators associated to the model. The mathematical study of this new formulation is easier and the discretization is also more efficient than for the previous formulation. We prove that both formulations are well-posed and show that the solution to the first formulation can be obtained from the solution to the second. Finally, the formulations are discretized using a spectral element method, and we simulate tsunamis generation from submarine earthquakes and landslides.

**Résumé.** Nous étudions un modèle linéaire pour la propagation des ondes acoustiques et des vagues dans un océan stratifié à surface libre. Une première formulation, déduite des équations d’Euler compressible linéarisées, a été proposée dans un précédent article. Dans cet article, nous développons une nouvelle formulation écrite avec un potentiel généralisé. Cette nouvelle formulation est obtenue en étudiant les espaces fonctionnels et les opérateurs associés au modèle. L’étude mathématique et la discrétisation de cette nouvelle formulation sont plus simples que celles de la précédente formulation. Nous montrons que les deux formulations sont bien posées, et que la solution de la première formulation peut être calculée avec la solution de la seconde. Les deux formulations sont ensuite discrétisées avec une méthode d’éléments finis spectraux. Nous présentons des simulations de tsunamis générés par des tremblement de terre et des glissements de terrain sous-marins.

**2020 Mathematics Subject Classification.** 35L10, 65M70, 76N30, 35Q86.

...

## INTRODUCTION

We present and analyze a model describing the propagation of hydro-acoustic waves and tsunami. Hydro-acoustic waves are acoustic waves propagating in water. They are increasingly used to investigate movements

<footnote>*Keywords and phrases:* Hydro-acoustic waves, Tsunami, Spectral element method, Linearized Euler equation</footnote>

<footnote n="1">¹ RWTH Aachen University, Math 111810, Templergraben 55, 52062 Aachen, Germany</footnote>

<footnote n="2">² Inria, M3DISIM Team, Inria Saclay, 1 rue Honoré d'Estienne d'Orves, 91120 Palaiseau, France</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">⁴ Inria, Team ANGE, Inria Paris, 2 rue Simone Iff, 75012 Paris, France</footnote>

<footnote n="5">⁵ Laboratoire Jacques-Louis Lions, Sorbonne Université, 75005 Paris, France</footnote>
a set of PDE, involving two scalar potentials, that can be solved instead of finding a solution to the linearized compressible Euler equations. The mathematical study is easier and the discretization more efficient that the velocity-field formulation. Finally, we present several numerical results and simulation to assess the properties of the introduced models. The simulations reproduce hydro-acoustic waves and tsunamis generated by earthquakes and submarine landslides. We show that the interference pattern characteristics of landslide-generated hydro-acoustic waves is reproduced by the model, providing a unique tool to detect and characterize the landslide source.

In more details, the article is organized as follows. We first continue this introduction by recalling briefly the **velocity-field formulation** obtained in [10] and present the new formulation (denoted here **potential-based formulation**) just after. In Section 1, the functional spaces and the operators used for studying the two mentioned formulations are introduced. In Section 2, we study the **potential-based formulation**. Even though this formulation is new, we start with this one because the analysis is much more direct than the analysis of **the velocity-field formulation**. The latter is then studied in Section 3. The main difficulty of this section lies in finding an adapted lifting operator for the boundary term. The correspondence between the two formulations is discussed in Section 4; we prove that the solution to the **velocity-field formulation** can be computed using the solution to **the potential-based formulation** (the reciprocal is formally true but of lesser interest and is not studied in detail). In Section 5 we describe how the two formulations are discretized using spectral finite elements [8, 22] and present several numerical experiments. A first set of simulations helps validate the model by comparison with the literature. It also illustrates the correspondence between the velocity-field formulation and the potential-based formulation. A second set of simulations investigates the classical hypothesis of irrotational flow using the correspondence between the velocity-field and potential-based formulations. The last set of simulations is a preliminary work towards the simulation of hydro-acoustic waves generated by landslides. We show that the proposed model is able to recover an interference pattern in the simulated hydro-acoustic waves, characteristic of the Lloyd-mirror effect that occurs during submarine landslides.

Readers more interested in the applications than in the mathematical details regarding the **potential-based formulation** could only read this introductory section and directly go to Section 5.

## Velocity-field formulation

We denote by $\Omega$ the domain representing an ocean at rest. The coordinates of $\Omega$ are written $(\mathbf{x}, z)$, with the horizontal coordinate $\mathbf{x} \in \mathbb{R}^{d-1}$, $d = 2$ or $d = 3$. The domain is unbounded in the horizontal direction and bounded in the vertical direction, with a fixed surface at $z = H$ and a time-independent topography $z_b(\mathbf{x})$ at the bottom, see Figure 1. The domain is written

$$ \Omega = \{(\mathbf{x}, z) | \mathbf{x} \in \mathbb{R}^{d-1}, z_b(\mathbf{x}) \leq z \leq H\}, $$

and its boundary is denoted $\Gamma = \partial\Omega$. The topography is assumed to satisfy the following conditions: $z_b \in W^{1,\infty}(\mathbb{R}^{d-1})$ and there exists positive scalars $(H_-, R_+)$ such that,

$$ 0 \leq z_b(\mathbf{x}) \leq H_- < H \quad \text{and} \quad \nabla z_b(\mathbf{x}) = 0 \text{ for } |\mathbf{x}| \geq R_+. $$

These properties ensure that the domain $\Omega$ does not degenerate and is Lipschitz. The surface and bottom boundaries are respectively

$$ \Gamma_s = \{(\mathbf{x}, H) | \mathbf{x} \in \mathbb{R}^{d-1}\} \quad \text{and} \quad \Gamma_b = \{(\mathbf{x}, z) | \mathbf{x} \in \mathbb{R}^{d-1}, z = z_b(\mathbf{x})\}. $$

The considered set of partial differential equation has been introduced in [10]; it consists in a linear system of equations for the fluid velocity $\mathbf{U}(\mathbf{x}, z, t)$ and reads

$$ \rho_0 \frac{\partial^2 \mathbf{U}}{\partial t^2} - \nabla (\rho_0 c_0^2 \nabla \cdot \mathbf{U} - \rho_0 g \mathbf{U} \cdot \mathbf{e}_z) - \nabla \cdot (\rho_0 g \mathbf{U}) \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T]. \quad (0.1) $$

Equation (0.1) represents the propagation of the acoustic, internal gravity and surface gravity waves in a stratified fluid. The vector $\mathbf{e}_z = (0 \ 0 \ 1)^t$ is the unit vector along the z-axis. The constant scalar $g > 0$ is the acceleration
of the seabed, as they can complete information brought by seismic waves [6, 7, 17]. Hydro-acoustic waves could improve early-warning systems of tsunamis generated by submarine earthquakes or landslides. To this aim, it is relevant to study models coupling acoustic waves and tsunamis.

The traditional model for hydro-acoustic waves and tsunamis combines a linear acoustic equation in the domain and the linearized free-surface equation of an incompressible fluid [13]. Another linear model based on an irrotational flow assumption was obtain with an Eulerian-Lagrangian approach [25]. However, the mathematical analysis for those models is lacking. In particular, it is not clear whether the systems preserve an energy. Energy preservation is a key element for ensuring stable numerical schemes. In a previous work [11] we derived a new formulation from the compressible Euler equation in Lagrangian coordinates. The resulting model consists in a wave-like equation and the system preserves an energy under a realistic condition on the background stratification.

The system can be seen as a particular case of the Galbrun equation with a vanishing mean flow and a non-homogeneous boundary condition of Dirichlet type. Most studies on the Galbrun equation focus on the harmonic regime, and the functional framework for its well-posedness is still subject to studies (see the review paper [26] and [19]). The literature on the analysis of the evolution problem is scarcer, but we can mention [2,20]. The Galbrun equation with a uniform mean flow and with homogeneous boundary conditions is studied in [2] using a regularization method. In [20], it is shown that the solution to the Galbrun equation can be deduced from the solution to the linearized Euler equation.

In this work, we show existence and uniqueness for the velocity field solution to the linear compressible Euler equations with gravity, non-homogeneous boundary conditions and no mean flow. Moreover we propose a set of PDE, involving two scalar potentials, that can be solved instead of finding a solution to the linearized compressible Euler equations. The mathematical study is easier and the discretization more efficient that the velocity-field formulation. Finally, we present several numerical results and simulation to assess the properties of the introduced models. The simulations reproduce hydro-acoustic waves and tsunamis generated by earthquakes and submarine landslides. We show that the interference pattern characteristics of landslide-generated hydro-acoustic waves is reproduced by the model, providing a unique tool to detect and characterize the landslide source.

In more details, the article is organized as follows. We first continue this introduction by recalling briefly **the velocity-field formulation** obtained in [11] and present the new formulation (named here **potential-based formulation**) just after. In Section 1, the functional spaces and the operators used for studying the two mentioned formulations are introduced. In Section 2, we study **the potential-based formulation**. Even though this formulation is new, we start with this one because the analysis is much more direct than the analysis of **the velocity-field formulation**. The latter is then studied in Section 3. The main difficulty of this section lies in finding an adapted lifting operator for the boundary term. The correspondence between the two formulations is discussed in Section 4; we prove that the solution to **the velocity-field formulation** can be computed using the solution to **the potential-based formulation** (the reciprocal is formally true but of lesser interest and is not studied in detail). In Section 5 we describe how the two formulations are discretized using spectral finite elements [8, 23] and present several numerical experiments. A first set of simulations helps validate the model by comparison with the literature. It also illustrates the correspondence between the velocity-field formulation and the potential-based formulation. A second set of simulations investigates the classical hypothesis of irrotational flow using the correspondence between the velocity-field and potential-based formulations. The last set of simulations is a preliminary work towards the simulation of hydro-acoustic waves generated by landslides. We show that the proposed model is able to recover an interference pattern in the simulated hydro-acoustic waves, characteristic of the Lloyd-mirror effect that occurs during submarine landslides.

Readers more interested in the applications than in the mathematical details regarding the **potential-based formulation** could only read this introductory section and directly go to Section 5.
of gravity. The parameters $\rho_0$ and $c_0$ are respectively the fluid density at equilibrium and the sound speed. They depend on $z$ only, and we make the following assumptions,

$$ \rho_0 \in C^1([0, H]) \quad \text{and} \quad c_0 \in C^0([0, H]). \qquad (0.2) $$

For the well-posedness of the problem, the velocity and the density must satisfy some usual non-degeneracy properties. We assume that there exists positive scalars ($\rho_-$, $\rho_+$) and ($c_-$, $c_+$) such that

$$ \rho_- \le \rho_0(z) \le \rho_+ \quad \text{and} \quad c_- \le c_0(z) \le c_+. $$

Following the literature [15], we also defined the scalar field $N^2$ as

$$ N^2(z) = - \frac{g}{\rho_0(z)} \frac{d\rho_0(z)}{dz} - \frac{g^2}{c_0^2(z)}. \qquad (0.3) $$

Such field is called the Brunt-Väisälä frequency, or buoyancy frequency. The case $N^2 < 0$ corresponds to a fluid that is denser above and lighter below, hence it is an unstable equilibrium. Since the equations were obtained by linearizing around a stable equilibrium, we assume in the following $N^2 \ge 0$. The case $N^2 = 0$ corresponds to the case of a constant temperature, in which case the fluid is called barotropic.

Equation (0.1) is completed with boundary conditions. On the bottom $\Gamma_b$ we consider a localized displacement of amplitude $u_b$ of the seabed, caused for example by a submarine earthquake, landslide or caldera collapse. On the surface $\Gamma_s$ a stress-free boundary condition is applied. Here the stress is a pressure [10], and is proportional to the divergence of the velocity field; hence the boundary conditions read

$$ \mathbf{U} \cdot \mathbf{n}_b = u_b \quad \text{on } \Gamma_b \times [0, T], \quad \nabla \cdot \mathbf{U} = 0 \quad \text{on } \Gamma_s \times [0, T], \qquad (0.4) $$

where $\mathbf{n}_b$ is the outward unitary normal of the domain $\Omega$ on $\Gamma_b$. The regularity of the displacement $u_b$ will be stated later. Finally, in our context it is relevant to choose vanishing initial conditions, which also simplifies the forthcoming analysis,

$$ \mathbf{U}(\mathbf{x}, z, 0) = 0, \quad \frac{\partial \mathbf{U}}{\partial t}(\mathbf{x}, z, 0) = 0 \quad \text{on } \Omega. \qquad (0.5) $$

## Potential-based formulation

Problem (0.1)-(0.4) can be written as an abstract wave equation using an unbounded linear operator $G$,

$$ \frac{d^2 \mathbf{U}}{dt^2} + G^* \tilde{G} \mathbf{U} = 0, \qquad (0.6) $$

where $G^*$ is the adjoint of $G$ and $\tilde{G}$ is an extension of $G$. The introduction of this non-symmetric formulation using the extension $\tilde{G}$ is motivated by the presence of a non-homogeneous essential condition in the boundary conditions ((0.4)). The expression (0.6) will be obtained rigorously in Section 1. One originality of this work is to construct and analyze a “dual” or “adjoint” wave-like problem for a new unknown $\Phi$, satisfying

$$ \frac{d^2 \Phi}{dt^2} + \tilde{G} G^* \Phi = 0. \qquad (0.7) $$

Such problem is shown to be equivalent – in a sense given rigorously later – to the problem (0.6). Moreover, it presents several advantages from mathematical and numerical perspectives. The new unknown $\Phi$ has three scalar components, $\Phi = (\varphi, \psi, \gamma)^t$, where $\varphi$ and $\psi$ are scalar fields and $\gamma = \varphi|_{\Gamma_s}$. Later, we show that $\Phi$ is related to the velocity $\mathbf{U}$ by the formula

$$ \mathbf{U} = G^* \Phi = -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z, \qquad (0.8) $$
Velocity-field formulation

We denote by Ω the domain representing an ocean at rest. The coordinates of Ω are written (x, z), with the horizontal coordinate x ∈ ℝ$^{d-1}$, d = 2 or d = 3. The domain is unbounded in the horizontal direction and bounded in the vertical direction, with a fixed surface at z = H and a time-independent topography z<sub>b</sub>(x) at the bottom; see Figure 1. The domain is written

$$
\Omega = \{(\mathbf{x}, z) | \mathbf{x} \in \mathbb{R}^{d-1}, z_b(\mathbf{x}) \leq z \leq H\},
$$

and its boundary is denoted by Γ = ∂Ω. The topography is assumed to satisfy the following conditions:
z<sub>b</sub> ∈ W<sup>1,∞</sup>(ℝ<sup>d-1</sup>) and there exists positive scalars (H<sub>-</sub>, R<sub>+</sub>) such that,

$$
0 \le z_b(\mathbf{x}) \le H_- < H \quad \text{and} \quad \nabla z_b(\mathbf{x}) = 0 \text{ for } |\mathbf{x}| \ge R_+.
$$

These properties ensure that the domain Ω does not degenerate and is Lipschitz. The surface and bottom boundaries are respectively

$$
\Gamma_s = \{(\mathbf{x}, H) | \mathbf{x} \in \mathbb{R}^{d-1}\} \quad \text{and} \quad \Gamma_b = \{(\mathbf{x}, z) | \mathbf{x} \in \mathbb{R}^{d-1}, z = z_b(\mathbf{x})\}.
$$

The considered set of partial differential equations has been introduced in [11]; it consists in a linear system of equations for the fluid velocity U(x, z, t) and reads

$$
\rho_0 \frac{\partial^2 \mathbf{U}}{\partial t^2} - \nabla (\rho_0 c_0^2 \nabla \cdot \mathbf{U} - \rho_0 g \mathbf{U} \cdot \mathbf{e}_z) - \nabla \cdot (\rho_0 g \mathbf{U}) \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T]. \quad (0.1)
$$

Equation (0.1) represents the propagation of the acoustic, internal gravity and surface gravity waves in a stratified fluid. The vector e<sub>z</sub> = (0 0 1)<sup>T</sup> is the unit vector along the z-axis. The constant scalar g > 0 is the acceleration of gravity. The parameters ρ<sub>0</sub> and c<sub>0</sub> are respectively the fluid density at equilibrium and the sound speed. They depend on z only, and we make the following assumptions:

$$
\rho_0 \in C^1([0, H]) \quad \text{and} \quad c_0 \in C^0([0, H]). \tag{0.2}
$$

For the well-posedness of the problem, the velocity and the density must satisfy some usual non-degeneracy properties. We assume that there exists positive scalars (ρ−, ρ+), and (c−, c+), such that

$$
\rho_{-} \leq \rho_{0}(z) \leq \rho_{+} \quad \text{and} \quad c_{-} \leq c_{0}(z) \leq c_{+}.
$$

Following the literature [15], we also defined the scalar field N² as

$$
N^2(z) = - \frac{g}{\rho_0(z)} \frac{d\rho_0(z)}{dz} - \frac{g^2}{c_0^2(z)}. \qquad (0.3)
$$

Such field is called the Brunt-Väisälä frequency, or buoyancy frequency. The case N² < 0 corresponds to a fluid
that is denser above and lighter below, hence it is an unstable equilibrium. Since the equations were obtained
by linearizing around a stable equilibrium, we assume in the following N² ≥ 0. The case N² = 0 corresponds
to the case of a constant temperature, in which case the fluid is called barotropic.

Equation (0.1) is completed with boundary conditions. On the bottom Γ<sub>b</sub> we consider a localized displacement of amplitude u<sub>b</sub> of the seabed, caused for example by a submarine earthquake, landslide or caldera collapse. On the surface Γ<sub>s</sub> a stress-free boundary condition is applied. Here the stress is a pressure [11], and is proportional to the divergence of the velocity field; hence the boundary conditions read

$$
\mathbf{U} \cdot \mathbf{n}_b = u_b \quad \text{on } \Gamma_b \times [0, T], \quad \nabla \cdot \mathbf{U} = 0 \quad \text{on } \Gamma_s \times [0, T], \qquad (0.4)
$$
hence it can be seen as a generalized potential. We also show that $\varphi$ and $\psi$ satisfy a second set of coupled partial differential equations describing the same physical system as (0.1)-(0.4), more precisely

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) + g \frac{\partial \varphi}{\partial z} - gN \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \quad (0.9) $$

$$ \frac{\partial^2 \psi}{\partial t^2} - N \frac{\partial \varphi}{\partial z} + N^2 \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T]. \quad (0.10) $$

The system is completed with boundary conditions. Using (0.4) and (0.8), we find on $\Gamma_b$

$$ \mathbf{U} \cdot \mathbf{n}_b = -\nabla \varphi \cdot \mathbf{n}_b + N \left( \psi + \frac{N}{g} \varphi \right) (\mathbf{e}_z \cdot \mathbf{n}_b) = u_b, \quad \text{on } \Gamma_b \times [0, T]. \quad (0.11) $$

This condition is a Robin boundary conditions for the potentials $(\varphi, \psi)$. It is easier to consider than condition (0.4) for both the analysis and the discretization. Finally, to obtain the boundary condition on $\Gamma_s$ we observe that (0.9) can be rewritten, using (0.8),

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \mathbf{U} - g\mathbf{U} \cdot \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T]. $$

Formally, evaluating the equation above on the boundary $\Gamma_s$ and using (0.4) we obtain a boundary condition involving a second-order time derivative that accounts for surface gravity waves,

$$ \frac{\partial^2 \varphi}{\partial t^2} - g\mathbf{U} \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T], \quad (0.12) $$

where $\mathbf{n}_s$ is the outward unitary normal of the domain $\Omega$ on $\Gamma_s$. The initial conditions are deduced from (0.5), they read

$$ \varphi(\mathbf{x}, z, 0) = \frac{\partial \varphi}{\partial t}(\mathbf{x}, z, 0) = 0, \quad \psi(\mathbf{x}, z, 0) = \frac{\partial \psi}{\partial t}(\mathbf{x}, z, 0) = 0 \quad \text{on } \Omega. \quad (0.13) $$

When $d = 3$, the system (0.9)-(0.13) involves only two scalar fields $(\varphi$ and $\psi)$, compared to three for the velocity based formulation (0.1). Moreover, this system is a generalization to the one commonly introduced in hydrodynamics [15]. Indeed, when considering the barotropic case, we have $N = 0$ and the system of partial differential equations (0.9)-(0.10) reduces to

$$ \frac{\partial^2 \varphi}{\partial t^2} - c_0^2 \Delta \varphi + g \frac{\partial \varphi}{\partial z} = 0, \quad \text{in } \Omega \times [0, T], \quad (0.14) $$

$$ \frac{\partial^2 \psi}{\partial t^2} = 0, \quad \text{in } \Omega \times [0, T]. \quad (0.15) $$

Because of the vanishing initial conditions, we have that $\psi = 0$. The system is then described by the function $\varphi$ only. The boundary conditions are also simplified,

$$ \frac{\partial^2 \varphi}{\partial t^2} + g \nabla \varphi \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T] \quad \text{and} \quad \nabla \varphi \cdot \mathbf{n}_b = u_b, \quad \text{on } \Gamma_b \times [0, T]. \quad (0.16) $$

This simplified system corresponds to the system of equations that was already introduced in the literature [10, 25].

# 1 Preliminary definitions

In this section, we define functional spaces and operators that will be used throughout the paper.
where $\mathbf{n}_b$ is the outward unitary normal of the domain $\Omega$ on $\Gamma_b$. The regularity of the displacement $u_b$ will be stated later. Finally, in our context it is relevant to choose vanishing initial conditions, which also simplifies the forthcoming analysis,

$$ \mathbf{U}(\mathbf{x}, z, 0) = 0, \qquad \frac{\partial \mathbf{U}}{\partial t}(\mathbf{x}, z, 0) = 0 \quad \text{on } \Omega. \qquad (0.5) $$

## Potential-based formulation

Problem (0.1)-(0.4) can be written as an abstract wave equation using an unbounded linear operator $G$,

$$ \frac{d^2 \mathbf{U}}{dt^2} + G^* \tilde{G} \mathbf{U} = 0, \qquad (0.6) $$

where $G^*$ is the adjoint of $G$ and $\tilde{G}$ is an extension of $G$. The introduction of this non-symmetric formulation using the extension $\tilde{G}$ is motivated by the presence of a non-homogeneous essential condition in the boundary conditions (0.4). The expression (0.6) will be obtained rigorously in Section 1. One originality of this work is to construct and analyze a “dual” or “adjoint” wave-like problem for a new unknown $\Phi$, satisfying

$$ \frac{d^2 \Phi}{dt^2} + \tilde{G} G^* \Phi = 0. \qquad (0.7) $$

Such problem is shown to be equivalent – in a sense given rigorously later – to the problem (0.6). Moreover, it presents several advantages from mathematical and numerical perspectives. The new unknown $\Phi$ has three scalar components, $\Phi = (\varphi, \psi, \gamma)^{\top}$, where $\varphi$ and $\psi$ are scalar fields and $\gamma = \varphi|_{\Gamma_s}$. Later, we show that $\Phi$ is related to the velocity $\mathbf{U}$ by the formula

$$ \mathbf{U} = G^* \Phi = -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z, \qquad (0.8) $$

hence it can be seen as a generalized potential. We also show that $\varphi$ and $\psi$ satisfy a second set of coupled partial differential equations describing the same physical system as (0.1)-(0.4), more precisely

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) + g \frac{\partial \varphi}{\partial z} - gN \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \qquad (0.9) $$

$$ \frac{\partial^2 \psi}{\partial t^2} - N \frac{\partial \varphi}{\partial z} + N^2 \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T]. \qquad (0.10) $$

The system is completed with boundary conditions. Using (0.4) and (0.8), we find on $\Gamma_b$

$$ \mathbf{U} \cdot \mathbf{n}_b = \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \cdot \mathbf{n}_b = u_b, \quad \text{on } \Gamma_b \times [0, T]. \qquad (0.11) $$

This condition is a generalized Neumann boundary conditions for the potentials $(\varphi, \psi)$. It is easier to consider than condition (0.4) for both the analysis and the discretization. Finally, to obtain the boundary condition on $\Gamma_s$ we observe that (0.9) can be rewritten, using (0.8),

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \mathbf{U} - g \mathbf{U} \cdot \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T]. $$

Formally, evaluating the equation above on the boundary $\Gamma_s$ and using (0.4) we obtain a boundary condition involving a second-order time derivative that accounts for surface gravity waves,

$$ \frac{\partial^2 \varphi}{\partial t^2} - g \mathbf{U} \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T], \qquad (0.12) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: The domain $\Omega$</figcaption></figure>

## 1.1 Hilbert spaces and trace operators

We start by introducing the space of $H^1$-functions,

$$H^1(\Omega) = \{\varphi \in L^2(\Omega) \mid \nabla\varphi \in L^2(\Omega)^d\},$$

as well as the usual surjective trace operator

$$\gamma_0 : H^1(\Omega) \to H^{1/2}(\Gamma),$$

that is the extension to function in $H^1(\Omega)$ of the trace operator $\varphi \mapsto \varphi|_{\Gamma}$, defined for smooth functions.

The forthcoming analysis requires the use of the standard space of square integrable functions with square integrable divergence,

$$H(\operatorname{div}, \Omega) = \{\mathbf{U} \in L^2(\Omega)^d \mid \nabla \cdot \mathbf{U} \in L^2(\Omega)\}.$$

The space $H(\operatorname{div}, \Omega)$ is equipped with the usual scalar product

$$\forall (\mathbf{U}, \tilde{\mathbf{U}}) \in H(\operatorname{div}, \Omega) \times H(\operatorname{div}, \Omega), \quad (\mathbf{U}, \tilde{\mathbf{U}})_{H(\operatorname{div}, \Omega)} = (\mathbf{U}, \tilde{\mathbf{U}})_{L^2(\Omega)^d} + (\nabla \cdot \mathbf{U}, \nabla \cdot \tilde{\mathbf{U}})_{L^2(\Omega)}.$$

In the following we use trace operators acting either on $\Gamma_s$ or on $\Gamma_b$ only, namely

$$\gamma_{0,s} : H^1(\Omega) \to H^{1/2}(\Gamma_s), \quad \gamma_{0,s}(\varphi) = \gamma_0(\varphi)|_{\Gamma_s},$$

$$\gamma_{0,b} : H^1(\Omega) \to H^{1/2}(\Gamma_b), \quad \gamma_{0,b}(\varphi) = \gamma_0(\varphi)|_{\Gamma_b}.$$

Since $\Gamma_s$ and $\Gamma_b$ are “well-separated” – i.e., the distance between the two boundary is at least $H - H_-$ that is strictly positive – one can also define surjective normal trace operators acting either on $\Gamma_s$ or on $\Gamma_b$ only. We introduce the normal trace operator

$$\gamma_{1,s} : H(\operatorname{div}, \Omega) \to H^{-1/2}(\Gamma_s).$$

When applied to smooth functions, it corresponds to the operator $\gamma_{1,s}(\mathbf{U}) = (\mathbf{U} \cdot \mathbf{n})|_{\Gamma_s}$. For functions in $H(\operatorname{div}, \Omega)$ the operator $\gamma_{1,s}$ is defined as follows. We introduce a function $\chi$ depending on $z$ only, and satisfying

$$\chi \in C^1([0, H]), \quad \chi(H) = 1, \quad \text{and} \quad \chi(z) = 0 \text{ for } z < H_-.$$

The normal trace operator on the surface is then defined by

$$\forall (\mathbf{U}, \varphi) \in H(\operatorname{div}, \Omega) \times H^1(\Omega), \quad \langle \gamma_{1,s}(\mathbf{U}), \gamma_{0,s}(\varphi) \rangle_{\Gamma_s} = (\nabla \cdot (\chi \mathbf{U}), \varphi)_{L^2(\Omega)} + (\chi \mathbf{U}, \nabla \varphi)_{L^2(\Omega)^d},$$

where the duality product between $H^{-1/2}(\Gamma_s)$ and $H^{1/2}(\Gamma_s)$ is denoted $\langle \cdot, \cdot \rangle_{\Gamma_s}$. In a similar way, we introduce $\gamma_{1,b}$, the normal trace operator acting on $\Gamma_b$ only. We denote by $\langle \cdot, \cdot \rangle_{\Gamma_b}$ the duality product between $H^{-1/2}(\Gamma_b)$ and $H^{1/2}(\Gamma_b)$ and define the normal trace for functions in $H(\operatorname{div}, \Omega)$ by

$$\forall (\mathbf{U}, \varphi) \in H(\operatorname{div}, \Omega) \times H^1(\Omega), \quad \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_s} = (\nabla \cdot ((1-\chi)\mathbf{U}), \varphi)_{L^2(\Omega)} + ((1-\chi)\mathbf{U}, \nabla \varphi)_{L^2(\Omega)^d}.$$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. The domain $\Omega$</figcaption></figure>

where $\mathbf{n}_s$ is the outward unitary normal of the domain $\Omega$ on $\Gamma_s$. The initial conditions are deduced from (0.5), they read

$$ \varphi(\mathbf{x}, z, 0) = \frac{\partial \varphi}{\partial t}(\mathbf{x}, z, 0) = 0, \quad \psi(\mathbf{x}, z, 0) = \frac{\partial \psi}{\partial t}(\mathbf{x}, z, 0) = 0 \quad \text{on } \Omega. \qquad (0.13) $$

If $d = 3$, the system (0.9) - (0.13) is solved numerically for two scalar fields ($\varphi$ and $\psi$), compared to three for the velocity based formulation (0.1). Moreover, this system extends common equations for hydrodynamics [15] to the non-barotropic case. In the barotropic case, we have $N = 0$ and the system of partial differential equations (0.9)-(0.10) reduces to

$$ \frac{\partial^2 \varphi}{\partial t^2} - c_0^2 \Delta \varphi + g \frac{\partial \varphi}{\partial z} = 0, \quad \text{in } \Omega \times [0, T], \qquad (0.14) $$

$$ \frac{\partial^2 \psi}{\partial t^2} = 0, \quad \text{in } \Omega \times [0, T]. \qquad (0.15) $$

Because of the vanishing initial conditions, we have $\psi = 0$. The system is then described by the function $\varphi$ alone. The boundary conditions are also simplified,

$$ \frac{\partial^2 \varphi}{\partial t^2} + g \nabla \varphi \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T] \quad \text{and} \quad \nabla \varphi \cdot \mathbf{n}_b = u_b, \quad \text{on } \Gamma_b \times [0, T]. \qquad (0.16) $$

The system (0.14)-(0.16) was previously introduced in the literature [11, 25].

**Remark 0.1.** If non vanishing initial data were considered, typically,

$$ \mathbf{U} = \mathbf{U}_0, \quad \frac{\partial \mathbf{U}}{\partial t} = \mathbf{U}_1 \quad \text{and} \quad \varphi = \varphi_0, \quad \frac{\partial \varphi}{\partial t} = \varphi_1, \quad \psi = \psi_0, \quad \frac{\partial \psi}{\partial t} = \psi_1, \quad \text{on } \Omega \times \{0\} $$

then some compatibility conditions should be satisfied by $\mathbf{U}_0$ and $\mathbf{U}_1$ in order to show equivalence between the two formulations. These compatibility conditions are naturally deduced from (0.8) and reads

$$ \mathbf{U}_0 = -\nabla\varphi_0 + N \left(\psi_0 + \frac{N}{g}\varphi_0\right) \mathbf{e}_z \quad \text{and} \quad \mathbf{U}_1 = -\nabla\varphi_1 + N \left(\psi_1 + \frac{N}{g}\varphi_1\right) \mathbf{e}_z. $$

# 1. PRELIMINARY DEFINITIONS

In this section, we define functional spaces and operators that will be used throughout the paper.
Finally, for any function $\mathbf{U} \in H(\mathrm{div}, \Omega)$, when stating that $\gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s)$ we mean that $\gamma_{1,s}(\mathbf{U})$ is a function in $H^{-1/2}(\Gamma_s)$ which can be identified with a function in $L^2(\Gamma_s)$, and that the duality product reduces to the scalar product in $L^2(\Gamma_s)$. This can be written as follows,

$$ \forall \mathbf{U} \in H(\mathrm{div}, \Omega), \quad \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s) \Rightarrow \exists f \in L^2(\Gamma_s) / \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_s} = \int_{\Gamma_s} f \gamma_{0,b}(\varphi) \, \mathrm{d}s. $$

When $\mathbf{U}$ is smooth, the function $f$ is given by $f = (\mathbf{U} \cdot \mathbf{n})|_{\Gamma_s}$.

## 1.2 The operator $\tilde{G}$ and its extension $\tilde{G}$

To introduce the operator associated to the evolution problem (0.1) and the abstract wave equation (0.6), we introduce the Hilbert space $\mathcal{H} = L^2(\Omega)^d$ equipped with the following weighted scalar product,

$$ (\mathbf{U}, \tilde{\mathbf{U}})_{\mathcal{H}} = \int_{\Omega} \rho_0 \mathbf{U} \cdot \tilde{\mathbf{U}} \, \mathrm{d}s. \quad (1.1) $$

we also define the space $\mathcal{G}$,

$$ \mathcal{G} = L^2(\Omega) \times L^2(\Omega) \times L^2(\Gamma_s), \quad (1.2) $$

equipped with the weighted scalar product

$$ \forall \Phi = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix} \in \mathcal{G}, \quad \forall \tilde{\Phi} = \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \\ \tilde{\gamma} \end{pmatrix} \in \mathcal{G}, \quad (\Phi, \tilde{\Phi})_{\mathcal{G}} = \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, \mathrm{d}s + \int_{\Omega} \rho_0 \psi \tilde{\psi} \, \mathrm{d}s + \int_{\Gamma_s} \frac{\rho_0}{g} \gamma \tilde{\gamma} \, \mathrm{d}s. \quad (1.3) $$

We first define the operator $\tilde{G}$ used in Equation (0.6). The domain of $\tilde{G}$ is denoted $\mathcal{D}(\tilde{G}) \subset \mathcal{H}$ and is defined by

$$ \mathcal{D}(\tilde{G}) = \{ \mathbf{U} \in H(\mathrm{div}, \Omega) \mid \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s) \}. \quad (1.4) $$

The operator $\tilde{G} : \mathcal{D}(\tilde{G}) \subset \mathcal{H} \to \mathcal{G}$ is then defined by

$$ \forall U \in D(\tilde{G}), \quad \tilde{G} \mathbf{U} = \begin{pmatrix} c_0^2 \left( \nabla \cdot \mathbf{U} - \frac{g}{c_0^2} \mathbf{U} \cdot \mathbf{e}_z \right) \\ N \mathbf{U} \cdot \mathbf{e}_z \\ -g \gamma_{1,s}(\mathbf{U}) \end{pmatrix}. \quad (1.5) $$

It can be shown that the operator $\tilde{G}$ is closed and densely defined. As already mentioned, it is useful to see the operator $\tilde{G}$ as an extension of an operator $G$, defined on the domain $\mathcal{D}(G) \subset \mathcal{H}$, given by

$$ \mathcal{D}(G) = \{ \mathbf{U} \in H(\mathrm{div}, \Omega) \mid \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s), \gamma_{1,b}(\mathbf{U}) = 0 \}. \quad (1.6) $$

We have $\mathcal{D}(G) \subset \mathcal{D}(\tilde{G})$ and by definition $G$ satisfies, for all $\mathbf{U} \in \mathcal{D}(G)$, $G\mathbf{U} = \tilde{G}\mathbf{U}$, therefore

$$ \forall U \in D(G), \quad G \mathbf{U} = \begin{pmatrix} c_0^2 \left( \nabla \cdot \mathbf{U} - \frac{g}{c_0^2} \mathbf{U} \cdot \mathbf{e}_z \right) \\ N \mathbf{U} \cdot \mathbf{e}_z \\ -g \gamma_{1,s}(\mathbf{U}) \end{pmatrix}. \quad (1.7) $$

The operator $G$ is also closed and densely defined.

## 1.3 The adjoint operators $G^*$ and $\tilde{G}^*$ and a Green's formula

Since the operators $G$ and $\tilde{G}$ are densely defined and closed, their adjoint – denoted $G^*$ and $\tilde{G}^*$ respectively – exist and are also densely defined and closed. We give their expression in this section. In the following, the space of smooth functions with compact support in $\Omega$ is denoted $\mathcal{D}(\Omega)$ and $\mathcal{D}(\tilde{\Omega})$ corresponds to the space of smooth functions up to the boundary.
1.1. Hilbert spaces and trace operators

We start by introducing the space of H¹-functions,

$$
H^1(\Omega) = \{\varphi \in L^2(\Omega) \mid \nabla\varphi \in L^2(\Omega)^d\},
$$

as well as the usual surjective trace operator

$$
\gamma_0 : H^1(\Omega) \to H^{1/2}(\Gamma),
$$

that is the extension to function in H¹(Ω) of the trace operator ϕ ↦→ ϕ|Γ, defined for smooth functions.

The forthcoming analysis requires the use of the standard space of square integrable functions with square integrable divergence,

$$
H(\mathrm{div}, \Omega) = \{ \mathbf{U} \in L^2(\Omega)^d \mid \nabla \cdot \mathbf{U} \in L^2(\Omega) \}.
$$

In the following we use continuous trace operators acting either on Γs or on Γb only, namely

$$
\gamma_{0,s} : H^1(\Omega) \to H^{1/2}(\Gamma_s), \quad \gamma_{0,s}(\varphi) = \gamma_0(\varphi)|_{\Gamma_s},
$$

$$
\gamma_{0,b} : H^1(\Omega) \to H^{1/2}(\Gamma_b), \quad \gamma_{0,b}(\varphi) = \gamma_0(\varphi)|_{\Gamma_b}.
$$

Since Γs and Γb are “well-separated” – i.e., the distance between the two boundaries is at least H − H − that is strictly positive – one can also define continuous surjective normal trace operators acting either on Γs or on Γb only. We introduce the normal trace operator

$$
\gamma_{1,s} : H(\mathrm{div}, \Omega) \to H^{-1/2}(\Gamma_s).
$$

When applied to smooth functions, it corresponds to the operator γ₁,s(U) = (U · n)|Γs. For functions in H(div, Ω) the operator γ₁,s is defined as follows. We introduce a function χ depending on z only, and satisfying

$$
\chi \in C^1([0, H]), \quad \chi(H) = 1, \quad \text{and} \quad \chi(z) = 0 \text{ for } z < H_-
$$

The normal trace operator on the surface is then defined by

$$
\forall (\mathbf{U}, \varphi) \in H(\mathrm{div}, \Omega) \times H^1(\Omega), \quad \langle \gamma_{1,s}(\mathbf{U}), \gamma_{0,s}(\varphi) \rangle_{\Gamma_s} = (\nabla \cdot (\chi(\mathbf{U}), \varphi))_{L^2(\Omega)} + (\chi(\mathbf{U}), \nabla \varphi)_{L^2(\Omega)^d},
$$

where 〈·,·〉Γs denotes the duality product between H−1/2(Γs) and H1/2(Γs). In a similar way, we introduce
γ1,b, the normal trace operator acting on Γb. We denote by 〈·,·〉Γb the duality product between H−1/2(Γb) and
H1/2(Γb) and define the normal trace for functions in H(div, Ω) by

$$
\forall (\mathbf{U}, \varphi) \in H(\mathrm{div}, \Omega) \times H^1(\Omega), \quad \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_s} = (\nabla \cdot ((1-\chi)\mathbf{U}), \varphi)_{L^2(\Omega)} + ((1-\chi)\mathbf{U}, \nabla\varphi)_{L^2(\Omega)^d}.
$$

Finally, when stating that $\gamma_{1,s}(\mathbf{U})$ belongs to $L^2(\Gamma_s)$, we mean that $\gamma_{1,s}(\mathbf{U})$ is a function in $H^{-1/2}(\Gamma_s)$ which can be identified with a function in $L^2(\Gamma_s)$, and that the duality product reduces to the scalar product in $L^2(\Gamma_s)$.
This can be written as follows:

$$
\forall \mathbf{U} \in H(\mathrm{div}, \Omega), \quad \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s) \Rightarrow \exists f \in L^2(\Gamma_s) / \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_s} = \int_{\Gamma_s} f \gamma_{0,b}(\varphi) \, ds.
$$

When **U** is smooth, the function *f* is given by *f* = (**U** · **n**)∣<sub>Γs</sub>.
**Theorem 1.1.** The operator $G^*: \mathcal{D}(G^*) \subset G \to \mathcal{H}$ is defined by

$$ \mathcal{D}(G^*) = \{\Phi = (\varphi, \psi, \gamma)^t \in G \mid \varphi \in H^1(\Omega), \gamma = \gamma_{0,s}(\varphi)\}, $$

and, for all $\Phi = (\varphi, \psi, \gamma)^t \in \mathcal{D}(G^*)$,

$$ G^*\Phi = -\nabla\varphi + N\left(\psi + \frac{N}{g}\varphi\right)e_z. \qquad (1.8) $$

*Proof.* Let $\mathbf{U}$ be a function in $\mathcal{D}(\Omega)^d$, and let $\Phi = (\varphi \ \psi \ \gamma)^t$ belong to $\mathcal{D}(G^*)$. It holds, by definition of the adjoint,

$$ (G\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_\mathcal{H} = (\mathbf{U}, \tilde{\mathbf{U}})_\mathcal{H}, $$

for some $\tilde{\mathbf{U}} \in \mathcal{H}$. The equality above is developed using the definition of $G$,

$$ (G\mathbf{U}, \Phi)_G = -\langle\nabla(\rho_0\varphi), \mathbf{U}\rangle_\Omega + \int_\Omega \rho_0 \left(N\psi - \frac{g}{c_0^2}\varphi\right) \mathbf{U} \cdot e_z \,dx = (\mathbf{U}, \tilde{\mathbf{U}})_\mathcal{H}, \qquad (1.9) $$

where $\langle \cdot, \cdot \rangle_\Omega$ correspond to the duality product in $\mathcal{D}(\Omega)^d$. The equality (1.9) shows that $\nabla(\rho_0\varphi)$ belongs to $L^2(\Omega)^3$, hence, since $\rho_0$ is smooth, $\varphi \in H^1(\Omega)$. Equation (1.9) also shows that

$$ G^*\Phi = \tilde{\mathbf{U}} = -\rho_0^{-1}\nabla(\rho_0\varphi) + \left(N\psi - \frac{g}{c_0^2}\varphi\right)e_z. $$

The simpler expression (1.8) is obtained by, first, distributing the gradient, then, using the definition of the scalar field $N$. $\square$

The following Green formula holds

**Lemma 1.2.** For all $\mathbf{U} \in \mathcal{D}(\tilde{\mathbf{G}})$ and $\Phi = (\varphi, \psi, \gamma)^t \in \mathcal{D}(G^*)$, it holds

$$ (\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_\mathcal{H} + \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \qquad (1.10) $$

*Proof.* For $\mathbf{U} \in \mathcal{D}(\bar{\Omega})^d$ and $\Phi = (\varphi, \psi, \varphi|_{\Gamma_s})^t$ with $\varphi$ and $\psi$ in $\mathcal{D}(\bar{\Omega})$, we have

$$ (\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_\mathcal{H} + \int_{\Gamma_b} \varphi|_{\Gamma_b} (\mathbf{U} \cdot \mathbf{n}_b)|_{\Gamma_b} \,ds = \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. $$

We conclude by using the density of $\mathcal{D}(\bar{\Omega})^d$ in $H(\text{div}, \Omega)$ and the density of $\mathcal{D}(\bar{\Omega})$ in $H^1(\Omega)$ and in $L^2(\Omega)$ (see [16] for details on these standard density results). $\square$

Thanks to the Lemma 1.2 we deduce the expression of $\tilde{\mathbf{G}}^*$.

**Corollary 1.3.** The operator $\tilde{\mathbf{G}}^*: \mathcal{D}(\tilde{\mathbf{G}}^*) \subset G \to \mathcal{H}$ is defined by

$$ \mathcal{D}(\tilde{\mathbf{G}}^*) = \{\Phi = (\varphi, \psi, \gamma)^t \in G \mid \varphi \in H^1(\Omega), \gamma_{0,s}(\varphi) = \gamma, \gamma_{0,b}(\varphi) = 0\} \subset \mathcal{D}(G^*), \qquad (1.11) $$

and for all $\Phi \in \mathcal{D}(\tilde{\mathbf{G}}^*)$, $\tilde{\mathbf{G}}^*\Phi = G^*\Phi$.

*Proof.* Since $\mathcal{D}(G) \subset \mathcal{D}(\tilde{\mathbf{G}})$, we have the inclusion $\mathcal{D}(\tilde{\mathbf{G}}^*) \subset \mathcal{D}(G^*)$ and $G^*$ is an extension of $\tilde{\mathbf{G}}^*$. Therefore Lemma 1.2 can be used as follows: for all $\mathbf{U} \in \mathcal{D}(\tilde{\mathbf{G}})$ and $\Phi = (\varphi, \psi, \gamma)^t \in \mathcal{D}(\tilde{\mathbf{G}}^*)$, it holds

$$ (\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G - (\mathbf{U}, \tilde{\mathbf{G}}^*\Phi)_\mathcal{H} = \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}, \qquad (1.12) $$

which shows that $\langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b} = 0$. Using the surjectivity of the normal trace operator $\gamma_{1,b}$, we deduce that $\gamma_{0,b}(\varphi) = 0$. $\square$
## 1.2. The operator $G$ and its extension $\tilde{G}$

To introduce the operator associated to the evolution problem (0.1) and the abstract wave equation (0.6), we introduce the Hilbert space $\mathcal{H} = L^2(\Omega)^d$ equipped with the weighted scalar product

$$ (\mathbf{U}, \tilde{\mathbf{U}})_{\mathcal{H}} = \int_{\Omega} \rho_0 \mathbf{U} \cdot \tilde{\mathbf{U}} \, dx. \qquad (1.1) $$

We also define the space $\mathcal{G}$,

$$ \mathcal{G} = L^2(\Omega) \times L^2(\Omega) \times L^2(\Gamma_s), \qquad (1.2) $$

equipped with the weighted scalar product

$$ \forall \mathbf{\Phi} = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix} \in \mathcal{G}, \quad \forall \tilde{\mathbf{\Phi}} = \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \\ \tilde{\gamma} \end{pmatrix} \in \mathcal{G}, \quad (\mathbf{\Phi}, \tilde{\mathbf{\Phi}})_{\mathcal{H}} = \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, dx + \int_{\Omega} \rho_0 \psi \tilde{\psi} \, dx + \int_{\Gamma_s} \frac{\rho_0}{g} \gamma \tilde{\gamma} \, ds. \qquad (1.3) $$

We first define the operator $\tilde{G}$ used in Equation (0.6). The domain of $\tilde{G}$ denoted by $\mathcal{D}(\tilde{G}) \subset \mathcal{H}$ is defined by

$$ \mathcal{D}(\tilde{G}) = \{ \mathbf{U} \in H(\mathrm{div}, \Omega) \mid \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s) \}. \qquad (1.4) $$

The operator $\tilde{G} : \mathcal{D}(\tilde{G}) \subset \mathcal{H} \to \mathcal{G}$ is then defined by

$$ \forall U \in \mathcal{D}(\tilde{G}), \quad \tilde{G} \mathbf{U} = \begin{pmatrix} c_0^2 \left( \nabla \cdot \mathbf{U} - \frac{g}{c_0^2} \mathbf{U} \cdot \mathbf{e}_z \right) \\ N \mathbf{U} \cdot \mathbf{e}_z \\ -g \gamma_{1,s}(\mathbf{U}) \end{pmatrix}. \qquad (1.5) $$

It can be shown that the operator $\tilde{G}$ is closed and densely defined (see Appendix 7.1.1). As already mentioned, it is useful to see the operator $\tilde{G}$ as an extension of an operator $G$, defined on the domain $\mathcal{D}(G) \subset \mathcal{H}$, given by

$$ \mathcal{D}(G) = \{ \mathbf{U} \in H(\mathrm{div}, \Omega) \mid \gamma_{1,s}(\mathbf{U}) \in L^2(\Gamma_s), \gamma_{1,b}(\mathbf{U}) = 0 \}. \qquad (1.6) $$

We have $\mathcal{D}(G) \subset \mathcal{D}(\tilde{G})$ and by definition $G$ satisfies, for all $\mathbf{U} \in \mathcal{D}(G)$, $G\mathbf{U} = \tilde{G}\mathbf{U}$, therefore

$$ \forall U \in \mathcal{D}(G), \quad G U = \begin{pmatrix} c_0^2 \left( \nabla \cdot \mathbf{U} - \frac{g}{c_0^2} \mathbf{U} \cdot \mathbf{e}_z \right) \\ N \mathbf{U} \cdot \mathbf{e}_z \\ -g \gamma_{1,s}(\mathbf{U}) \end{pmatrix}. \qquad (1.7) $$

The operator $G$ is also closed and densely defined (see Appendix 7.1.1).

## 1.3. The adjoint operators $G^*$ and $\tilde{G}^*$ and Green's formula

Since the operators $G$ and $\tilde{G}$ are densely defined and closed, their adjoint – denoted by $G^*$ and $\tilde{G}^*$ respectively – exist and are also densely defined and closed. We give their expression in this section. In the following, we denote by $\mathcal{D}(\Omega)$ the space of smooth functions with compact support in $\Omega$, and let $\mathcal{D}(\bar{\Omega})$ be the space of smooth functions up to the boundary.

**Theorem 1.1.** The operator $G^* : \mathcal{D}(G^*) \subset \mathcal{G} \to \mathcal{H}$ is defined by

$$ \mathcal{D}(G^*) = \{ \Phi = (\varphi, \psi, \gamma)^{\top} \in \mathcal{G} \mid \varphi \in H^1(\Omega), \gamma = \gamma_{0,s}(\varphi) \}, $$
The space $\mathcal{D}(G^*)$ is equipped with the graph norm,

$$ \|\Phi\|_{\mathcal{D}(G^*)}^2 = \|\Phi\|_{\mathcal{G}}^2 + \|G^*\Phi\|_{\mathcal{H}}^2. $$

And we have the following result,

**Proposition 1.4.** There exists a constant $C_c > 0$ such that

$$ \forall \Phi = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix} \in \mathcal{D}(G^*), \quad \|\Phi\|_{\mathcal{D}(G^*)} \ge C_c \|\varphi\|_{H^1(\Omega)}. $$

*Proof.* In the proof, we use the symbol $\lesssim$ for inequalities that hold up to a constant independent of $\Phi$. For $\Phi \in \mathcal{D}(G^*)$, we have

$$ \|\varphi\|_{L^2(\Omega)^3} \lesssim \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \, dx \le \|\Phi\|_{\mathcal{G}}^2 \quad \text{and} \quad \|\nabla \varphi\|_{L^2(\Omega)^d} \lesssim \|\nabla \varphi\|_{\mathcal{H}}^2. $$

It holds, by the triangular inequality,

$$ \|\nabla \varphi\|_{\mathcal{H}}^2 \lesssim \left\| -\nabla \varphi + N\left(\psi + \frac{N}{g}\varphi\right)\mathbf{e}_z \right\|_{\mathcal{H}}^2 + \left\| N\left(\psi + \frac{N}{g}\varphi\right)\mathbf{e}_z \right\|_{\mathcal{H}}^2, $$

hence the norm of the gradient is bounded by $\|G^*\Phi\|_{\mathcal{H}}^2 + \|\Phi\|_{\mathcal{G}}^2$, which concludes the proof. $\square$

In this section, we have introduced all the necessary operators for the study of the potential-based problem (0.9)-(0.13) and the velocity-field problem (0.1)-(0.5).

## 2 Analysis of the potential-based formulation

This section is dedicated to the study of the problem (0.9)-(0.13). We show that the problem is well posed and that its solution satisfies an energy equality. We recall the system of PDE satisfied by the potentials $\varphi$ and $\psi$,

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) + g \frac{\partial \varphi}{\partial z} - gN \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \quad (2.1) $$

$$ \frac{\partial^2 \psi}{\partial t^2} - N \frac{\partial \varphi}{\partial z} + N^2 \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \quad (2.2) $$

with boundary conditions

$$ -\nabla \varphi \cdot \mathbf{n}_b + N \left( \psi + \frac{N}{g} \varphi \right) (\mathbf{e}_z \cdot \mathbf{n}_b) = u_b, \quad \text{on } \Gamma_b \times [0, T], \quad (2.3) $$

and

$$ \frac{\partial^2 \varphi}{\partial t^2} - g \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T]. \quad (2.4) $$

The system (2.1)-(2.4) is completed with vanishing initial conditions.

### 2.1 Variational formulation

The natural idea for writing the variational formulation associated to (2.1)-(2.2) consists in testing (2.1)-(2.2) against a function $(\tilde{\varphi}, \tilde{\psi}) \in H^1(\Omega) \times L^2(\Omega)$. After integrating by parts, using the boundary conditions (2.3)-(2.4) and the definition of $N$ given by (0.3), we obtain the problem: given $u_b$ regular enough, find

$$ \begin{pmatrix} \varphi \\ \psi \end{pmatrix} \in L^2(0, T; H^1(\Omega) \times L^2(\Omega)), \quad \frac{d}{dt} \begin{pmatrix} \varphi \\ \psi \end{pmatrix} \in L^2(0, T; L^2(\Omega)^2), \quad (2.5) $$
and, for all $\Phi = (\varphi, \psi, \gamma)^{\top} \in \mathcal{D}(G^*)$,

$$G^*\Phi = -\nabla\varphi + N\left(\psi + \frac{N}{g}\varphi\right)\mathbf{e}_z. \qquad (1.8)$$

*Proof.* Let $\mathbf{U}$ be a function in $\mathcal{D}(\Omega)^d$, and let $\Phi = (\varphi, \psi, \gamma)^{\top}$ belong to $\mathcal{D}(G^*)$. By definition of the adjoint, we have

$$(G\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_H = (\mathbf{U}, \tilde{\mathbf{U}})_H,$$

for some $\tilde{\mathbf{U}} \in \mathcal{H}$. The equality above is developed using the definition of $G$,

$$(G\mathbf{U}, \Phi)_G = -\langle \nabla(\rho_0\varphi), \mathbf{U} \rangle_{\Omega} + \int_{\Omega} \rho_0 \left(N\psi - \frac{g}{c_0^2}\varphi\right) \mathbf{U} \cdot \mathbf{e}_z \, dx = (\mathbf{U}, \tilde{\mathbf{U}})_H, \qquad (1.9)$$

where $\langle \cdot, \cdot \rangle_{\Omega}$ correspond to the duality product in $\mathcal{D}(\Omega)^d$. The equality (1.9) shows that $\nabla(\rho_0\varphi)$ belongs to $L^2(\Omega)^3$, hence, since $\rho_0$ is smooth, $\varphi$ belongs to $H^1(\Omega)$. Equation (1.9) also shows that

$$G^*\Phi = \tilde{\mathbf{U}} = -\rho_0^{-1}\nabla(\rho_0\varphi) + \left(N\psi - \frac{g}{c_0^2}\varphi\right)\mathbf{e}_z.$$

The expression (1.8) is obtained by distributing the gradient and using the definition of the scalar field $N$. $\square$

The following Green formula holds:

**Lemma 1.2.** For all $\mathbf{U} \in \mathcal{D}(\tilde{\mathbf{G}})$ and $\Phi = (\varphi, \psi, \gamma)^{\top} \in \mathcal{D}(G^*)$, we have

$$(\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_H + \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \qquad (1.10)$$

*Proof.* For $\mathbf{U} \in \mathcal{D}(\bar{\Omega})^d$ and $\Phi = (\varphi, \psi, \varphi|_{\Gamma_s})^{\top}$ with $\varphi$ and $\psi$ in $\mathcal{D}(\bar{\Omega})$, we have

$$(\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G = (\mathbf{U}, G^*\Phi)_H + \int_{\Gamma_b} \varphi|_{\Gamma_b} (\mathbf{U} \cdot \mathbf{n}_b)|_{\Gamma_b} \, ds = \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}.$$

We conclude by using the density of $\mathcal{D}(\bar{\Omega})^d$ in $H(\text{div}, \Omega)$ and the density of $\mathcal{D}(\bar{\Omega})$ in $H^1(\Omega)$ and in $L^2(\Omega)$ (see [16] for details on these standard density results). $\square$

Thanks to the Lemma 1.2 we deduce the expression of $\tilde{G}^*$.

**Corollary 1.3.** The operator $\tilde{G}^*$: $\mathcal{D}(\tilde{G}^*) \subset \mathcal{G} \to \mathcal{H}$ is defined by

$$\mathcal{D}(\tilde{G}^*) = \{\Phi = (\varphi, \psi, \gamma)^{\top} \in \mathcal{G} \mid \varphi \in H^1(\Omega), \gamma_{0,s}(\varphi) = \gamma, \gamma_{0,b}(\varphi) = 0\} \subset \mathcal{D}(G^*), \qquad (1.11)$$

and for all $\Phi \in \mathcal{D}(\tilde{G}^*)$, $\tilde{G}^*\Phi = G^*\Phi$.

*Proof.* Since $\mathcal{D}(G) \subset \mathcal{D}(\tilde{\mathbf{G}})$, we have the inclusion $\mathcal{D}(\tilde{G}^*) \subset \mathcal{D}(G^*)$ and $G^*$ is an extension of $\tilde{G}^*$. Therefore Lemma 1.2 can be used as follows: for all $\mathbf{U} \in \mathcal{D}(\tilde{\mathbf{G}})$ and $\Phi = (\varphi, \psi, \gamma)^{\top} \in \mathcal{D}(\tilde{G}^*)$, we have

$$(\tilde{\mathbf{G}}\mathbf{U}, \Phi)_G - (\mathbf{U}, \tilde{G}^*\Phi)_H = \langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}, \qquad (1.12)$$

which shows that $\langle \gamma_{1,b}(\mathbf{U}), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b} = 0$. Using the surjectivity of the normal trace operator $\gamma_{1,b}$, we deduce that $\gamma_{0,b}(\varphi) = 0$. $\square$
solution to

$$ \frac{d^2}{dt^2} \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, dx + \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \psi \tilde{\psi} \, dx \\ + \int_{\Omega} \rho_0 \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \left( -\nabla \tilde{\varphi} + N \left( \tilde{\psi} + \frac{N}{g} \tilde{\varphi} \right) \mathbf{e}_z \right) \, dx \\ + \frac{d^2}{dt^2} \int_{\Gamma_s} \frac{\rho_0}{g} \varphi \tilde{\varphi} \, ds + \int_{\Gamma_b} \rho_0 u_b \gamma_{0,b}(\tilde{\varphi}) \, ds = 0, \quad \forall \left( \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \end{pmatrix} \right) \in H^1(\Omega) \times L^2(\Omega). \quad (2.6) $$

The formulation (2.6) will be useful for the numerical approximation. Indeed, the natural spaces for the discretization are classical: $H^1(\Omega)$ for $\varphi$ and $L^2(\Omega)$ for $\psi$. Moreover, for $d=3$, the velocity formulation has three scalar unknowns whereas the potential formulation requires only two scalars unknowns. Finally, the source $u_b$ appears naturally as a Neumann condition in this formulation.

## 2.2 Existence and uniqueness results

The existence of a solution to (2.1)-(2.2) cannot be directly proved by standard methods (such as [24]), because of the surface condition (2.4) involving the second-order time derivative of $\varphi$. Therefore, we need to introduce a new unknown to the problem, denoted $\gamma \in L^2(\Gamma_s)$, and we define the vector of unknowns

$$ \Phi(t) = \begin{pmatrix} \varphi(t) \\ \psi(t) \\ \gamma(t) \end{pmatrix} \in \mathcal{G}. \quad (2.7) $$

From the variational formulation (2.6) we deduce that the problem (2.1) - (2.4) reduces to the following abstract formulation: assume $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$ given, find

$$ \Phi \in L^2(0, T; \mathcal{D}(G^*)), \quad \frac{d}{dt} \Phi \in L^2(0, T; \mathcal{G}), $$

solution to

$$ \frac{d^2}{dt^2}(\Phi(t), \tilde{\Phi})_{\mathcal{G}} + (G^*\Phi(t), G^*\tilde{\Phi})_{\mathcal{H}} = \ell_b(t, \tilde{\Phi}), \quad \forall \tilde{\Phi} \in \mathcal{D}(G^*), \quad \text{in } \mathcal{D}'(0, T), \quad (2.8) $$

$$ \Phi(0) = \frac{d}{dt} \Phi(0) = 0, \quad (2.9) $$

where $\ell_b : (0, T) \times \mathcal{D}(G^*) \to \mathbb{R}$ is the linear form

$$ \ell_b(t, \Phi) = \langle u_b(t), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \quad (2.10) $$

Note that if $\Phi \in \mathcal{D}(G^*)$, then $\gamma$ is the surface trace of $\varphi$, and the equation (2.8) is exactly the equation (2.6). For the formulation (2.8)-(2.9) we have the following result,

**Proposition 2.1.** Assume that $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$. Then the problem (2.8)-(2.9) has a unique solution and, up to a modification on zero measure sets,

$$ \Phi \in C^0([0, T]; \mathcal{D}(G^*)) \cap C^1([0, T]; \mathcal{G}). $$

*Proof.* First we show that, for almost all $t \in (0, T)$, the form $\ell_b(t)$ is a bounded linear functional on $\mathcal{D}(G^*)$. Let $\Phi = (\varphi, \psi, \gamma_{0,s}(\varphi))^t \in \mathcal{D}(G^*)$. From the continuity of the trace operator, there exists a scalar $C_H > 0$ depending only on $H$, such that

$$ |\langle \ell_b(t), \Phi \rangle| = |\langle u_b(t), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}| \le C_H \|u_b(t)\|_{H^{-1/2}(\Gamma_b)} \|\varphi\|_{H^1(\Omega)}, $$

and from Proposition 1.4 we obtain

$$ |\langle \ell_b(t), \Phi \rangle| \le C_H C_c^{-1} \|u_b(t)\|_{H^{-1/2}(\Gamma_b)} \|\Phi\|_{\mathcal{D}(G^*)}, $$
The space $\mathcal{D}(G^*)$ is equipped with the graph norm,

$$ \|\Phi\|_{\mathcal{D}(G^*)}^2 = \|\Phi\|_{\mathcal{G}}^2 + \|G^*\Phi\|_{\mathcal{H}}^2. $$

And we have the following result:

**Proposition 1.4.** There exists a constant $C_c > 0$ such that

$$ \forall \Phi = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix} \in \mathcal{D}(G^*), \quad \|\Phi\|_{\mathcal{D}(G^*)} \ge C_c \|\varphi\|_{H^1(\Omega)}. $$

*Proof.* We use the symbol $\lesssim$ for inequalities that hold up to a multiplicative constant independent of $\Phi$. For $\Phi \in \mathcal{D}(G^*)$, we have

$$ \|\varphi\|_{L^2(\Omega)^3} \lesssim \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi dx \le \|\Phi\|_{\mathcal{G}}^2 \quad \text{and} \quad \|\nabla \varphi\|_{L^2(\Omega)^d} \lesssim \|\nabla \varphi\|_{\mathcal{H}}^2. $$

Using the triangular inequality, we get

$$ \|\nabla \varphi\|_{\mathcal{H}}^2 \lesssim \left\| -\nabla \varphi + N\left(\psi + \frac{N}{g}\varphi\right)\mathbf{e}_z \right\|_{\mathcal{H}}^2 + \left\| N\left(\psi + \frac{N}{g}\varphi\right)\mathbf{e}_z \right\|_{\mathcal{H}}^2, $$

hence the norm of the gradient is bounded by $\|G^*\Phi\|_{\mathcal{H}}^2 + \|\Phi\|_{\mathcal{G}}^2$, which concludes the proof. $\square$

In this section, we have introduced all the necessary operators for the study of the potential-based problem (0.9)-(0.13) and the velocity-field problem (0.1)-(0.5).

## 2. ANALYSIS OF THE POTENTIAL-BASED FORMULATION

This section is dedicated to the study of the problem (0.9)-(0.13). We show that the problem is well-posed and that its solution satisfies an energy equality. We recall the system satisfied by the potentials $\varphi$ and $\psi$

$$ \frac{\partial^2 \varphi}{\partial t^2} + c_0^2 \nabla \cdot \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) + g \frac{\partial \varphi}{\partial z} - gN \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \quad (2.1) $$

$$ \frac{\partial^2 \psi}{\partial t^2} - N \frac{\partial \varphi}{\partial z} + N^2 \left( \psi + \frac{N}{g} \varphi \right) = 0, \quad \text{in } \Omega \times [0, T], \quad (2.2) $$

with boundary conditions

$$ \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \cdot \mathbf{n}_b = u_b, \quad \text{on } \Gamma_b \times [0, T], \quad (2.3) $$

$$ \frac{\partial^2 \varphi}{\partial t^2} - g \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \cdot \mathbf{n}_s = 0, \quad \text{on } \Gamma_s \times [0, T]. \quad (2.4) $$

The system (2.1)-(2.4) is completed with vanishing initial conditions.
hence $l_b(t)$ is bounded. From the continuity in time of $u_b$, we have $l_b \in H^1(0, T; D(G^*)')$. From this property and (2.8), we deduce that

$$ \frac{d^2}{dt^2} \Phi \in L^2(0, T; D(G^*)'), $$

hence the initial conditions (2.9) make sense. The existence and uniqueness of a solution to the problem (2.8)-(2.9) follows then directly from standard results, see e.g. [24] and [21]. □

## 2.3 Energy identity

An energy identity for the system (2.8)-(2.9) is obtained following the usual approach [24]. We define the energy

$$ \mathcal{E}(t) = \frac{1}{2} (\|\partial_t \Phi\|_{G}^{2} + \|G^{*} \Phi\|_{H}^{2}). \qquad (2.11) $$

Taking formally $\tilde{\Phi} = \partial_t \Phi$ in the weak formulation (2.8) yields

$$ (\partial_{tt}^2 \Phi, \partial_t \Phi)_G + (G^* \Phi, G^* \partial_t \Phi)_H = \langle u_b, \gamma_{0,b}(\partial_t \varphi) \rangle_{\Gamma_b}, $$

which is equivalent to

$$ \frac{d\mathcal{E}}{dt} = \frac{d}{dt} \langle u_b, \gamma_{0,b}(\varphi) \rangle_{\Gamma_b} - \langle \partial_t u_b, \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \qquad (2.12) $$

Integrating the above equation from 0 to t and using the vanishing initial conditions yields

$$ \mathcal{E}(t) = \langle u_b(t), \gamma_{0,b}(\varphi(t)) \rangle_{\Gamma_b} - \int_0^t \langle \partial_t u_b(s), \gamma_{0,b}(\varphi(s)) \rangle_{\Gamma_b} ds. \qquad (2.13) $$

One can show that the identity (2.13) holds for the solutions given by Proposition 2.1 (see [21]). The inequality (2.13) is the starting point to derive an estimate of the solution. We give below such estimate.

**Proposition 2.2.** There exists $C > 0$ such that, for any $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$, the solution $\Phi$ to (2.8)-(2.9) satisfies,

$$ \sup_{s \in [0,t]} \mathcal{E}(s) \le C (t^2 + 1) B^2(t), \qquad (2.14) $$

Where B is given by

$$ B(t) = \sup_{s \in [0, t]} \|u_b(s)\|_{H^{-1/2}(\Gamma_b)} + \int_0^t \|\partial_t u_b(s)\|_{H^{-1/2}(\Gamma_b)} \mathrm{d}s. \qquad (2.15) $$

*Proof.* Starting with the equation (2.13), we have for the right-hand side

$$ \mathcal{E}(t) = \langle u_b(t), \gamma_{0,b}(\varphi(t)) \rangle_{\Gamma_b} - \int_0^t \langle \partial_t u_b(s), \gamma_{0,b}(\varphi(s)) \rangle_{\Gamma_b} \mathrm{d}s \le B(t) \sup_{s \in [0,t]} \|\gamma_{0,b}(\varphi)(s)\|_{H^{1/2}(\Gamma_b)}. \qquad (2.16) $$

To estimate the norm on $H^{1/2}(\Gamma_b)$, we use the continuity of the trace and Proposition 1.4,

$$ \sup_{s \in [0,t]} \|\gamma_{0,b}(\varphi)(s)\|_{H^{1/2}(\Gamma_b)} \lesssim \sup_{s \in [0,t]} \|\Phi(s)\|_{D(G^*)}, \qquad (2.17) $$

where we use the symbol $\lesssim$ for inequalities that hold up to a constant independent of $u_b$, $\Phi$ and $t$. We show now that the graph norm in the right-hand side of (2.17) can be bounded by the energy. From the definition of the scalar product in $D(G^*)$ and the energy (2.11), we have

$$ \|\Phi(t)\|_{D(G^*)}^2 = \|\Phi(t)\|_{G}^{2} + \|G^{*}\Phi(t)\|_{H}^{2} \le \|\Phi(t)\|_{G}^{2} + 2\mathcal{E}(t). \qquad (2.18) $$

Since the initial conditions vanish, it holds

$$ \Phi(t) = \int_{0}^{t} \partial_{t} \Phi(r) \mathrm{d}r \quad \Rightarrow \quad \|\Phi(t)\|_{G} \le \int_{0}^{t} \|\partial_{t} \Phi(r)\|_{G} \mathrm{d}r \le \int_{0}^{t} \sqrt{2\mathcal{E}(r)} \mathrm{d}r. $$
## 2.1. Variational formulation

The natural idea for writing the variational formulation associated to (2.1)-(2.2) consists in testing (2.1)-(2.2) against a function $(\tilde{\varphi}, \tilde{\psi}) \in H^1(\Omega) \times L^2(\Omega)$. After integrating by parts, using the boundary conditions (2.3)-(2.4) and the definition of $N$ given by (0.3), we obtain the problem: given $u_b$ regular enough, find

$$ \begin{pmatrix} \varphi \\ \psi \end{pmatrix} \in L^2(0, T; H^1(\Omega) \times L^2(\Omega)), \quad \frac{d}{dt} \begin{pmatrix} \varphi \\ \psi \end{pmatrix} \in L^2(0, T; L^2(\Omega)^2), \qquad (2.5) $$

solution to

$$ \begin{aligned} & \frac{d^2}{dt^2} \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, dx + \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \psi \tilde{\psi} \, dx \\ & \qquad + \int_{\Omega} \rho_0 \left( -\nabla \varphi + N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z \right) \cdot \left( -\nabla \tilde{\varphi} + N \left( \tilde{\psi} + \frac{N}{g} \tilde{\varphi} \right) \mathbf{e}_z \right) \, dx \\ & \qquad + \frac{d^2}{dt^2} \int_{\Gamma_s} \frac{\rho_0}{g} \varphi \tilde{\varphi} \, ds + \int_{\Gamma_b} \rho_0 u_b \gamma_{0,b}(\tilde{\varphi}) \, ds = 0, \quad \forall \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \end{pmatrix} \in H^1(\Omega) \times L^2(\Omega). \end{aligned} \qquad (2.6) $$

The formulation (2.6) will be useful for the numerical approximation. Indeed, the natural spaces for the discretization are classical: $H^1(\Omega)$ for $\varphi$ and $L^2(\Omega)$ for $\psi$. Moreover, for $d = 3$, the velocity formulation has three scalar unknowns whereas the potential formulation requires only two scalar unknowns. Finally, the source $u_b$ appears naturally as a Neumann condition in this formulation.

## 2.2. Existence and uniqueness results

The existence of a solution to (2.1)-(2.2) cannot be directly proven by standard methods (such as [24]), because of the surface condition (2.4) involving the second-order time derivative of $\varphi$. Therefore, we need to introduce a new unknown to the problem, denoted by $\gamma(t) \in L^2(\Gamma_s)$, and we define the vector of unknowns

$$ \Phi(t) = \begin{pmatrix} \varphi(t) \\ \psi(t) \\ \gamma(t) \end{pmatrix} \in \mathcal{G}. \qquad (2.7) $$

From the variational formulation (2.6) we deduce that the problem (2.1) - (2.4) reduces to the following abstract formulation: assume $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$ given, find

$$ \Phi \in L^2(0, T; \mathcal{D}(G^*)), \quad \frac{d}{dt} \Phi \in L^2(0, T; \mathcal{G}), $$

solution to

$$ \frac{d^2}{dt^2} (\Phi(t), \tilde{\Phi})_{\mathcal{G}} + (G^* \Phi(t), G^* \tilde{\Phi})_{\mathcal{H}} = \ell_b(t, \tilde{\Phi}), \quad \forall \tilde{\Phi} \in \mathcal{D}(G^*), \quad \text{in } \mathcal{D}'(0, T), \qquad (2.8) $$

$$ \Phi(0) = \frac{d}{dt} \Phi(0) = 0, \qquad (2.9) $$

where $\ell_b: (0, T) \times \mathcal{D}(G^*) \to \mathbb{R}$ is the linear form

$$ \ell_b(t, \Phi) = \langle u_b(t), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \qquad (2.10) $$

Note that if $\Phi \in \mathcal{D}(G^*)$, then $\gamma$ is the surface trace of $\varphi$, and the equation (2.8) is exactly the equation (2.6). For the formulation (2.8)-(2.9) we have the following result:
Using this inequality to simplify (2.18) yields

$$ \|\Phi(t)\|_{D(G^*)}^2 \le \left( \int_0^t \sqrt{2\mathcal{E}(r)} dr \right)^2 + 2\mathcal{E}(t) \le \left( t \sup_{r \in [0,t]} \sqrt{2\mathcal{E}(r)} \right)^2 + 2\mathcal{E}(t). $$

We finally obtain the bound

$$ \sup_{r \in [0,t]} \|\Phi(r)\|_{D(G^*)}^2 \le 2(t^2 + 1) \sup_{r \in [0,t]} \mathcal{E}(r), $$

hence (2.17), (2.16) and the inequality just above give,

$$ \mathcal{E}(t) \lesssim \sqrt{(t^2 + 1)B} \sup_{r \in [0,t]} \sqrt{\mathcal{E}(r)}, \qquad (2.19) $$

which allows to deduce (2.14) using the Young inequality.

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

The well-posedness of the potential-based formulation (2.8) was obtained by standard tools [24]. The study of the velocity-field formulation is more involved, and is the subject of the next section.

# 3 Analysis of the velocity-field formulation

This section is devoted to prove an existence and uniqueness result for the system of partial differential equations for the velocity-field formulation. We recall the problem at hand: for $u_b$ given, find $\mathbf{U}$ solution to

$$ \rho_0 \frac{\partial^2 \mathbf{U}}{\partial t^2} - \nabla (\rho_0 c_0^2 \nabla \cdot \mathbf{U} - \rho_0 g \mathbf{U} \cdot \mathbf{e}_z) - \nabla \cdot (\rho_0 g \mathbf{U}) \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T], \qquad (3.1) $$

with the boundary conditions

$$ \mathbf{U} \cdot \mathbf{n}_b = u_b \quad \text{on } \Gamma_b \times [0, T], \quad \nabla \cdot \mathbf{U} = 0 \quad \text{on } \Gamma_s \times [0, T], \qquad (3.2) $$

and with vanishing initial conditions.

## 3.1 Variational formulation and uniqueness result

The variational formulation associated to the evolution problem (3.1)-(3.2) is obtained by testing the system (3.1) against a function $\tilde{\mathbf{U}}$ and integrating over $\Omega$. The test function $\tilde{\mathbf{U}}$ is chosen such that its normal trace on $\Gamma_b$ vanishes. Using the boundary conditions (3.2) we obtain:

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \mathbf{U}(t) \tilde{\mathbf{U}} dx + \int_{\Omega} \rho_0 c_0^2 \left( \nabla \cdot \mathbf{U}(t) - \frac{g}{c_0^2} \mathbf{U}(t) \cdot \mathbf{e}_z \right) \left( \nabla \cdot \tilde{\mathbf{U}} - \frac{g}{c_0^2} \tilde{\mathbf{U}} \cdot \mathbf{e}_z \right) dx \\ + \int_{\Omega} \rho_0 N^2 \mathbf{U}(t) \cdot \mathbf{e}_z \tilde{\mathbf{U}} \cdot \mathbf{e}_z dx + \int_{\Gamma_s} \rho_0 g \mathbf{U}(t) \cdot \mathbf{n}_s \tilde{\mathbf{U}} \cdot \mathbf{n}_s ds = 0. \qquad (3.3) $$

The formulation above is completed with the non-homogeneous boundary condition $\mathbf{U} \cdot \mathbf{n}_b = u_b$ on $\Gamma_b$. These formal computations show that the adequate variational formulation to study is the following: assume $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ given, and find

$$ \mathbf{U} \in L^2(0, T; \mathcal{D}(\tilde{G})), \quad \frac{d}{dt} \mathbf{U} \in L^2(0, T; \mathcal{H}), \qquad (3.4) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}, \tilde{\mathbf{U}})_{\mathcal{H}} + (\tilde{G}\mathbf{U}, G\tilde{\mathbf{U}})_g = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.5) $$

$$ \gamma_{1,b}(\mathbf{U}) = u_b, \quad \text{in } (0, T), \qquad (3.6) $$

$$ \mathbf{U}(0) = \frac{d}{dt} \mathbf{U}(0) = 0. \qquad (3.7) $$
**Proposition 2.1.** Assume that $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$. Then the problem (2.8)-(2.9) has a unique solution and, up to a modification on zero measure sets,

$$ \Phi \in C^0([0, T]; \mathcal{D}(G^*)) \cap C^1([0, T]; \mathcal{G}). $$

*Proof.* First we show that, for almost all $t \in (0, T)$, the form $\ell_b(t)$ is a bounded linear functional on $\mathcal{D}(G^*)$. Let $\Phi = (\varphi, \psi, \gamma_{0,s}(\varphi))^{\top} \in \mathcal{D}(G^*)$. From the continuity of the trace operator, there exists a scalar $C_H > 0$ depending only on $H$, such that

$$ |\langle \ell_b(t), \Phi \rangle| = |\langle u_b(t), \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}| \le C_H \|u_b(t)\|_{H^{-1/2}(\Gamma_b)} \|\varphi\|_{H^1(\Omega)}, $$

and from Proposition 1.4 we obtain

$$ |\langle \ell_b(t), \Phi \rangle| \le C_H C_c^{-1} \|u_b(t)\|_{H^{-1/2}(\Gamma_b)} \|\Phi\|_{\mathcal{D}(G^*)}, $$

hence $\ell_b(t)$ is bounded. From the continuity in time of $u_b$, we have $\ell_b \in H^1(0, T; \mathcal{D}(G^*)')$. From this property and (2.8), we deduce that

$$ \frac{d^2}{dt^2} \Phi \in L^2(0, T; \mathcal{D}(G^*)'), $$

hence the initial conditions (2.9) make sense. The existence and uniqueness of a solution to the problem (2.8)-(2.9) follows then directly from standard results; see e.g. Theorem 9.3 in [24] for a result valid for $u_b \in L^2(0, T; H^{-1/2}(\Gamma_b))$ and a proof done by transposition. See also [22] (without proof) for an adapted version of this theorem with more regular source term. $\square$

## 2.3. Energy identity

An energy identity for the system (2.8)-(2.9) is obtained following the usual approach [24]. We define the energy

$$ \mathcal{E}(t) = \frac{1}{2} (\|\partial_t \Phi\|_{\mathcal{G}}^2 + \|G^* \Phi\|_{\mathcal{H}}^2). \quad (2.11) $$

Taking formally $\tilde{\Phi} = \partial_t \Phi$ in the weak formulation (2.8) yields

$$ (\partial_{tt}^2 \Phi, \partial_t \Phi)_\mathcal{G} + (G^* \Phi, G^* \partial_t \Phi)_\mathcal{H} = \langle u_b, \gamma_{0,b}(\partial_t \varphi) \rangle_{\Gamma_b}, $$

which is equivalent to

$$ \frac{d\mathcal{E}}{dt} = \frac{d}{dt} \langle u_b, \gamma_{0,b}(\varphi) \rangle_{\Gamma_b} - \langle \partial_t u_b, \gamma_{0,b}(\varphi) \rangle_{\Gamma_b}. \quad (2.12) $$

Integrating the above equation from 0 to $t$ and using the vanishing initial conditions yields

$$ \mathcal{E}(t) = \langle u_b(t), \gamma_{0,b}(\varphi(t)) \rangle_{\Gamma_b} - \int_0^t \langle \partial_t u_b(s), \gamma_{0,b}(\varphi(s)) \rangle_{\Gamma_b} ds. \quad (2.13) $$

One can show that the identity (2.13) holds for the solutions given by Proposition 2.1 (see Theorem 2.2 in [22]). The equality (2.13) is the starting point to derive an estimate of the solution. We give below such estimate.

**Proposition 2.2.** There exists $C > 0$ such that, for any $u_b \in H^1(0, T; H^{-1/2}(\Gamma_b))$, the solution $\Phi$ to (2.8)-(2.9) satisfies,

$$ \sup_{s \in [0,t]} \mathcal{E}(s) \le C(t^2+1) B^2(t), \quad (2.14) $$

where $B$ is given by

$$ B(t) = \sup_{s \in [0,t]} \|u_b(s)\|_{H^{-1/2}(\Gamma_b)} + \int_0^t \|\partial_t u_b(s)\|_{H^{-1/2}(\Gamma_b)} ds. \quad (2.15) $$
Note that compared to Section 2.2 we have assumed slightly more regularity in time for the source term $u_b$. This is another drawback of the velocity-field formulation and is due to the nature of the condition: here the inhomogeneous boundary condition in (3.2) is of essential type (similar to an inhomogeneous Dirichlet boundary condition). Of course one could weaken this regularity assumption, however, from our current analysis this would also weaken the regularity of the solution. It is rather direct – using Lions-Magenes theory [24] – to prove the following result.

**Theorem 3.1.** The solution to (3.5)-(3.7) is unique.

*Proof.* The problem (3.5)-(3.7) with $u_b \equiv 0$ amounts to find $\mathbf{U}$ with the regularity (3.4) solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}, \tilde{\mathbf{U}})_H + (G\mathbf{U}, G\tilde{\mathbf{U}})_G = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.8) $$

with vanishing initial data. It follows from [24] that this problem has a unique solution and it is zero. $\square$

Existence and stability results with respect to the data $u_b$ are more difficult to obtain because of the essential inhomogeneous boundary condition. The common approach consists in decomposing the solution $\mathbf{U}$ as $\mathbf{U} = \mathbf{U}_0 + L(u_b)$, where the function $\mathbf{U}_0 \in \mathcal{D}(G)$ is solution to a homogeneous problem, and the operator $L$ is a lifting operator. We aim to define a lifting operator in a way that preserves the symmetry between the potential-based and the velocity-field problems. Hence the lifting should be defined as $L(u_b) = -G^*\Phi_b$, where $\Phi_b \in \mathcal{D}(G_\alpha^*)$ is solution to the elliptic problem

$$ \forall \tilde{\Phi} = \begin{pmatrix} \tilde{\Phi} \\ \tilde{\mathbf{V}} \end{pmatrix} \in \mathcal{D}(G^*) \text{ with } \tilde{\Phi} = \begin{pmatrix} \tilde{\Phi} \\ \tilde{\psi} \\ \tilde{\gamma} \end{pmatrix}, \quad (G^*\Phi_b, G^*\tilde{\Phi})_H = \langle u_b, \gamma_{0,b}(\tilde{\Phi}) \rangle_{\Gamma_b}. \qquad (3.9) $$

However, in our case, defining such a lifting is not trivial because of the following result.

**Theorem 3.2.** The range of the operators $G$ and $G^*$ are not closed.

*Proof.* We consider only the case $d = 2$, and show that there is no scalar $C > 0$ such that, for all function $\mathbf{U}(t) \in \mathcal{D}(G) \cap \text{Ker}(G)^\perp$, it holds

$$ \|G\mathbf{U}\|_{G}^{2} \geq C \|\mathbf{U}\|_{H}^{2}. \qquad (3.10) $$

This property implies that the range of $G$ is not closed, hence the range of $G^*$ is also not closed. First note that $\text{Ker}(G) = \{0\}$. Indeed, the kernel of $G$ is defined by

$$ \text{Ker}(G) = \{\mathbf{U} \in \mathcal{D}(G) \mid U_z = 0, \partial_x U_x = 0, \gamma_{1,s}(\mathbf{U}) = 0\}, $$

Hence, for every $\mathbf{U} = (U_x, U_z) \in \text{Ker}(G)$, it holds $U_z = 0$ and $U_x$ is constant in the $x$ direction. Since the domain is infinite in the $x$ direction and $U_x \in L^2(\Omega)$, it is equal to zero. Hence $\mathcal{D}(G^*) \cap \text{Ker}(G)^\perp = \mathcal{D}(G^*)$. Now, let $u \in \mathcal{D}(\Omega)$ be a function with compact support in $\mathbb{R} \times (H_-, H)$, and for each integer $n > 1$ let the function $u_n \in \mathcal{D}(\Omega)$ be defined by $u_n(x) = u(x, n(z - z_0))$, with $z_0 = (H - H_-)/2$. We have

$$ \|\partial_x u_n\|_{L^2(\Omega)} = \frac{1}{n} \|\partial_x u\|_{L^2(\Omega)}, \quad \|\partial_z u_n\|_{L^2(\Omega)} = \|\partial_z u\|_{L^2(\Omega)}. $$

Then let $\mathbf{U}_n \in \mathcal{D}(G)$ be the divergence-free function defined by

$$ \mathbf{U}_n(x, z) = \begin{pmatrix} \partial_z u_n(x, z) \\ -\partial_x u_n(x, z) \end{pmatrix} = \begin{pmatrix} n(\partial_z u)(x, n(z - z_0)) \\ -(\partial_x u)(x, n(z - z_0)) \end{pmatrix}. $$

We have $G\mathbf{U}_n = (-g\partial_x u_n, N\partial_x u_n, 0)^t$ and

$$ \|G\mathbf{U}_n\|_{G}^{2} = \frac{1}{n} \left\| \sqrt{\rho_0 \frac{g^2 + N^2}{c_0^2}} \partial_x u \right\|_{L^2(\Omega)}^{2}, \quad \|\mathbf{U}_n\|_{H}^{2} = \frac{1}{n} \left\| \sqrt{\rho_0} \partial_x u \right\|_{L^2(\Omega)}^{2} + \left\| \sqrt{\rho_0} \partial_z u \right\|_{L^2(\Omega)}^{2}. \qquad (3.11) $$

For $n$ going to infinity we have that $\|G\mathbf{U}_n\|_{G}^{2}$ goes to zero and $\|\mathbf{U}_n\|_{H}^{2}$ converges to $\|\sqrt{\rho_0} \partial_z u\|_{L^2(\Omega)}^{2}$, hence (3.10) can not hold. $\square$
Proof. Starting with the equation (2.13), we have for the right-hand side

$$ \mathcal{E}(t) = \langle u_b(t), \gamma_{0,b}(\varphi(t)) \rangle_{\Gamma_b} - \int_0^t \langle \partial_t u_b(s), \gamma_{0,b}(\varphi(s)) \rangle_{\Gamma_b} ds \le B(t) \sup_{s \in [0,t]} \|\gamma_{0,b}(\varphi)(s)\|_{H^{1/2}(\Gamma_b)}. \quad (2.16) $$

To estimate the norm on $H^{1/2}(\Gamma_b)$, we use the continuity of the trace and Proposition 1.4,

$$ \sup_{s \in [0,t]} \|\gamma_{0,b}(\varphi)(s)\|_{H^{1/2}(\Gamma_b)} \lesssim \sup_{s \in [0,t]} \|\Phi(s)\|_{D(G^*)}, \quad (2.17) $$

where the symbol $\lesssim$ denotes inequalities that hold up to a multiplicative constant independent of $u_b$, $\Phi$ and $t$. We show now that the graph norm in the right-hand side of (2.17) is bounded by the energy. From the definition of the scalar product in $D(G^*)$ and the energy (2.11), we get

$$ \|\Phi(t)\|_{D(G^*)}^2 = \|\Phi(t)\|_{G}^2 + \|G^*\Phi(t)\|_{H}^2 \le \|\Phi(t)\|_{G}^2 + 2\mathcal{E}(t). \quad (2.18) $$

Since the initial conditions vanish, we have

$$ \Phi(t) = \int_0^t \partial_t \Phi(r) dr \quad \Rightarrow \quad \|\Phi(t)\|_G \le \int_0^t \|\partial_t \Phi(r)\|_G dr \le \int_0^t \sqrt{2\mathcal{E}(r)} dr. $$

Using this inequality to simplify (2.18) yields

$$ \|\Phi(t)\|_{D(G^*)}^2 \le \left( \int_0^t \sqrt{2\mathcal{E}(r)} dr \right)^2 + 2\mathcal{E}(t) \le \left( t \sup_{r \in [0,t]} \sqrt{2\mathcal{E}(r)} \right)^2 + 2\mathcal{E}(t). $$

We finally obtain the bound

$$ \sup_{r \in [0,t]} \|\Phi(r)\|_{D(G^*)}^2 \le 2(t^2 + 1) \sup_{r \in [0,t]} \mathcal{E}(r), $$

hence (2.17), (2.16) and the inequality just above give

$$ \mathcal{E}(s) \lesssim \sqrt{2(s^2 + 1)} B(s) \sup_{r \in [0,s]} \sqrt{\mathcal{E}(r)}. $$

The Young inequality applied to the right-hand side of the above equation yields

$$ \mathcal{E}(s) \lesssim (s^2 + 1)B^2(s) + \frac{1}{2} \sup_{r \in [0,s]} \mathcal{E}(r). $$

Taking the supremum over $[0, t]$ of both side of the equation and using the fact that $t \mapsto B(t)$ is increasing, we deduce (2.14). □

The well-posedness of the potential-based formulation (2.8) was obtained by standard tools [24]. The study of the velocity-field formulation is more involved, and is the subject of the next section.

### 3. ANALYSIS OF THE VELOCITY-FIELD FORMULATION

This section is devoted to prove an existence and uniqueness result for the system of partial differential equations for the velocity-field formulation. We recall the problem at hand: for $u_b$ given, find $\mathbf{U}$ solution to

$$ \rho_0 \frac{\partial^2 \mathbf{U}}{\partial t^2} - \nabla (\rho_0 c_0^2 \nabla \cdot \mathbf{U} - \rho_0 g \mathbf{U} \cdot \mathbf{e}_z) - \nabla \cdot (\rho_0 g \mathbf{U}) \mathbf{e}_z = 0, \quad \text{in } \Omega \times [0, T], \quad (3.1) $$
Since the range of G is not closed, the bilinear form of the problem (3.9) is not coercive and the existence of a solution is not ensured. As a result, we cannot define in a straightforward way a lifting operator of the form $L(u_b) = -G^*\Phi_b$. Other lifting choices are possible but are not compatible with the potential-based formulation and are also not straightforward to analyze. To circumvent this problem, we introduce in what follows a dissipative version of the problem (3.5)-(3.7) that will be easier to analyze. Existence and uniqueness results will be obtained using a limit process, namely, by letting the dissipation go to zero.

## 3.2 Existence results

### 3.2.1 A formulation with artificial dissipation

Instead of studying the existence of a solution for problem (3.5)-(3.7), we introduce a modified problem for a new unknown satisfying – assuming that the solution $U(t)$ to (3.5)-(3.7) exists – for $\alpha > 0$

$$ U_{\alpha}(t) = e^{-\alpha t} U(t), \quad \text{in } \mathcal{D}'(0, T). $$

The variational formulation for this new unknown is then: for $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ given, find

$$ U_{\alpha} \in L^2(0, T; \mathcal{D}(\tilde{G})), \quad \frac{d}{dt} U_{\alpha} \in L^2(0, T; \mathcal{H}), \qquad (3.12) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + \alpha^2(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + (\tilde{G}\mathbf{U}_{\alpha}, G\tilde{\mathbf{U}})_{\tilde{G}} = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.13) $$

$$ \gamma_{1,b}(\mathbf{U}_{\alpha}) = e^{-\alpha t} u_b, \quad \text{in } (0, T), \qquad (3.14) $$

$$ \mathbf{U}_{\alpha}(0) = \frac{d}{dt} \mathbf{U}_{\alpha}(0) = 0. \qquad (3.15) $$

The lemma below is straightforward to prove.

**Lemma 3.3.** Let $U_{\alpha}$ be a solution to (3.13)-(3.15). Then $e^{\alpha t} U_{\alpha}$ is the unique solution to (3.5)-(3.7).

As an immediate consequence of the lemma above, we have that problem (3.13)-(3.15) has a unique solution. The relaxed formulation can be written in a more compact form by introducing the space

$$ \mathcal{G} = \mathcal{G} \times \mathcal{H}, \quad (\cdot, \cdot)_{\mathcal{G}} = (\cdot, \cdot)_{\mathcal{H}} + (\cdot, \cdot)_{\mathcal{H}}, $$

and the operator $G_{\alpha}: \mathcal{D}(G_{\alpha}) \subset \mathcal{H} \to \mathcal{G}$, defined by $\mathcal{D}(G_{\alpha}) = \mathcal{D}(G)$ and

$$ \forall \mathbf{U} \in \mathcal{D}(G_{\alpha}), \quad G_{\alpha}\mathbf{U} = \begin{pmatrix} G\mathbf{U} \\ \alpha\mathbf{U} \end{pmatrix}. $$

The operator $\tilde{G}_{\alpha}: \mathcal{D}(\tilde{G}_{\alpha}) \subset \mathcal{H} \to \mathcal{G}$ is defined similarly, and it holds $\mathcal{D}(\tilde{G}_{\alpha}) = \mathcal{D}(\tilde{G})$. The variational formulation (3.13)-(3.15) is then equivalent to: find $U_{\alpha}$ with the regularity (3.12) and solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + (\tilde{G}_{\alpha}\mathbf{U}_{\alpha}, G_{\alpha}\tilde{\mathbf{U}})_{\tilde{G}} = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G_{\alpha}), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.16) $$

$$ \gamma_{1,b}(\mathbf{U}_{\alpha}) = e^{-\alpha t} u_b, \quad \text{in } (0, T), \qquad (3.17) $$

$$ \mathbf{U}_{\alpha}(0) = \frac{d}{dt} \mathbf{U}_{\alpha}(0) = 0. \qquad (3.18) $$
with the boundary conditions

$$ \mathbf{U} \cdot \mathbf{n}_b = u_b \quad \text{on } \Gamma_b \times [0, T], \quad \nabla \cdot \mathbf{U} = 0 \quad \text{on } \Gamma_s \times [0, T], \qquad (3.2) $$

and with vanishing initial conditions.

### 3.1. Variational formulation and uniqueness result

The variational formulation associated to the evolution problem (3.1)-(3.2) is obtained by testing the system (3.1) against a function $\tilde{\mathbf{U}}$ and integrating over $\Omega$. The test function $\tilde{\mathbf{U}}$ is chosen such that its normal trace on $\Gamma_b$ vanishes. Using the boundary conditions (3.2) we obtain:

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \mathbf{U}(t) \tilde{\mathbf{U}} \, dx + \int_{\Omega} \rho_0 c_0^2 \left( \nabla \cdot \mathbf{U}(t) - \frac{g}{c_0^2} \mathbf{U}(t) \cdot \mathbf{e}_z \right) \left( \nabla \cdot \tilde{\mathbf{U}} - \frac{g}{c_0^2} \tilde{\mathbf{U}} \cdot \mathbf{e}_z \right) \, dx \\ + \int_{\Omega} \rho_0 N^2 \mathbf{U}(t) \cdot \mathbf{e}_z \tilde{\mathbf{U}} \cdot \mathbf{e}_z \, dx + \int_{\Gamma_s} \rho_0 g \mathbf{U}(t) \cdot \mathbf{n}_s \tilde{\mathbf{U}} \cdot \mathbf{n}_s \, ds = 0. \quad (3.3) $$

The formulation above is completed with the non-homogeneous boundary condition $\mathbf{U} \cdot \mathbf{n}_b = u_b$ on $\Gamma_b$. These formal computations show that the adequate variational formulation to study is the following: assume $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ given, and find

$$ \mathbf{U} \in L^2(0, T; \mathcal{D}(\tilde{G})), \quad \frac{d}{dt} \mathbf{U} \in L^2(0, T; \mathcal{H}), \qquad (3.4) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}, \tilde{\mathbf{U}})_{\mathcal{H}} + (\tilde{G}\mathbf{U}, G\tilde{\mathbf{U}})_{\mathcal{G}} = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.5) $$

$$ \gamma_{1,b}(\mathbf{U}) = u_b, \quad \text{in } (0, T), \qquad (3.6) $$

$$ \mathbf{U}(0) = \frac{d}{dt} \mathbf{U}(0) = 0. \qquad (3.7) $$

Note that compared to Section 2.2 we have assumed slightly more regularity in time for the source term $u_b$. This is another drawback of the velocity-field formulation and is due to the nature of the condition: here the inhomogeneous boundary condition in (3.2) is of essential type (similar to an inhomogeneous Dirichlet boundary condition). Of course one could weaken this regularity assumption, however, from our current analysis this would also weaken the regularity of the solution.

It is rather direct using Lions-Magenes theory (Theorem 8.1 in [24] for example) to prove that the solution to (3.5)-(3.7) is unique. Existence and stability results with respect to the data $u_b$ are more difficult to obtain because of the essential inhomogeneous boundary condition. The common approach consists in decomposing the solution $\mathbf{U}$ as $\mathbf{U} = \mathbf{U}_0 + L(u_b)$, where the function $\mathbf{U}_0 \in \mathcal{D}(G)$ is solution to a homogeneous problem, and the operator $L$ is a lifting operator. We aim to define a lifting operator in a way that preserves the symmetry between the potential-based and the velocity-field problems. Hence the lifting should be defined as $L(u_b) = -G^*\Phi_b$, where $\Phi_b \in \mathcal{D}(G^*_{\alpha})$ is solution to the elliptic problem

$$ \forall \tilde{\Phi} = \begin{pmatrix} \tilde{\Phi} \\ \tilde{V} \end{pmatrix} \in \mathcal{D}(G^*) \text{ with } \tilde{\Phi} = \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \\ \tilde{\gamma} \end{pmatrix}, \quad (G^*\Phi_b, G^*\tilde{\Phi})_{\mathcal{H}} = \langle u_b, \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_b}. \qquad (3.8) $$

However, in our case, defining such a lifting is not trivial because of the following result.
### 3.2.2 Lifting operator for the dissipative problem

The operators $G_\alpha$ and $\tilde{G}_\alpha$ are densely defined and closed, so are their adjoints $G_\alpha^*$ and $\tilde{G}_\alpha^*$. We have, in particular,

$$ \mathcal{D}(G_\alpha^*) = \mathcal{D}(G^*) \times \mathcal{H}, \quad \text{and} \quad \forall (\Phi, \mathbf{V}) \in \mathcal{D}(G_\alpha^*), \quad G_\alpha^* \begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} = G^*\Phi + \alpha\mathbf{V} \in \mathcal{H}. \quad (3.19) $$

Moreover $\mathcal{D}(\tilde{G}_\alpha^*) = \mathcal{D}(\tilde{G}^*) \times \mathcal{H}$ and $G_\alpha^*$ is an extension of $\tilde{G}_\alpha^*$. In the following the identity operators in $\mathcal{H}$ and in $\mathbf{G}$ will be used. By abuse of notation, both are denoted by $\mathbf{1}$.

**Theorem 3.4.** The range of the operators $G_\alpha$ and $G_\alpha^*$ are closed.

*Proof.* Since $G_\alpha = \begin{pmatrix} G \\ \alpha \mathbf{1} \end{pmatrix}$, we have $\text{Ker } G_\alpha = \emptyset$. Then, for $\mathbf{U} \in \mathcal{D}(G_\alpha) \cap \text{Ker}(G_\alpha)^\perp = \mathcal{D}(G_\alpha)$, it holds

$$ \|G_\alpha \mathbf{U}\|_{\mathcal{G}}^2 = \|G\mathbf{U}\|_{\mathcal{G}}^2 + \|\alpha \mathbf{U}\|_{\mathcal{H}}^2 \geq \alpha^2 \|\mathbf{U}\|_{\mathcal{H}}^2, $$

which concludes the proof.

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

This result is key for constructing a lifting operator. We recall that the kernel of the operator $G_\alpha^*$ is a closed subspace of $\mathbf{G}$ and we introduce $Q_\alpha \in \mathcal{L}(\mathbf{G})$, the orthogonal projection on $\text{Ker}(G_\alpha^*)$. Then $Q_\alpha^2 = Q_\alpha$ and $Q_\alpha$ is self-adjoint.

**Theorem 3.5.** For $\alpha \leq 1$ we have

$$ \forall \mathbf{\Phi} = \begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} \in \mathcal{D}(G_\alpha^*), \quad \|Q_\alpha \mathbf{\Phi}\|_{\mathcal{G}}^2 + \|G_\alpha^* \mathbf{\Phi}\|_{\mathcal{H}}^2 \geq \alpha^2 \|\mathbf{\Phi}\|_{\mathcal{G}}^2. $$

*Proof.* Let $\mathbf{\Phi} \in \mathcal{D}(G_\alpha^*)$. We have $\|Q_\alpha \mathbf{\Phi}\|_{\mathcal{G}}^2 + \|(I - Q_\alpha)\mathbf{\Phi}\|_{\mathcal{G}}^2 = \|\mathbf{\Phi}\|_{\mathcal{G}}^2$, and therefore

$$ \|G_\alpha^* \mathbf{\Phi}\|_{\mathcal{H}}^2 = \|G_\alpha^* (I - Q_\alpha) \mathbf{\Phi}\|_{\mathcal{H}}^2 \geq \alpha^2 \|(I - Q_\alpha) \mathbf{\Phi}\|_{\mathcal{G}}^2. $$

Hence we have

$$ \|Q_\alpha \mathbf{\Phi}\|_{\mathcal{G}}^2 + \|G_\alpha^* \mathbf{\Phi}\|_{\mathcal{H}}^2 \geq \|Q_\alpha \mathbf{\Phi}\|_{\mathcal{G}}^2 + \alpha^2 \|(I - Q_\alpha) \mathbf{\Phi}\|_{\mathcal{G}}^2 \geq \alpha^2 \|\mathbf{\Phi}\|_{\mathcal{G}}^2. $$

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

The exact expression of $Q_\alpha$ is of no practical interest in what follows. It has to be noted however that it is non trivial, due to the fact that $\text{Ker } G_\alpha^*$ is an infinite dimensional space. As a direct consequence of Proposition 1.4, we have the following bound for the $H^1(\Omega)$-norm,

**Lemma 3.6.** Let $\alpha < 1$, there exists a scalar $C_c > 0$ such that

$$ \forall \mathbf{\Phi} = \begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} \in \mathcal{D}(G_\alpha^*) \text{ with } \mathbf{\Phi} = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix}, \quad \|\mathbf{\Phi}\|_{\mathcal{G}}^2 + \|G_\alpha^* \mathbf{\Phi}\|_{\mathcal{H}}^2 \geq C_c \|\varphi\|_{H^1(\Omega)}^2. $$

*Proof.* By definition and the triangle inequality,

$$ \|\mathbf{\Phi}\|_{\mathcal{G}}^2 + \|G_\alpha^* \mathbf{\Phi}\|_{\mathcal{H}}^2 \geq \|\mathbf{\Phi}\|_{\mathcal{G}}^2 + \|\mathbf{V}\|_{\mathcal{H}}^2 + (\|G^*\mathbf{\Phi}\|_{\mathcal{H}} - \alpha \|\mathbf{V}\|_{\mathcal{H}})^2. $$

The conclusion follows from the Young's inequality, the assumption on $\alpha$ and Proposition 1.4.

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

Aiming to define the aforementioned lifting operator of a data $u_b \in H^{-1/2}(\Gamma_b)$ (the dependence in time is omitted), we introduce the following problem: for all $u_b$, find $\mathbf{\Phi}_b \in \mathcal{D}(G_\alpha^*)$ solution to

$$ \forall \tilde{\mathbf{\Phi}} = \begin{pmatrix} \tilde{\Phi} \\ \tilde{\mathbf{V}} \end{pmatrix} \in \mathcal{D}(G_\alpha^*) \text{ with } \tilde{\mathbf{\Phi}} = \begin{pmatrix} \tilde{\varphi} \\ \tilde{\psi} \\ \tilde{\gamma} \end{pmatrix}, \quad (Q_\alpha \mathbf{\Phi}_b, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Phi}_b, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = \langle u_b, \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_b}. \quad (3.20) $$

Thanks to Theorem 3.5 and Lemma 3.6, it is classical to show that Equation (3.20) admits a unique solution $\mathbf{\Phi}_b \in \mathcal{D}(G_\alpha^*)$ that depends continuously on $u_b$ (it is a standard application of Lax-Milgram Lemma). We then construct the lifting operator $L_\alpha \in \mathcal{L}(H^{-1/2}(\Gamma_b), \mathcal{H})$ by setting

$$ \forall u_b \in H^{-1/2}(\Gamma_b), \quad L_\alpha(u_b) = -G_\alpha^* \mathbf{\Phi}_b. $$
**Theorem 3.1.** The range of the operators G and G* are not closed.

For conciseness, the proof of Theorem 3.1 is given in Appendix 7.1.2. Since the range of G is not closed, the bilinear form of the problem (3.8) is not coercive and the existence of a solution is not ensured. Hence, one cannot define in a straightforward way a lifting operator of the form $L(u_b) = -G^*\Phi_b$. Other lifting choices are possible but are not compatible with the potential-based formulation in the sense that our strategy to prove the equivalence between the velocity and potential-based formualtion relies on the existence of such lifting.

To circumvent this problem, we introduce in what follows a dissipative version of the problem (3.5)-(3.7) that will be easier to analyze. Existence and uniqueness results for the non-dissipative problem will be obtained using a limit process, by letting the dissipation go to zero.

## 3.2. Existence results

### 3.2.1. A formulation with artificial dissipation

Instead of studying the existence of a solution for problem (3.5)-(3.7), we introduce a modified problem for a new unknown satisfying – assuming that the solution $U(t)$ to (3.5)-(3.7) exists – for $\alpha > 0$

$$U_{\alpha}(t) = e^{-\alpha t} U(t), \quad \text{in } \mathcal{D}'(0, T).$$

The variational formulation for this new unknown is then: for $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ given, find

$$U_{\alpha} \in L^2(0, T; \mathcal{D}(\tilde{G})), \quad \frac{d}{dt} U_{\alpha} \in L^2(0, T; \mathcal{H}), \qquad (3.9)$$

solution to

$$\frac{d^2}{dt^2}(\mathbf{U}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(\mathbf{U}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + \alpha^2(\mathbf{U}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + (\tilde{G}\mathbf{U}_\alpha, G\tilde{\mathbf{U}})_\mathcal{G} = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G), \quad \text{in } \mathcal{D}'(0, T), \qquad (3.10)$$

$$\gamma_{1,b}(\mathbf{U}_{\alpha}) = e^{-\alpha t} u_{b}, \quad \text{in } (0,T), \qquad (3.11)$$

$$\mathbf{U}_{\alpha}(0) = \frac{d}{dt} \mathbf{U}_{\alpha}(0) = 0. \qquad (3.12)$$

The lemma below is straightforward to prove.

**Lemma 3.2.** We let $U_{\alpha}$ denote a solution to (3.10)-(3.12). Then $e^{\alpha t} U_{\alpha}$ is the unique solution to (3.5)-(3.7).

As an immediate consequence of the lemma above, we have that problem (3.10)-(3.12) has a unique solution. The relaxed formulation can be written in a more compact form by introducing the space

$$G = G \times H, \quad (\cdot, \cdot)_G = (\cdot, \cdot)_G + (\cdot, \cdot)_H,$$

and the operator $G_{\alpha}: \mathcal{D}(G_{\alpha}) \subset \mathcal{H} \to G$, defined by $\mathcal{D}(G_{\alpha}) = \mathcal{D}(G)$ and

$$\forall \mathbf{U} \in \mathcal{D}(G_{\alpha}), \quad G_{\alpha}\mathbf{U} = \begin{pmatrix} G\mathbf{U} \\ \alpha\mathbf{U} \end{pmatrix}.$$
**Proposition 3.7.** The function $L_\alpha(u_b) \in \mathcal{H}$ has the following properties:

$$L_\alpha(u_b) \in \mathcal{D}(\tilde{G}_\alpha) = \mathcal{D}(\tilde{G}), \quad \tilde{G}_\alpha L_\alpha(u_b) = Q_\alpha \Phi_b \in \text{Ker}(G_\alpha^*), \quad \text{and} \quad \gamma_{1,b}(L_\alpha(u_b)) = u_b.$$

*Proof.* In the equation (3.20) we choose a test function $\tilde{\Phi}$ in the space $\mathcal{D}(\tilde{G}_\alpha^*)$. Since $\gamma_{0,b}(\tilde{\varphi}) = 0$, we obtain

$$(G_{\alpha}^{*} \mathbf{\Phi}_{b}, \tilde{G}_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}}, \quad \forall \tilde{\mathbf{\Phi}} \in \mathcal{D}(\tilde{G}_{\alpha}^{*}). \qquad (3.21)$$

This implies

$$G_{\alpha}^{*} \mathbf{\Phi}_{b} \in \mathcal{D}(\tilde{G}_{\alpha}^{**}) = \mathcal{D}(\tilde{G}_{\alpha}) = \mathcal{D}(\tilde{G}) \quad \text{and} \quad \tilde{G}_{\alpha}^{**} G_{\alpha}^{*} \mathbf{\Phi}_{b} = \tilde{G}_{\alpha} G_{\alpha}^{*} \mathbf{\Phi}_{b} = -Q_{\alpha} \mathbf{\Phi}_{b} = -\tilde{G}_{\alpha} L_{\alpha}(u_{b}), \qquad (3.22)$$

where we have used that $\tilde{G}_\alpha^{**} = \tilde{G}_\alpha$ since $\tilde{G}_\alpha$ is closed and densely defined. The first two properties of the proposition are proved. To prove the last property, we use the abstract Green's formula of Lemma 1.2: for all $\tilde{\Phi} \in \mathcal{D}(G_\alpha^*)$ we have,

$$(G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = (G_{\alpha}^{*} \mathbf{\Phi}_{b}, G^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} + (G_{\alpha}^{*} \mathbf{\Phi}_{b}, \alpha \tilde{\mathbf{V}})_{\mathcal{H}} \qquad (3.23)$$

$$= (\tilde{G}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \tilde{\mathbf{\Phi}})_{\mathcal{H}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_{b}} + (G_{\alpha}^{*} \mathbf{\Phi}_{b}, \alpha \tilde{\mathbf{V}})_{\mathcal{H}}. \qquad (3.24)$$

Hence we have

$$(G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = (\tilde{G}_{\alpha}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \tilde{\mathbf{\Phi}})_{\mathcal{G}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_{b}},$$

and with the equations (3.20) and (3.22) we obtain

$$(G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_{b}} \qquad (3.25)$$

$$= -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} + \langle u_{b}, \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_{b}}, \qquad (3.26)$$

from which we deduce that $\gamma_{1,b}(L_\alpha(u_b)) = -\gamma_{1,b}(G_\alpha^* \Phi_b) = u_b$. $\square$

**Remark 3.8.** Assume now that $u_b \in L^2(0, T; H^{-1/2}(\Gamma_b))$. The operator $L_\alpha$ being continuous from $H^{-1/2}(\Gamma_b)$ to $\mathcal{H}$, one can deduce that

$$L_{\alpha}(u_{b}) \in L^{2}(0, T; \mathcal{H}).$$

Proposition 3.7 can be extended: we have, for almost all $t \in (0, T)$,

$$L_{\alpha}(u_{b}(t)) \in \mathcal{D}(\tilde{G}_{\alpha}), \quad \tilde{G}_{\alpha} L_{\alpha}(u_{b}(t)) \in \text{Ker} \, G_{\alpha}^{*}, \quad \text{and} \quad \gamma_{1,s}(L_{\alpha}(u_{b}(t))) = u_{b}(t).$$

Moreover, the regularity in time of $L_\alpha(u_b)$ depends straightforwardly on the regularity in time of $u_b$. In particular,

$$u_b \in H^k(0, T; H^{-1/2}(\Gamma_b)) \Rightarrow L_\alpha(u_b) \in H^k(0, T; \mathcal{D}(\tilde{G})). \qquad (3.27)$$

### 3.2.3 Existence result

We show now the existence of solution to the dissipative problem.

**Proposition 3.9.** If $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ and $u_b(0) = \frac{d}{dt}u_b(0) = 0$, the problem (3.16)-(3.18) admits a unique solution $\mathbf{U}_\alpha$. It satisfies, up to modifications on zero-measure sets,

$$\mathbf{U}_{\alpha} \in C^{1}([0, T]; \mathcal{H}) \cap C^{0}([0, T]; \mathcal{D}(\tilde{G})). \qquad (3.28)$$

*Proof.* Using the lifting operator, we look first for $\mathbf{U}_{\alpha,0}(t) \in \mathcal{D}(G)$ such that, for all $\tilde{\mathbf{U}} \in \mathcal{D}(G)$,

$$\begin{aligned} & \frac{d^2}{dt^2}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_{\mathcal{H}} + (G_{\alpha} \mathbf{U}_{\alpha,0}, G_{\alpha} \tilde{\mathbf{U}})_{\mathcal{G}} \\ & = \left( \frac{d^2}{dt^2}(e^{-\alpha t} L_{\alpha}(u_b)) + 2\alpha \frac{d}{dt}(e^{-\alpha t} L_{\alpha}(u_b)), \tilde{\mathbf{U}} \right)_{\mathcal{H}} \\ & = e^{-\alpha t} \left( \frac{d^2}{dt^2} L_{\alpha}(u_b) - \alpha^2 L_{\alpha}(u_b), \tilde{\mathbf{U}} \right)_{\mathcal{H}} \quad \text{in } \mathcal{D}'(0, T), \end{aligned} \qquad (3.29)$$
The operator $\tilde{G}_{\alpha} : \mathcal{D}(\tilde{G}_{\alpha}) \subset \mathcal{H} \rightarrow \mathcal{G}$ is defined similarly, and we have $\mathcal{D}(\tilde{G}_{\alpha}) = \mathcal{D}(\tilde{G})$. The variational formulation (3.10)-(3.12) is then equivalent to: find $\mathbf{U}_{\alpha}$ with the regularity (3.9) and solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + (\tilde{G}_{\alpha}\mathbf{U}_{\alpha}, G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} = 0, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G_{\alpha}), \quad \text{in } \mathcal{D}'(0, T), \quad (3.13) $$

$$ \gamma_{1,b}(\mathbf{U}_{\alpha}) = e^{-\alpha t}u_b, \quad \text{in } (0,T), \quad (3.14) $$

$$ \mathbf{U}_{\alpha}(0) = \frac{d}{dt}\mathbf{U}_{\alpha}(0) = 0. \quad (3.15) $$

### 3.2.2. Lifting operator for the dissipative problem

The operators $G_{\alpha}$ and $\tilde{G}_{\alpha}$ are densely defined and closed, so are their adjoints $G_{\alpha}^{*}$ and $\tilde{G}_{\alpha}^{*}$. We have, in particular,

$$ \mathcal{D}(G_{\alpha}^{*}) = \mathcal{D}(G^{*}) \times \mathcal{H}, \quad \text{and} \quad \forall (\Phi, \mathbf{V}) \in \mathcal{D}(G_{\alpha}^{*}), \quad G_{\alpha}^{*}\begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} = G^{*}\Phi + \alpha \mathbf{V} \in \mathcal{H}. \quad (3.16) $$

Moreover $\mathcal{D}(\tilde{G}_{\alpha}^{*}) = \mathcal{D}(\tilde{G}^{*}) \times \mathcal{H}$ and $G_{\alpha}^{*}$ is an extension of $\tilde{G}_{\alpha}^{*}$. As a direct consequence of Proposition 1.4, we have the following bound for the $H^1(\Omega)$-norm.

**Lemma 3.3.** Let $\alpha < 1$, there exists a scalar $C_c > 0$ such that

$$ \forall \mathbf{\Phi} = \begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} \in \mathcal{D}(G_{\alpha}^{*}) \text{ with } \mathbf{\Phi} = \begin{pmatrix} \varphi \\ \psi \\ \gamma \end{pmatrix}, \quad \|\mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|G_{\alpha}^{*}\mathbf{\Phi}\|_{\mathcal{H}}^{2} \geq C_c \|\varphi\|_{H^1(\Omega)}^{2}. $$

*Proof.* From the definition of $G_{\alpha}^{*}$ and by the triangle inequality,

$$ \|\mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|G_{\alpha}^{*}\mathbf{\Phi}\|_{\mathcal{H}}^{2} \geq \|\mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|\mathbf{V}\|_{\mathcal{H}}^{2} + (\|G^{*}\mathbf{\Phi}\|_{\mathcal{H}} - \alpha\|\mathbf{V}\|_{\mathcal{H}})^{2}. $$

The conclusion follows from Young's inequality, the assumption on $\alpha$ and Proposition 1.4. $\square$

The following result is key for constructing a lifting operator.

**Theorem 3.4.** The range of the operators $G_{\alpha}$ and $G_{\alpha}^{*}$ are closed.

*Proof.* We have $\ker(G_{\alpha}) = \{0\}$ by the definition of $G_{\alpha}$. Hence, for $\mathbf{U} \in \mathcal{D}(G_{\alpha}) \cap \ker(G_{\alpha})^{\perp} = \mathcal{D}(G_{\alpha})$, we have

$$ \|G_{\alpha}\mathbf{U}\|_{\mathcal{G}}^{2} = \|G\mathbf{U}\|_{\mathcal{G}}^{2} + \|\alpha\mathbf{U}\|_{\mathcal{H}}^{2} \geq \alpha^{2}\|\mathbf{U}\|_{\mathcal{H}}^{2}, $$

which concludes the proof. $\square$

Next the orthogonal projection on $\ker G_{\alpha}^{*}$, denoted by $Q_{\alpha} \in \mathcal{L}(\mathcal{G})$ is introduced. The projection is well defined since the kernel of the operator $G_{\alpha}^{*}$ is a closed subspace of $\mathcal{G}$. Moreover, we have $Q_{\alpha}^{2} = Q_{\alpha}$ and $Q_{\alpha}$ is self-adjoint. We also have the inequality

**Theorem 3.5.** For $\alpha \leq 1$ we have

$$ \forall \mathbf{\Phi} = \begin{pmatrix} \Phi \\ \mathbf{V} \end{pmatrix} \in \mathcal{D}(G_{\alpha}^{*}), \quad \|Q_{\alpha}\mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|G_{\alpha}^{*}\mathbf{\Phi}\|_{\mathcal{H}}^{2} \geq \alpha^{2}\|\mathbf{\Phi}\|_{\mathcal{G}}^{2}. $$

*Proof.* Let $\mathrm{Id}$ denote the identity operators both in $\mathcal{H}$ and in $\mathcal{G}$, and take $\mathbf{\Phi} \in \mathcal{D}(G_{\alpha}^{*})$. Since $Q_{\alpha}$ is an orthogonal projection operator, we have $\|Q_{\alpha}\mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|(\mathrm{Id} - Q_{\alpha})\mathbf{\Phi}\|_{\mathcal{G}}^{2} = \|\mathbf{\Phi}\|_{\mathcal{G}}^{2}$, and therefore, using Lemma 3.4,

$$ \|G_{\alpha}^{*}\mathbf{\Phi}\|_{\mathcal{H}}^{2} = \|G_{\alpha}^{*}(\mathrm{Id} - Q_{\alpha})\mathbf{\Phi}\|_{\mathcal{H}}^{2} \geq \alpha^{2} \|(\mathrm{Id} - Q_{\alpha})\mathbf{\Phi}\|_{\mathcal{G}}^{2}. $$
with vanishing initial conditions. From the remark above, the data is in $L^2(0, T; \mathcal{D}(\tilde{G}))$. The existence and uniqueness of $\mathbf{U}_{\alpha,0}$ is then obtained by application of standard results [9], and $\mathbf{U}_{\alpha,0}$ has the following regularity,

$$ \mathbf{U}_{\alpha,0} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(\tilde{G})). \qquad (3.30) $$

Let $\mathbf{U}_{\alpha} = \mathbf{U}_{\alpha,0} + L_{\alpha}(e^{-\alpha t}u_b)$. Since we have – up to modifications on zero-measure sets – that the inclusion of $H^2(0, T; \mathcal{H})$ into $C^1([0, T]; \mathcal{H})$ (see [11]), we deduce that

$$ \mathbf{U}_{\alpha} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(\tilde{G})). \qquad (3.31) $$

Finally we show that $\mathbf{U}_{\alpha}$ is solution to (3.16). The computations are mostly straightforward, but we give some details for the following term,

$$ (\tilde{G}_{\alpha}\mathbf{U}_{\alpha}, G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} = (\tilde{G}_{\alpha}\mathbf{U}_{\alpha,0}, G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} + (\tilde{G}_{\alpha}L_{\alpha}(e^{-\alpha t}), G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} \qquad (3.32) $$

Since $G_{\alpha} \subset \tilde{G}_{\alpha}$, it holds $\tilde{G}_{\alpha}\mathbf{U}_{\alpha,0} = G_{\alpha}\mathbf{U}_{\alpha,0}$, hence the first term of (3.32) is replaced using (3.29). For the second term of (3.32), from the properties of the lifting operator given in Proposition 3.7, it holds

$$ (\tilde{G}_{\alpha}L_{\alpha}(e^{-\alpha t}u_b), G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} = e^{-\alpha t}(Q_{\alpha}\Phi_b, G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} = e^{-\alpha t}(G_{\alpha}^*Q_{\alpha}\Phi_b, \tilde{\mathbf{U}})_{\mathcal{H}} = 0, \qquad (3.33) $$

which concludes the proof. □

# 4 From potential-based solutions to velocity-field solutions

In Sections 2 and 3, we have proved that both the velocity-field problem (3.5)-(3.7) and the potential-based problem (2.8)-(2.9) have a unique solution. The aim of this section is to prove that velocity-fields solution to (3.5)-(3.7) can be constructed from potential-based solutions. More precisely, we will prove the following theorem,

**Theorem 4.1.** Let $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ and $u_b(0) = \frac{d}{dt}u_b(0) = 0$, and let $\Phi$ be the solution to the problem (2.8)-(2.9) with source term $u_b$. Then $\mathbf{U} = G^*\Phi$ is the unique solution to the problem (3.5)-(3.7) with the same source term.

A key ingredient in the following analysis is the Von Neumann theorem [18], that we recall below for the sake of completeness.

**Theorem 4.2** (Von Neumann). If $T: \mathcal{D}(T) : \mathcal{H} \to \mathcal{G}$ is a closed densely defined operator, then $T^*T$ is self-adjoint and $\mathcal{D}(T^*T)$ is dense in $\mathcal{D}(T)$.

As already shown in Sec. 3.1, the solution to the problem (3.5)-(3.7) should be defined as the sum of the solution to a homogeneous problem and the lifting of the source term. As presented in Section 3.2.1, the definition of the lifting requires in turn to consider the dissipative problem (3.16)-(3.18). For this reason, it is more convenient to prove first the equivalence between the dissipative problem (3.16)-(3.18) and the corresponding potential-based formulation, that we introduce in the following section. The result stated in Theorem 4.1 for the non-dissipative problem (2.8)-(2.9) and (3.5)-(3.7) is then deduced through a limit process carried out at the end of this section.

## 4.1 The case with artificial dissipation

In this section we recall the dissipative problem (3.16), and define the associated “dual” or “adjoint” problem. For simplicity, we first consider volume sources, and show then how to deduce a result similar to Theorem 4.1 for problems with boundary sources. The dissipative formulation with a volume source is: let $F_{\mathbf{U}}$ regular enough, find

$$ \mathbf{U}_{\alpha} \in L^2(0, T; \mathcal{D}(\tilde{G}_{\alpha})), \quad \frac{d}{dt}\mathbf{U}_{\alpha} \in L^2(0, T; \mathcal{H}), \qquad (4.1) $$
Hence we have

$$ \|Q_{\alpha} \mathbf{\Phi}\|_{\mathcal{G}}^{2} + \|G_{\alpha}^{*} \mathbf{\Phi}\|_{\mathcal{H}}^{2} \geq \|Q_{\alpha} \mathbf{\Phi}\|_{\mathcal{G}}^{2} + \alpha^{2} \|(\mathrm{Id} - Q_{\alpha}) \mathbf{\Phi}\|_{\mathcal{G}}^{2} \geq \alpha^{2} \|\mathbf{\Phi}\|_{\mathcal{G}}^{2}. \quad \square $$

The exact expression of $Q_\alpha$ is of no practical interest in what follows. It has to be noted however that it is non trivial, due to the fact that Ker $G_\alpha^*$ is an infinite dimensional space.

Aiming to define the aforementioned lifting operator of a data $u_b \in H^{-1/2}(\Gamma_b)$ (the dependence in time is omitted), we introduce the following problem: for all $u_b$, find $\mathbf{\Phi}_b \in \mathcal{D}(G_\alpha^*)$ solution to

$$ \forall \tilde{\mathbf{\Phi}} = \begin{pmatrix} \tilde{\mathbf{\Phi}} \\ \tilde{\mathbf{V}} \end{pmatrix} \in \mathcal{D}(G_{\alpha}^{*}) \text{ with } \tilde{\mathbf{\Phi}} = \begin{pmatrix} \tilde{\mathbf{\varphi}} \\ \tilde{\mathbf{\psi}} \\ \tilde{\gamma} \end{pmatrix}, \quad (Q_{\alpha} \mathbf{\Phi}_b, \tilde{\mathbf{\Phi}})_{\mathcal{G}} + (G_{\alpha}^{*} \mathbf{\Phi}_b, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = \langle u_b, \gamma_{0,b}(\tilde{\mathbf{\varphi}}) \rangle_{\Gamma_b}. \quad (3.17) $$

Thanks to Theorem 3.5 and Lemma 3.3, it is classical to show that Equation (3.17) admits a unique solution $\mathbf{\Phi}_b \in \mathcal{D}(G_\alpha^*)$ that depends continuously on $u_b$ (it is a standard application of the Lax-Milgram Lemma). We then construct the lifting operator $L_\alpha \in \mathcal{L}(H^{-1/2}(\Gamma_b), \mathcal{H})$ by setting

$$ \forall u_b \in H^{-1/2}(\Gamma_b), \quad L_\alpha(u_b) = -G_\alpha^* \mathbf{\Phi}_b. $$

**Proposition 3.6.** The function $L_\alpha(u_b) \in \mathcal{H}$ has the following properties:

$$ L_{\alpha}(u_b) \in \mathcal{D}(\tilde{G}_{\alpha}) = \mathcal{D}(\tilde{G}), \quad \tilde{G}_{\alpha} L_{\alpha}(u_b) = Q_{\alpha} \mathbf{\Phi}_{b} \in \mathrm{Ker}(G_{\alpha}^{*}), \quad \text{and} \quad \gamma_{1,b}(L_{\alpha}(u_b)) = u_b. $$

*Proof.* In the equation (3.17) we choose a test function $\tilde{\mathbf{\Phi}}$ in the space $\mathcal{D}(\tilde{G}_\alpha^*)$. Since $\gamma_{0,b}(\tilde{\mathbf{\varphi}}) = 0$, we obtain

$$ (G_{\alpha}^{*} \mathbf{\Phi}_{b}, \tilde{G}_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}}, \quad \forall \tilde{\mathbf{\Phi}} \in \mathcal{D}(\tilde{G}_{\alpha}^{*}). \quad (3.18) $$

This implies

$$ G_{\alpha}^{*} \mathbf{\Phi}_{b} \in \mathcal{D}(\tilde{G}_{\alpha}^{**}) = \mathcal{D}(\tilde{G}_{\alpha}) = \mathcal{D}(\tilde{G}) \quad \text{and} \quad \tilde{G}_{\alpha}^{**} G_{\alpha}^{*} \mathbf{\Phi}_{b} = \tilde{G}_{\alpha} G_{\alpha}^{*} \mathbf{\Phi}_{b} = -Q_{\alpha} \mathbf{\Phi}_{b} = -\tilde{G}_{\alpha} L_{\alpha}(u_{b}), \quad (3.19) $$

where we have used that $\tilde{G}_\alpha^{**} = \tilde{G}_\alpha$ since $\tilde{G}_\alpha$ is closed and densely defined. The first two properties of the proposition are proved. To prove the last property, we use the abstract Green formula of Lemma 1.2: for all $\tilde{\mathbf{\Phi}} \in \mathcal{D}(G_\alpha^*)$ we have,

$$ (G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = (G_{\alpha}^{*} \mathbf{\Phi}_{b}, G^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} + (G_{\alpha}^{*} \mathbf{\Phi}_{b}, \alpha \tilde{\mathbf{V}})_{\mathcal{H}} \quad (3.20) $$

$$ = (\tilde{G}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \tilde{\mathbf{\Phi}})_{\mathcal{H}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\mathbf{\varphi}}) \rangle_{\Gamma_b} + (G_{\alpha}^{*} \mathbf{\Phi}_{b}, \alpha \tilde{\mathbf{V}})_{\mathcal{H}}. \quad (3.21) $$

Hence we have

$$ (G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = (\tilde{G}_{\alpha}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \tilde{\mathbf{\Phi}})_{\mathcal{G}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\mathbf{\varphi}}) \rangle_{\Gamma_b}, $$

and with the equations (3.17) and (3.19) we obtain

$$ (G_{\alpha}^{*} \mathbf{\Phi}_{b}, G_{\alpha}^{*} \tilde{\mathbf{\Phi}})_{\mathcal{H}} = -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} - \langle \gamma_{1,b}(G_{\alpha}^{*} \mathbf{\Phi}_{b}), \gamma_{0,b}(\tilde{\mathbf{\varphi}}) \rangle_{\Gamma_b} \quad (3.22) $$

$$ = -(Q_{\alpha} \mathbf{\Phi}_{b}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} + \langle u_{b}, \gamma_{0,b}(\tilde{\mathbf{\varphi}}) \rangle_{\Gamma_b}, \quad (3.23) $$

from which we deduce that $\gamma_{1,b}(L_\alpha(u_b)) = -\gamma_{1,b}(G_\alpha^* \mathbf{\Phi}_b) = u_b$. $\square$
solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(\mathbf{U}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + (G_\alpha \mathbf{U}_\alpha, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = (F_\alpha, \tilde{\mathbf{U}})_\mathcal{H}, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G_\alpha), \quad \text{in } \mathcal{D}'(0, T), \quad (4.2) $$

$$ \mathbf{U}_{\alpha}(0) = \frac{d}{dt} \mathbf{U}_{\alpha}(0) = 0. \quad (4.3) $$

The potential-based formulation is defined using the adjoint $G_\alpha^*$, and reads: given $F_\Phi$ regular enough, find

$$ \Phi_{\alpha} \in L^2(0, T; \mathcal{D}(G_{\alpha}^{*})), \quad \frac{d}{dt} \Phi_{\alpha} \in L^2(0, T; \mathcal{G}) $$

solution to: for all $\tilde{\Phi} \in \mathcal{D}(G_\alpha^*)$

$$ \frac{d^2}{dt^2}(\mathbf{\Phi}_{\alpha}, \tilde{\mathbf{\Phi}})_\mathcal{G} + 2\alpha \frac{d}{dt}(\mathbf{\Phi}_{\alpha}, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Phi}_{\alpha}, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = (F_\Phi, \tilde{\mathbf{\Phi}})_\mathcal{G}, \quad \mathcal{D}'(0, T), \quad (4.4) $$

$$ \mathbf{\Phi}_{\alpha}(0) = \frac{d}{dt} \mathbf{\Phi}_{\alpha}(0) = 0. \quad (4.5) $$

The existence and uniqueness of the solution to the problems (4.2)-(4.3) and (4.4)-(4.5) are a direct application of well-known results, given e.g. in [9]. Whenever $F_\alpha \in L^2(0, T; \mathcal{H})$ and $F_\alpha \in L^2(0, T; \mathcal{G})$, there exists a unique solution $\mathbf{U}_\alpha$ to (4.2)-(4.3), and a unique solution $\mathbf{\Phi}_\alpha$ to (4.4)-(4.5), and they satisfy

$$ \mathbf{U}_{\alpha} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(G)), \quad \mathbf{\Phi}_{\alpha} \in C^1([0, T]; \mathcal{G}) \cap C^0([0, T]; \mathcal{D}(G_{\alpha}^{*})). $$

We state now a relation property between problems (4.2)-(4.3) and (4.4)-(4.5).

**Theorem 4.3.** Let $\mathbf{\Phi}_\alpha \in \mathcal{D}(G_\alpha^*)$ be the solution to the problem (4.4)-(4.5). If $F_\Phi \in L^2(0, T; \mathcal{D}(G_\alpha^*))$, then $\mathbf{U}_\alpha = G_\alpha^* \mathbf{\Phi}_\alpha$ is the unique solution to the problem (4.2)-(4.3), with source term defined by $F_\alpha = G_\alpha^* F_\Phi$.

*Proof.* We introduce the function $\Psi_\alpha$ defined by

$$ \Psi_{\alpha}(t) = \int_{0}^{t} \mathbf{\Phi}_{\alpha}(s) \, ds. $$

Integrating in time (4.4) yields an equation for $\Psi_\alpha$, for all $\tilde{\Phi} \in \mathcal{D}(G_\alpha^*)$,

$$ \frac{d^2}{dt^2}(\mathbf{\Psi}_{\alpha}, \tilde{\mathbf{\Phi}})_\mathcal{G} + 2\alpha \frac{d}{dt}(\mathbf{\Psi}_{\alpha}, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Psi}_{\alpha}, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = \left( \int_0^t F_\Phi(s) \, ds, \tilde{\mathbf{\Phi}} \right)_\mathcal{G}, \quad (4.6) $$

and since we have by construction $\Psi_\alpha \in C^2([0, T]; \mathcal{G}) \cap C^1([0, T]; \mathcal{D}(G_\alpha^*))$, we deduce from (4.6) that

$$ \Psi_{\alpha} \in C^0([0, T]; \mathcal{D}(G_{\alpha} G_{\alpha}^{*})). \quad (4.7) $$

Now, for all $\tilde{\mathbf{U}} \in \mathcal{D}(G_\alpha^* G_\alpha)$ we set $\tilde{\mathbf{\Phi}} = G_\alpha \tilde{\mathbf{U}}$ in (4.6). We obtain, using the definition of the adjoint of $G_\alpha$, Equation (4.7) and the assumption that $F_\Phi \in L^2(0, T; \mathcal{D}(G^*))$,

$$ \frac{d^2}{dt^2}(G_\alpha^* \mathbf{\Psi}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(G_\alpha^* \mathbf{\Psi}_\alpha, \tilde{\mathbf{U}})_\mathcal{H} + (G_\alpha G_\alpha^* \mathbf{\Psi}_\alpha, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = \left( \int_0^t G_\alpha^* F_\Phi(t') \, dt', \tilde{\mathbf{U}} \right)_\mathcal{H}. \quad (4.8) $$

Therefore, setting $\mathbf{V}_\alpha = G_\alpha^* \mathbf{\Psi}_\alpha$, we have that

$$ \mathbf{V}_{\alpha} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(G)) $$

is solution to: for all $\tilde{\mathbf{U}} \in \mathcal{D}(G_\alpha^* G_\alpha)$,

$$ \frac{d^2}{dt^2}(\mathbf{V}_{\alpha}, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(\mathbf{V}_{\alpha}, \tilde{\mathbf{U}})_\mathcal{H} + (G_\alpha \mathbf{V}_{\alpha}, G_\alpha \tilde{\mathbf{U}})_\mathcal{H} = \left( \int_0^t G_\alpha^* F_\Phi(t') \, dt', \tilde{\mathbf{U}} \right)_\mathcal{H}. \quad (4.9) $$
**Remark 3.7.** Assume now that $u_b \in L^2(0, T; H^{-1/2}(\Gamma_b))$. The operator $L_\alpha$ being continuous from $H^{-1/2}(\Gamma_b)$ to $\mathcal{H}$, one can deduce that

$$L_\alpha(u_b) \in L^2(0, T; \mathcal{H}).$$

Proposition 3.6 can be extended: we have, for almost all $t \in (0, T)$,

$$L_\alpha(u_b(t)) \in \mathcal{D}(\tilde{G}_\alpha), \quad \tilde{G}_\alpha L_\alpha(u_b(t)) \in \text{Ker } G_\alpha^*, \quad \text{and} \quad \gamma_{1,s}(L_\alpha(u_b(t))) = u_b(t).$$

Moreover, the regularity in time of $L_\alpha(u_b)$ depends straightforwardly on the regularity in time of $u_b$. In particular,

$$u_b \in H^k(0, T; H^{-1/2}(\Gamma_b)) \Rightarrow L_\alpha(u_b) \in H^k(0, T; \mathcal{D}(\tilde{G})). \quad (3.24)$$

### 3.2.3. Existence result

We show now the existence of solution to the dissipative problem.

**Proposition 3.8.** If $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ and $u_b(0) = \frac{d}{dt}u_b(0) = 0$, the problem (3.13)-(3.15) admits a unique solution $\mathbf{U}_\alpha$. It satisfies, up to modifications on zero-measure sets,

$$\mathbf{U}_\alpha \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(\tilde{G})). \quad (3.25)$$

*Proof.* Using the lifting operator, we look first for $\mathbf{U}_{\alpha,0}(t) \in \mathcal{D}(G)$ such that, for all $\tilde{\mathbf{U}} \in \mathcal{D}(G)$,

$$\begin{aligned} & \frac{d^2}{dt^2}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_\mathcal{H} + (G_\alpha \mathbf{U}_{\alpha,0}, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} \\ &= \left( \frac{d^2}{dt^2}(e^{-\alpha t} L_\alpha(u_b)) + 2\alpha \frac{d}{dt}(e^{-\alpha t} L_\alpha(u_b)), \tilde{\mathbf{U}} \right)_\mathcal{H} \\ &= e^{-\alpha t} \left( \frac{d^2}{dt^2} L_\alpha(u_b) - \alpha^2 L_\alpha(u_b), \tilde{\mathbf{U}} \right)_\mathcal{H} \quad \text{in } \mathcal{D}'(0, T), \end{aligned} \quad (3.26)$$

with vanishing initial conditions. From the remark above, the data is in $L^2(0, T; \mathcal{D}(\tilde{G}))$. The existence and uniqueness of $\mathbf{U}_{\alpha,0}$ is then obtained by application of standard results [10], and we have

$$\mathbf{U}_{\alpha,0} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(G)). \quad (3.27)$$

Let $\mathbf{U}_\alpha = \mathbf{U}_{\alpha,0} + L_\alpha(e^{-\alpha t}u_b)$. Since we have – up to modifications on zero-measure sets – that the inclusion of $H^2(0, T; \mathcal{H})$ into $C^1([0, T]; \mathcal{H})$ (see [12]), we deduce that

$$\mathbf{U}_\alpha \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(\tilde{G})). \quad (3.28)$$

Finally we show that $\mathbf{U}_\alpha$ is solution to (3.13). The computations are mostly straightforward; we give some details for the following term:

$$(\tilde{G}_\alpha \mathbf{U}_\alpha, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = (\tilde{G}_\alpha \mathbf{U}_{\alpha,0}, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} + (\tilde{G}_\alpha L_\alpha(e^{-\alpha t}), G_\alpha \tilde{\mathbf{U}})_\mathcal{G} \quad (3.29)$$

Since $\tilde{G}_\alpha$ is an extension of $G_\alpha$, the equality $\tilde{G}_\alpha \mathbf{U}_{\alpha,0} = G_\alpha \mathbf{U}_{\alpha,0}$ holds. Hence, the first term of (3.29) can be replaced using (3.26). From the Proposition 3.6, the second term of (3.29) satisfies

$$(\tilde{G}_\alpha L_\alpha(e^{-\alpha t} u_b), G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = e^{-\alpha t} (Q_\alpha \Phi_b, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = e^{-\alpha t} (G_\alpha^* Q_\alpha \Phi_b, \tilde{\mathbf{U}})_\mathcal{H} = 0. \quad (3.30)$$

□
Since $G_\alpha$ is a closed densely defined operator, $\mathcal{D}(G_\alpha^*G_\alpha)$ is dense in $\mathcal{D}(G_\alpha)$ by the Von Neumann theorem, hence (4.9) can be extended to functions in $\tilde{\mathbf{U}} \in \mathcal{D}(G_\alpha) = \mathcal{D}(G)$. At this point we have shown that

$$ \mathbf{V}_\alpha = \int_0^t G_\alpha^* \mathbf{\Phi}_\alpha(s) \mathrm{d}s $$

is solution to (4.2)-(4.3), with source term given by

$$ F_U(t) = \int_0^t G_\alpha^* F_\Phi(s) \mathrm{d}s. $$

Since $F_U$ is differentiable in time by construction, one can deduce that the time derivative of $\mathbf{V}_\alpha$ is solution to (4.2) – (4.3) with source term $G_\alpha^*F_\Phi$, which allows us to conclude. □

We consider now the case of a source term located at the bottom, namely the dissipative problem (3.16)-(3.17) and its potential-based formulation: find

$$ \mathbf{\Phi}_\alpha \in L^2(0, T; \mathcal{D}(G_\alpha^*)), \quad \frac{d}{dt} \mathbf{\Phi}_\alpha \in L^2(0, T; \mathbf{G}) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{\Phi}_\alpha, \tilde{\mathbf{\Phi}})_\mathcal{G} + 2\alpha \frac{d}{dt}(\mathbf{\Phi}_\alpha, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Phi}_\alpha, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = \langle e^{-\alpha t} u_b, \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_b}, \quad \forall \tilde{\mathbf{\Phi}} \in \mathcal{D}(G_\alpha^*), \quad \text{in } \mathcal{D}'(0, T), \qquad (4.10) $$

$$ \mathbf{\Phi}_\alpha(0) = \frac{d}{dt} \mathbf{\Phi}_\alpha(0) = 0. \qquad (4.11) $$

As soon as the right-hand side of (4.10) defines a functional in $H^2(0, T; \mathcal{D}(G_\alpha^*)')$ (see the proof of Proposition 2.1), the solution to (4.10)-(4.11) exists and is unique. The result below is deduced from Theorem 4.3 above and the property of the lifting operator given in Proposition 3.7.

**Proposition 4.4.** Let $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ and $u_b(0) = \frac{d}{dt}u_b(0) = 0$, and let $\mathbf{\Phi}_\alpha$ be the solution to the problem (4.10)-(4.11). Then $\mathbf{U}_\alpha = G_\alpha^*\mathbf{\Phi}_\alpha$ is the unique solution to the problem (3.16)-(3.18).

*Proof.* Let $\mathbf{\Phi}_b(t)$ be the solution to (3.20) associated to the data $u_b(t)$ for almost all time $t \in [0, T]$. We have $\mathbf{\Phi}_b \in H^2(0, T; \mathcal{D}(G_\alpha^*))$ and the equation (4.10) can be written

$$ \frac{d^2}{dt^2}(\mathbf{\Phi}_\alpha, \tilde{\mathbf{\Phi}})_\mathcal{G} + 2\alpha \frac{d}{dt}(\mathbf{\Phi}_\alpha, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Phi}_\alpha, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = e^{-\alpha t}(Q_\alpha \mathbf{\Phi}_b, \tilde{\mathbf{\Phi}})_\mathcal{G} + e^{-\alpha t}(G_\alpha^* \mathbf{\Phi}_b, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H}. \quad (4.12) $$

Let $\mathbf{\Phi}_0 = \mathbf{\Phi}_\alpha - e^{-\alpha t}\mathbf{\Phi}_b$ and define

$$ \begin{aligned} F_{\mathbf{\Phi}} &= -\frac{d^2}{dt^2}(e^{-\alpha t}\mathbf{\Phi}_b) - 2\alpha \frac{d}{dt}(e^{-\alpha t}\mathbf{\Phi}_b) + e^{-\alpha t}Q_\alpha\mathbf{\Phi}_b \\ &= e^{-\alpha t}(\alpha^2\mathbf{\Phi}_b - \frac{d^2}{dt^2}\mathbf{\Phi}_b + Q_\alpha\mathbf{\Phi}_b) \in L^2(0, T; \mathcal{D}(G_\alpha^*)), \end{aligned} $$

then $\mathbf{\Phi}_0$ is the unique solution to: for all $\tilde{\mathbf{\Phi}} \in \mathcal{D}(G_\alpha^*)$,

$$ \frac{d^2}{dt^2}(\mathbf{\Phi}_0, \tilde{\mathbf{\Phi}})_\mathcal{G} + 2\alpha \frac{d}{dt}(\mathbf{\Phi}_0, \tilde{\mathbf{\Phi}})_\mathcal{G} + (G_\alpha^* \mathbf{\Phi}_0, G_\alpha^* \tilde{\mathbf{\Phi}})_\mathcal{H} = (F_{\mathbf{\Phi}}, \tilde{\mathbf{\Phi}})_\mathcal{G}, $$

with vanishing initial data (since $u_b$ and $\frac{d}{dt}u_b$ vanish at the initial time). Using Theorem 4.3, the function $\mathbf{U}_{\alpha,0}$ defined by $\mathbf{U}_{\alpha,0} = G_\alpha^*\mathbf{\Phi}_0 \in \mathcal{D}(G_\alpha)$ is solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_\mathcal{H} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha,0}, \tilde{\mathbf{U}})_\mathcal{H} + (G_\alpha \mathbf{U}_{\alpha,0}, G_\alpha \tilde{\mathbf{U}})_\mathcal{G} = (G_\alpha^* F_{\mathbf{\Phi}}, \tilde{\mathbf{U}})_\mathcal{G}, $$
## 4. FROM POTENTIAL-BASED SOLUTIONS TO VELOCITY-FIELD SOLUTIONS

In Sections 2 and 3, we have proved that both the velocity-field problem (3.5)-(3.7) and the potential-based problem (2.8)-(2.9) have a unique solution. The aim of this section is to prove that velocity-field solution to (3.5)-(3.7) can be constructed from potential-based solutions. More precisely, we will prove the following theorem:

**Theorem 4.1.** Let $u_b \in H^2(0, T; H^{-1/2}(\Gamma_b))$ and $u_b(0) = \frac{d}{dt}u_b(0) = 0$, and let $\Phi$ be the solution to the problem (2.8)-(2.9) with source term $u_b$. Then $\mathbf{U} = G^*\Phi$ is the unique solution to the problem (3.5)-(3.7) with the same source term.

A key ingredient in the following analysis is the Von Neumann theorem [18], that we recall below for the sake of completeness.

**Theorem 4.2** (Von Neumann). If $T: \mathcal{D}(T): \mathcal{H} \to \mathcal{G}$ is a closed densely defined operator, then $T^*T$ is self-adjoint and $\mathcal{D}(T^*T)$ is dense in $\mathcal{D}(T)$.

As already presented in Section 3.1, the solution to the problem (3.5)-(3.7) should be defined as the sum of the solution to a homogeneous problem and the lifting of the source term. We showed in Section 3.2.1 that the definition of the lifting requires in turn to consider the dissipative problem (3.13)-(3.15). For this reason, we first prove the equivalence between the dissipative problem (3.13)-(3.15) and the corresponding potential-based problem, that we introduce in the next section. The Theorem 4.1 for the non-dissipative problems is then deduced from a limit process carried out at the end of this section.

### 4.1. The case with artificial dissipation

In this section we recall the dissipative problem (3.13), and define the associated potential-based problem. For simplicity, we first consider volume sources, and show then how to deduce a result similar to Theorem 4.1 for problems with boundary sources. The dissipative formulation with a volume source reads: For $F_{\mathbf{U}} \in L^2(0, T; \mathcal{H})$, find

$$ \mathbf{U}_{\alpha} \in L^2(0, T; \mathcal{D}(\tilde{G}_{\alpha})), \quad \frac{d}{dt} \mathbf{U}_{\alpha} \in L^2(0, T; \mathcal{H}), \qquad (4.1) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(\mathbf{U}_{\alpha}, \tilde{\mathbf{U}})_{\mathcal{H}} + (G_{\alpha}\mathbf{U}_{\alpha}, G_{\alpha}\tilde{\mathbf{U}})_{\mathcal{G}} = (F_{\mathbf{U}}, \tilde{\mathbf{U}})_{\mathcal{H}}, \quad \forall \tilde{\mathbf{U}} \in \mathcal{D}(G_{\alpha}), \quad \text{in } \mathcal{D}'(0, T), \qquad (4.2) $$

$$ \mathbf{U}_{\alpha}(0) = \frac{d}{dt} \mathbf{U}_{\alpha}(0) = 0. \qquad (4.3) $$

The potential-based formulation is defined using the adjoint $G_{\alpha}^*$, and reads: Given $F_{\Phi} \in L^2(0, T; \mathcal{G})$, find

$$ \mathbf{\Phi}_{\alpha} \in L^2(0, T; \mathcal{D}(G_{\alpha}^*)), \quad \frac{d}{dt} \mathbf{\Phi}_{\alpha} \in L^2(0, T; \mathcal{G}) \qquad (4.4) $$

solution to: for all $\tilde{\mathbf{\Phi}} \in \mathcal{D}(G_{\alpha}^*)$

$$ \frac{d^2}{dt^2}(\mathbf{\Phi}_{\alpha}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} + 2\alpha \frac{d}{dt}(\mathbf{\Phi}_{\alpha}, \tilde{\mathbf{\Phi}})_{\mathcal{G}} + (G_{\alpha}^* \mathbf{\Phi}_{\alpha}, G_{\alpha}^* \tilde{\mathbf{\Phi}})_{\mathcal{H}} = (F_{\Phi}, \tilde{\mathbf{\Phi}})_{\mathcal{G}}, \quad \mathcal{D}'(0, T), \qquad (4.5) $$

$$ \mathbf{\Phi}_{\alpha}(0) = \frac{d}{dt} \mathbf{\Phi}_{\alpha}(0) = 0. \qquad (4.6) $$

The existence and uniqueness of solutions to the problems (4.2)-(4.3) and (4.5)-(4.6) are a direct application of well-known results, given e.g. in [10]. There exists a unique solution $\mathbf{U}_{\alpha}$ to (4.2)-(4.3), and a unique solution $\mathbf{\Phi}_{\alpha}$ to (4.5)-(4.6), and they satisfy

$$ \mathbf{U}_{\alpha} \in C^1([0, T]; \mathcal{H}) \cap C^0([0, T]; \mathcal{D}(G)), \quad \mathbf{\Phi}_{\alpha} \in C^1([0, T]; \mathcal{G}) \cap C^0([0, T]; \mathcal{D}(G_{\alpha}^*)). $$
where, by construction,

$$G_{\alpha}^{*} F_{\Phi} = e^{-\alpha t} \left( \frac{d^2}{dt^2} L_{\alpha}(u_b) - \alpha^2 L_{\alpha}(u_b) \right).$$

This shows that $U_{\alpha,0} = G_{\alpha}^{*} \Phi_{0}$ is indeed the solution to (3.29). It remains only to define, as in the proof of Theorem 3.9, the function $U_{\alpha} = U_{\alpha,0} + e^{-\alpha t} L_{\alpha}(u_b) = U_{\alpha,0} - e^{-\alpha t} G_{\alpha}^{*} \Phi_{b}$ and observe that it is solution to (3.16) with the required regularity. □

We have shown that for problems with artificial dissipation, the solution $U_{\alpha}$ to the velocity-field problem (3.16)-(3.18) is obtained from the solution $\Phi_{\alpha}$ to the potential-based problem (4.10)-(4.11), thanks to the relation $U_{\alpha} = G_{\alpha}^{*} \Phi_{\alpha}$. The next step is to prove a similar relation for the dissipative free problems (3.5)-(3.16) and (2.8)-(2.9). It should be noted that the relation between the dissipative potential-based problem (4.10)-(4.11) and its non dissipative version (2.8)-(2.9) are not related as obviously as for the velocity-field problems. In particular, one does not have a relation of the form $\Phi(t) = e^{\alpha t} \Phi_{\alpha}(t)$, where $\Phi$ is the solution to (2.8)-(2.9) and $\Phi_{\alpha}$ is the solution to (4.10)-(4.11). Deducing $\Phi$ from $\Phi_{\alpha}$ requires to study the limit process $\alpha \to 0$.

## 4.2 The dissipation free case

At this point we have shown that $U$, the solution to (3.5)-(3.7), can be expressed as

$$U = e^{\alpha t} U_{\alpha} = e^{\alpha t} G_{\alpha}^{*} \Phi_{\alpha}.$$

The main idea is now to pass to the limit when $\alpha$ goes to zero in the equality above, to conclude that $U = G^{*} \Phi$, where $\Phi$ is solution to (2.8)-(2.9). Note that the weak convergence

$$G_{\alpha}^{*} \Phi_{\alpha} \xrightarrow{L^{2}(0,T;H)} G^{*} \Phi$$

would imply that $U = G^{*} \Phi$. It is then sufficient to prove that the weak convergence holds and that the limit $\Phi$ is solution to (2.8)-(2.9).

To prove the mentioned convergence results we first provide preliminary uniform (with respect to $\alpha$) energy estimates for the solution $\Phi_{\alpha}$. Following the same proof as Proposition 2.2, the following estimation can be proved,

$$\frac{1}{2} \| \frac{d}{dt} \Phi_{\alpha}(t) \|_{G}^{2} + 2\alpha \int_{0}^{t} \| \Phi_{\alpha}(t') \|_{G}^{2} dt' + \frac{1}{2} \| G_{\alpha}^{*} \Phi_{\alpha}(t) \|_{H}^{2} \lesssim B_{\alpha}^{2}(t), \quad (4.13)$$

where the inequality above holds up to a constant independent of $\Phi_{\alpha}$ or the source term $u_{b}$, and with

$$B_{\alpha}(t) = \sup_{s \in [0, t]} \|e^{-\alpha s} u_{b}(t)\|_{H^{-1/2}(\Gamma_{b})} + \int_{0}^{t} \|\partial_{t}(e^{-\alpha s} u_{b}(s))\|_{H^{-1/2}(\Gamma_{b})} \mathrm{d}s.$$

For $\alpha \le 1$, it holds

$$B_{\alpha}(t) \le D(t) := C(t) + \int_{0}^{t} \|u_{b}(s)\|_{H^{-1/2}(\Gamma_{b})} \mathrm{d}s,$$

where $B(t)$ is defined in (2.15). Therefore, the right-hand side of (4.13) can be replaced by a positive function independent of $\alpha$. For $\alpha \le 1$, we have

$$\| \frac{d}{dt} \Phi_{\alpha}(t) \|_{G} + \| G_{\alpha}^{*} \Phi_{\alpha}(t) \|_{H} \lesssim C(t). \quad (4.14)$$

From (4.10), we also deduce an estimate of the second order time derivative in the Hilbert space $D(G_{\alpha}^{*})'$,

$$\| \frac{d^2}{dt^2} \Phi_{\alpha}(t) \|_{D(G_{\alpha}^{*})'} \lesssim C(t). \quad (4.15)$$

Those preliminary observations allow us to state the following result.
We state now a relation property between problems (4.2)-(4.3) and (4.5)-(4.6).

**Theorem 4.3.** Let $\Phi_{\alpha} \in \mathcal{D}(G_{\alpha}^{*})$ be the solution to the problem (4.5)-(4.6). If $F_{\Phi} \in L^{2}(0, T; \mathcal{D}(G_{\alpha}^{*}))$, then $U_{\alpha} = G_{\alpha}^{*}\Phi_{\alpha}$ is the unique solution to the problem (4.2)-(4.3), with source term defined by $F_{U} = G_{\alpha}^{*}F_{\Phi}$.

*Proof.* We introduce the function $\Psi_{\alpha}$ defined by

$$ \Psi_{\alpha}(t) = \int_{0}^{t} \Phi_{\alpha}(s) \, ds. $$

Integrating in time (4.5) yields an equation for $\Psi_{\alpha}$, for all $\tilde{\Phi} \in \mathcal{D}(G_{\alpha}^{*})$,

$$ \frac{d^2}{dt^2}(\Psi_{\alpha}, \tilde{\Phi})_{\mathcal{G}} + 2\alpha \frac{d}{dt}(\Psi_{\alpha}, \tilde{\Phi})_{\mathcal{G}} + (G_{\alpha}^{*}\Psi_{\alpha}, G_{\alpha}^{*}\tilde{\Phi})_{\mathcal{H}} = \left(\int_{0}^{t} F_{\Phi}(s) \, ds, \tilde{\Phi}\right)_{\mathcal{G}}, \quad (4.7) $$

and since we have by construction $\Psi_{\alpha} \in C^{2}([0, T]; \mathcal{G}) \cap C^{1}([0, T]; \mathcal{D}(G_{\alpha}^{*}))$, we deduce from (4.7) that

$$ \Psi_{\alpha} \in C^{0}([0, T]; \mathcal{D}(G_{\alpha}G_{\alpha}^{*})). \quad (4.8) $$

Now, for all $\tilde{U} \in \mathcal{D}(G_{\alpha}^{*}G_{\alpha})$ we set $\tilde{\Phi} = G_{\alpha}\tilde{U}$ in (4.7). Using Equation (4.8), the definition of the adjoint of $G_{\alpha}$ and the regularity of $F_{U}$ we obtain

$$ \frac{d^2}{dt^2}(G_{\alpha}^{*}\Psi_{\alpha}, \tilde{U})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(G_{\alpha}^{*}\Psi_{\alpha}, \tilde{U})_{\mathcal{H}} + (G_{\alpha}G_{\alpha}^{*}\Psi_{\alpha}, G_{\alpha}\tilde{U})_{\mathcal{G}} = \left(\int_{0}^{t} G_{\alpha}^{*}F_{\Phi}(t') \, dt', \tilde{U}\right)_{\mathcal{H}}. \quad (4.9) $$

Therefore, setting $V_{\alpha} = G_{\alpha}^{*}\Psi_{\alpha}$, we have $V_{\alpha} \in C^{1}([0, T]; \mathcal{H}) \cap C^{0}([0, T]; \mathcal{D}(G))$ is solution to

$$ \frac{d^2}{dt^2}(V_{\alpha}, \tilde{U})_{\mathcal{H}} + 2\alpha \frac{d}{dt}(V_{\alpha}, \tilde{U})_{\mathcal{H}} + (G_{\alpha}V_{\alpha}, G_{\alpha}\tilde{U})_{\mathcal{H}} = \left(\int_{0}^{t} G_{\alpha}^{*}F_{\Phi}(t') \, dt', \tilde{U}\right)_{\mathcal{H}}, \quad \forall \tilde{U} \in \mathcal{D}(G_{\alpha}^{*}G_{\alpha}). \quad (4.10) $$

Since $G_{\alpha}$ is a closed densely defined operator, $\mathcal{D}(G_{\alpha}^{*}G_{\alpha})$ is dense in $\mathcal{D}(G_{\alpha})$ by the Von Neumann theorem, hence (4.10) can be extended to functions in $\tilde{U} \in \mathcal{D}(G_{\alpha}) = \mathcal{D}(G)$. At this point we have shown that

$$ V_{\alpha} = \int_{0}^{t} G_{\alpha}^{*} \Phi_{\alpha}(s) \, ds $$

is solution to (4.2)-(4.3), with source term given by $F_{U}(t) = \int_{0}^{t} G_{\alpha}^{*}F_{\Phi}(s) \, ds$. By construction $F_{U}$ is differentiable in time, hence one can deduce that the time derivative of $V_{\alpha}$ is solution to (4.2) – (4.3) with source term $G_{\alpha}^{*}F_{\Phi}$. $\square$

We consider now the case of a source located at the bottom. The velocity-field problem is (3.13)-(3.15). The associated potential-based formulation is: Find $\Phi$ with the regularity (4.4) solution to

$$ \frac{d^2}{dt^2}(\Phi_{\alpha}, \tilde{\Phi})_{\mathcal{G}} + 2\alpha \frac{d}{dt}(\Phi_{\alpha}, \tilde{\Phi})_{\mathcal{G}} + (G_{\alpha}^{*}\Phi_{\alpha}, G_{\alpha}^{*}\tilde{\Phi})_{\mathcal{H}} = \langle e^{-\alpha t} u_{b}, \gamma_{0,b}(\tilde{\varphi}) \rangle_{\Gamma_{b}}, \quad \forall \tilde{\Phi} \in \mathcal{D}(G_{\alpha}^{*}), \text{ in } \mathcal{D}'(0, T), \quad (4.11) $$

and the vanishing initial conditions (4.6). The solution to (4.11) exists and is unique – see the proof of Proposition 2.1. The next result uses the lifting operator studied in Section 3.2.2.

**Proposition 4.4.** Let $u_{b} \in H^{2}(0, T; H^{-1/2}(\Gamma_{b}))$ and $u_{b}(0) = \frac{d}{dt}u_{b}(0) = 0$, and let $\Phi_{\alpha}$ be the solution to the problem (4.11)-(4.6). Then $U_{\alpha} = G_{\alpha}^{*}\Phi_{\alpha}$ is the unique solution to the problem (3.13)-(3.15).
**Lemma 4.5.** The functions $G_{\alpha}^{*}\Phi_{\alpha}$ for $\alpha \in \mathbb{R}^{+}$ converges weakly in $L^{2}(0, T; \mathcal{H})$ to $G^{*}\Phi$ when $\alpha \to 0$ with $\Phi$ solution to (2.8)-(2.9).

*Proof.* Using the estimates (4.14)-(4.15) and the fact that $\Phi_{\alpha}$ vanishes at the initial time, we have that $\Phi_{\alpha}$, when considered as a sequence in $\alpha$, is bounded in the Hilbert space

$$W := L^2(0, T; \mathcal{D}(G_{\alpha}^{*})).$$

Note that the domain of the operator $G_{\alpha}^{*}$ is independent of $\alpha$, from the definition $\mathcal{D}(G_{\alpha}^{*}) = \mathcal{D}(G^{*}) \times \mathcal{H}$. Hence, we deduce that $\Phi_{\alpha}$ converges weakly – up to a subsequence – to a function $\Phi$ in $W$. Decomposing $\Phi = (\Phi \mathbf{V})^{t}$ and passing to the limit, $\alpha \to 0$ in the formulation (4.10), we obtain

$$\frac{d^2}{dt^2}(\Phi, \tilde{\Phi})_G + (G^*\Phi, G^*\tilde{\Phi})_H = \langle u_b, \tilde{\phi} \rangle_{\Gamma_b} \quad \forall \tilde{\Phi} = (\Phi 0)^t \in \mathcal{D}(G^*) \times \mathcal{H}, \quad \text{in } \mathcal{D}'(0, T), \quad (4.16)$$

and

$$\frac{d^2}{dt^2}(\mathbf{V}, \tilde{\mathbf{V}})_H = 0 \quad \forall \tilde{\mathbf{\Phi}} = (0 \tilde{\mathbf{V}})^t \in \mathcal{D}(G^*) \times \mathcal{H}, \quad \text{in } \mathcal{D}'(0, T). \quad (4.17)$$

Hence $\Phi$ is solution to (2.8). This almost concludes the proof. Indeed, only initial conditions must be investigated to conclude. From the energy estimate (4.14) we also have the weak convergence in $H^1(0, T; \mathbf{G})$, therefore, for all $\tilde{\mathbf{\Phi}} \in \mathbf{G}$,

$$(\mathbf{\Phi}_{\alpha}(0), \tilde{\mathbf{\Phi}})_{\mathbf{G}} = \frac{1}{T} \int_{0}^{T} \frac{d}{dt}((t-T)(\mathbf{\Phi}_{\alpha}, \tilde{\mathbf{\Phi}})_{\mathbf{G}}) dt \underset{\alpha \to 0}{\longrightarrow} \frac{1}{T} \int_{0}^{T} \frac{d}{dt}((t-T)(\mathbf{\Phi}, \tilde{\mathbf{\Phi}})_{\mathbf{G}}) dt = (\mathbf{\Phi}(0), \tilde{\mathbf{\Phi}})_{\mathbf{G}},$$

hence $\mathbf{\Phi}(0) = \mathbf{\Phi}_{\alpha}(0) = 0$. Similarly, thanks to the weak convergence in $H^2(0, T; (\mathcal{D}(G^*) \times \mathcal{H})')$ we conclude that

$$\frac{d}{dt} \mathbf{\Phi}(0) = \frac{d}{dt} \mathbf{\Phi}_{\alpha}(0) = 0.$$

□

We have shown that $\mathbf{\Phi}_{\alpha}$ converges weakly to $\Phi$, that the limit $\Phi$ is solution to (2.8)-(2.9), and that $G_{\alpha}^{*}\mathbf{\Phi}_{\alpha}$ converges weakly to $G^{*}\Phi$. Then $\mathbf{U} = G^{*}\Phi$ is solution to (3.5)-(3.7), which proves Theorem 4.1.

# 5 Numerical illustration

We now use the two formulations (3.5)-(3.7) and (2.6) to give several numerical illustrations of tsunamis and hydro-acoustic waves generated by an earthquake or landslide source. Three sets of simulations are presented in this section. The first one reproduces a classical scenario of an submarine earthquake generating a tsunami and acoustic waves. This scenario is used for several validations: convergence analysis, validating the model by comparison with the literature, and illustrating the equivalence of the potential-based and velocity-field formulations. In the second set of simulations, the equation $\mathbf{U} = G^{*}\Phi$ is used to investigate the validity of the classical irrotational assumption for the velocity field. The third simulation is a preliminary work towards the study of acoustic waves generated by submarine landslides. We compute the spectrogram of pressure recorded at a given sensor for a source located on the seabed which emits a range of frequencies typical of those observed in the field. This novel simulation highlights the interference pattern caused by the reflection of the acoustic waves at the surface.

We first describe the choice of parameters. Then the discretization of the velocity-field and the potential-based simulations are presented. Finally, for each set of simulations, we describe the scenario and provide some illustrations.

As stated in the introduction, the background functions $\rho_0, c_0$ must satisfy the following positivity property,

$$-\frac{g}{\rho_0} \frac{d\rho_0}{dz} - \frac{g^2}{c_0^2} = N^2 \geq 0. \quad (5.1)$$
Proof. Let $\Phi_b(t)$ be the solution to (3.17) associated to the data $u_b(t)$ for almost all time $t \in [0, T]$. We have $\Phi_b \in H^2(0, T; D(G_\alpha^*))$ and the equation (4.11) can be written

$$ \frac{d^2}{dt^2}(\Phi_\alpha, \tilde{\Phi})_G + 2\alpha \frac{d}{dt}(\Phi_\alpha, \tilde{\Phi})_G + (G_\alpha^* \Phi_\alpha, G_\alpha^* \tilde{\Phi})_\mathcal{H} = e^{-\alpha t}(Q_\alpha \Phi_b, \tilde{\Phi})_G + e^{-\alpha t}(G_\alpha^* \Phi_b, G_\alpha^* \tilde{\Phi})_\mathcal{H}. \quad (4.12) $$

Let $\Phi_0 = \Phi_\alpha - e^{-\alpha t} \Phi_b$ and define

$$ \begin{aligned} F_{\tilde{\Phi}} &= -\frac{d^2}{dt^2}(e^{-\alpha t} \tilde{\Phi}_b) - 2\alpha \frac{d}{dt}(e^{-\alpha t} \tilde{\Phi}_b) + e^{-\alpha t} Q_\alpha \tilde{\Phi}_b \\ &= e^{-\alpha t} \left( \alpha^2 \tilde{\Phi}_b - \frac{d^2}{dt^2} \tilde{\Phi}_b + Q_\alpha \tilde{\Phi}_b \right) \in L^2(0, T; D(G_\alpha^*)) \end{aligned} $$

then $\Phi_0$ is the unique solution to

$$ \frac{d^2}{dt^2}(\Phi_0, \tilde{\Phi})_G + 2\alpha \frac{d}{dt}(\Phi_0, \tilde{\Phi})_G + (G_\alpha^* \Phi_0, G_\alpha^* \tilde{\Phi})_\mathcal{H} = (F_{\tilde{\Phi}}, \tilde{\Phi})_G, \quad \forall \tilde{\Phi} \in D(G_\alpha^*) $$

with vanishing initial data (since $u_b$ and $\frac{d}{dt}u_b$ vanish at the initial time). From Theorem 4.3, the function $U_{\alpha,0}$ defined by $U_{\alpha,0} = G_\alpha^* \Phi_0 \in D(G_\alpha)$ is solution to

$$ \frac{d^2}{dt^2}(U_{\alpha,0}, \tilde{U})_\mathcal{H} + 2\alpha \frac{d}{dt}(U_{\alpha,0}, \tilde{U})_\mathcal{H} + (G_\alpha U_{\alpha,0}, G_\alpha \tilde{U})_G = (G_\alpha^* F_{\tilde{\Phi}}, \tilde{U})_\mathcal{H}, $$

where, by construction,

$$ G_\alpha^* F_{\tilde{\Phi}} = e^{-\alpha t} \left( \frac{d^2}{dt^2} L_\alpha(u_b) - \alpha^2 L_\alpha(u_b) \right). $$

This shows that $U_{\alpha,0} = G_\alpha^* \Phi_0$ is indeed the solution to (3.26). It remains only to define, as in the proof of Theorem 3.8, the function $U_\alpha = U_{\alpha,0} + e^{-\alpha t} L_\alpha(u_b) = U_{\alpha,0} - e^{-\alpha t} G_\alpha^* \Phi_b$ and observe that it is solution to (3.13) with the required regularity. $\square$

We have shown that for problems with artificial dissipation, the solution $U_\alpha$ to the velocity-field problem (3.13)-(3.15) is obtained from the solution $\Phi_\alpha$ to the potential-based problem (4.11)-(4.6), thanks to the relation $U_\alpha = G_\alpha^* \Phi_\alpha$. The next step is to prove a similar relation for the dissipative-free problems (3.5)-(3.7) and (2.8)-(2.9).

It should be noted that relating the dissipative potential-based problem (4.11) and its non-dissipative version (2.8) is not as direct as for the velocity-field problems. In particular, there is no relation of the form $\Phi(t) = e^{\alpha t} \Phi_\alpha(t)$, where $\Phi$ is the solution to (2.8) and $\Phi_\alpha$ is the solution to (4.11). Deducing $\Phi$ from $\Phi_\alpha$ requires to study the limit process $\alpha \to 0$.

## 4.2. The dissipation-free case

At this point we have shown that $U$, the solution to (3.5)-(3.7), can be expressed as

$$ U = e^{\alpha t} U_{\alpha} = e^{\alpha t} G_{\alpha}^{*} \Phi_{\alpha}. $$

The main idea is now to pass to the limit when $\alpha$ goes to zero in the equality above, to conclude that $U = G^* \Phi$, where $\Phi$ is solution to (2.8)-(2.9). Note that the weak convergence

$$ G_{\alpha}^{*} \Phi_{\alpha} \xrightarrow{L^2(0,T; \mathcal{H})} G^{*} \Phi $$
Unless stated otherwise, we use the simplified case where $N$ and $c_0$ are chosen constant. This gives an ordinary differential equation satisfied by $\rho_0$,

$$ \frac{d\rho_0}{dz} + \rho_0 \left( \frac{N^2}{g} + \frac{g}{c_0^2} \right) = 0, \qquad (5.2) $$

hence the density profile has the form

$$ \rho_0(z) = \rho_0(0) \exp(-n^2 z), \quad n^2 = \frac{N^2}{g} + \frac{g}{c_0^2} > 0. \qquad (5.3) $$

The numerical values used in the simulations are $\rho_0 = 1000 \text{ kg m}^{-3}$, $c_0 = 1500 \text{ ms}^{-1}$ and $N = 0.001 \text{ s}^{-1}$. Note that the model can also handle more complex stratification, where the background functions are computed from a given temperature profile. This yields depth-dependent $N$ and $c_0$. For conciseness, the case with a given temperature profile is presented in Appendix 7.1.

## 5.1 Discretization

In this part, we introduce the variational formulations to be discretized and the notations for the finite element approximation. The space discretization, obtained with a high-order spectral element method [23, 8] is then written. We end this part by presenting the fully-discrete scheme, which uses a second-order discretization in time.

### 5.1.1 velocity-field formulation

For the continuous problem (3.5)-(3.7), the velocity is sought in a subspace of $H(\text{div}, \Omega)$. The discretization of $H(\text{div}, \Omega)$ can be done using e.g. Raviart-Thomas elements [4], but the problem is discretized in $V_h \subset H^1(\Omega)$ for the simplicity of implementation. Such a strategy is certainly not adequate for harmonic problems [3], or for transient problems with a mean flow [2]. However, for the current problem we have not observed poor behaviour of the solution with this choice. To impose the boundary condition we use a Lagrange multiplier, hence, a finite dimensional space $M_h \subset L^2(\Omega)$ is introduced. The following discrete variational formulation is considered: find

$$ \mathbf{U}_h \in C^2([0, T]; V_h^d) \quad \text{and} \quad \lambda_h \in C^0([0, T]; M_h) \qquad (5.4) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_h(t), \tilde{\mathbf{U}}_h)_H + (\tilde{G}\mathbf{U}_h(t), \tilde{G}\tilde{\mathbf{U}}_h)_G + (\lambda_h, \tilde{\mathbf{U}}_h \cdot \mathbf{n}_b)_{L^2(\Gamma_b)} = 0, \quad \forall \tilde{\mathbf{U}}_h \in V_h^d, \quad \forall t \in [0, T], \qquad (5.5) $$

$$ (\mathbf{U}_h(t) \cdot \mathbf{n}_b - u_b, \mu_h)_{L^2(\Gamma_h)} = 0 \quad \forall \mu_h \in M_h, \quad \forall t \in [0, T]. \qquad (5.6) $$

The spectral element method is used for the construction of $V_h$, and $M_h$ is constructed using traces of functions in $V_h$. Adequate quadrature formulae (that use Gauss-Lobatto quadrature rule) are used to compute integrals, leading in particular to mass lumping: the mass matrix is diagonal. The semi-discrete algebraic problem reads

$$ \frac{d^2}{dt^2} \mathbb{M}_U U_h + \mathbb{K}_U U_h + C_h \Lambda_h = 0, \qquad (5.7) $$

$$ C_h^T U_h - U_{b,h} = 0, \qquad (5.8) $$

where $U_h$ and $\Lambda_h$ are vectors of degrees of freedoms. The time interval $[0, T]$ is partitioned into $M$ equal intervals of length $\Delta t = T/M$. The finite-element approximation (5.7) is discretized in time with a leapfrog scheme. We consider the sequence $\{U_h^n \in V_h\}_{n \in \{1, \dots, M\}}$ solution to

$$ \mathbb{M}_U \frac{U_h^{n+1} - 2U_h^n + U_h^{n-1}}{\Delta t^2} + \mathbb{K}_U U_h^n + C_h \Lambda_h^n = 0, \qquad (5.9) $$

$$ C_h^T U_h^n - U_{b,h}(n\Delta t) = 0. \qquad (5.10) $$

Note that in (5.9) the mass matrix $\mathbb{M}_U$ is diagonal, hence easily invertible. Moreover, at each time step, the matrix $C_h^T C_h$ must also be inverted. It is also diagonal when adequate quadrature formulae are used and thanks to our choice of having $M_h$ as a space of traces of functions of $V_h$.
would imply that $\mathbf{U} = G^*\Phi$. It is then sufficient to prove that the weak convergence holds and that the limit $\Phi$ is solution to (2.8)-(2.9).

To prove the convergence, we start by giving preliminary energy estimates for the solution $\Phi_\alpha$. Those estimates are uniform with respect to $\alpha$. Following the same proof as Proposition 2.2, we have

$$ \frac{1}{2} \| \frac{d}{dt} \Phi_{\alpha}(t) \|_{\mathcal{G}}^{2} + 2\alpha \int_{0}^{t} \| \Phi_{\alpha}(t') \|_{\mathcal{G}}^{2} dt' + \frac{1}{2} \| G_{\alpha}^{*} \Phi_{\alpha}(t) \|_{\mathcal{H}}^{2} \lesssim B_{\alpha}^{2}(t), \quad (4.13) $$

up to a multiplicative constant – independent of $\Phi_\alpha$ and of the source term $u_b$ – and where $B_\alpha$ is defined by

$$ B_{\alpha}(t) = \sup_{s \in [0, t]} \|e^{-\alpha s} u_b(t)\|_{H^{-1/2}(\Gamma_b)} + \int_0^t \|\partial_t(e^{-\alpha s} u_b(s))\|_{H^{-1/2}(\Gamma_b)} \, ds. $$

For $\alpha \le 1$, we have

$$ B_{\alpha}(t) \le D(t) := C(t) + \int_{0}^{t} \|u_{b}(s)\|_{H^{-1/2}(\Gamma_{b})} \, ds, $$

where $B(t)$ is defined in (2.15). Therefore, the right-hand side of (4.13) can be replaced by a positive function independent of $\alpha$. For $\alpha \le 1$, we have

$$ \| \frac{d}{dt} \Phi_{\alpha}(t) \|_{\mathcal{G}} + \| G_{\alpha}^{*} \Phi_{\alpha}(t) \|_{\mathcal{H}} \lesssim C(t). \quad (4.14) $$

From (4.11), we also deduce an estimate of the second order time derivative in the Hilbert space $\mathcal{D}(G_\alpha^*)'$,

$$ \| \frac{d^2}{dt^2} \Phi_{\alpha}(t) \|_{D(G_{\alpha}^{*})'} \lesssim C(t). \quad (4.15) $$

Those preliminary observations allow us to state the following result.

**Lemma 4.5.** The functions $G_\alpha^* \Phi_\alpha$ for $\alpha \in \mathbb{R}^+$ converges weakly in $L^2(0, T; \mathcal{H})$ to $G^*\Phi$ when $\alpha \to 0$ with $\Phi$ solution to (2.8)-(2.9).

*Proof.* Using the estimates (4.14)-(4.15) and the fact that $\Phi_\alpha$ vanishes at the initial time, we have that $\Phi_\alpha$, when considered as a sequence in $\alpha$, is bounded in the Hilbert space

$$ \mathcal{W} := L^2(0, T; \mathcal{D}(G_\alpha^*)) $$

Note that the domain of the operator $G_\alpha^*$ is independent of $\alpha$, from the definition $\mathcal{D}(G_\alpha^*) = \mathcal{D}(G^*) \times \mathcal{H}$. Hence, we deduce that $\Phi_\alpha$ converges weakly – up to a subsequence – to a function $\tilde{\Phi}$ in $\mathcal{W}$. Decomposing $\tilde{\Phi} = (\tilde{\Phi} \mathbf{V})^\top$ and passing to the limit, $\alpha \to 0$ in the formulation (4.11), we obtain

$$ \frac{d^2}{dt^2}(\Phi, \tilde{\Phi})_G + (G^*\Phi, G^*\tilde{\Phi})_{\mathcal{H}} = \langle u_b, \tilde{\varphi} \rangle_{\Gamma_b} \quad \forall \tilde{\Phi} = (\Phi 0)^\top \in \mathcal{D}(G^*) \times \mathcal{H}, \quad \text{in } \mathcal{D}'(0, T), \quad (4.16) $$

and

$$ \frac{d^2}{dt^2}(\mathbf{V}, \tilde{\mathbf{V}})_{\mathcal{H}} = 0 \quad \forall \tilde{\mathbf{\Phi}} = (0 \ \tilde{\mathbf{V}})^\top \in \mathcal{D}(G^*) \times \mathcal{H}, \quad \text{in } \mathcal{D}'(0, T). \quad (4.17) $$

Hence $\Phi$ is solution to (2.8). It only remains to investigate the initial conditions to conclude the proof. From the energy estimate (4.14) we also have the weak convergence in $H^1(0, T; \mathcal{G})$, therefore, for all $\tilde{\Phi} \in \mathcal{G}$,

$$ (\Phi_{\alpha}(0), \tilde{\Phi})_{\mathcal{G}} = \frac{1}{T} \int_{0}^{T} \frac{d}{dt} ((t-T)(\Phi_{\alpha}, \tilde{\Phi})_{\mathcal{G}}) \, dt \underset{\alpha \to 0}{\longrightarrow} \frac{1}{T} \int_{0}^{T} \frac{d}{dt} ((t-T)(\Phi, \tilde{\Phi})_{\mathcal{G}}) \, dt = (\Phi(0), \tilde{\Phi})_{\mathcal{G}}, $$
### 5.1.2 Potential-based formulation

As presented in Section 2, we use the variational formulation (2.6). The problem is written as a coupled system for the variables $(\varphi(t), \psi(t)) \in H^1(\Omega) \times L^2(\Omega)$ for all time $t \in [0, T]$. By taking test function $\tilde{\psi} = 0$ in (2.6), the variational formulation reads

$$ \frac{d^2}{dt^2} \left( \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, dx + \int_{\Gamma_s} \frac{\rho_0}{g} \varphi \tilde{\varphi} \, dx + \int_{\Omega} \rho_0 \left( \nabla \varphi - \frac{N^2}{g} \varphi \mathbf{e}_z \right) \cdot \left( \nabla \tilde{\varphi} - \frac{N^2}{g} \tilde{\varphi} \mathbf{e}_z \right) \, dx \right. \\ \left. - \int_{\Omega} \rho_0 N \psi \mathbf{e}_z \cdot \left( \nabla \tilde{\varphi} - \frac{N^2}{g} \tilde{\varphi} \mathbf{e}_z \right) \, dx = - \int_{\Gamma_b} \frac{\rho_0}{g} u_b \tilde{\varphi}, \quad (5.11) $$

and by taking test function $\tilde{\varphi} = 0$ in (2.6), the variational formulation reads

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \psi \tilde{\psi} \, dx - \int_{\Omega} \rho_0 \left( \nabla \varphi - \frac{N^2}{g} \varphi \mathbf{e}_z \right) \cdot \left( N \tilde{\psi} \mathbf{e}_z \right) \, dx + \int_{\Omega} \rho_0 N^2 \psi \tilde{\psi} \, dx = 0. \quad (5.12) $$

We introduce a finite-dimensional space $L_h \subset L^2(\Omega)$. The following bilinear forms are defined for $\varphi_h, \tilde{\varphi}_h \in V_h$ and $\psi_h, \tilde{\psi}_h \in L_h$,

$$ m_{\varphi}(\varphi_h, \tilde{\varphi}_h) = \left(\frac{\rho_0}{c_0^2} \varphi_h, \tilde{\varphi}_h\right)_{L^2(\Omega)} + \left(\frac{\rho_0}{g} \varphi_h, \tilde{\varphi}_h\right)_{L^2(\Gamma_s)}, \quad m_{\psi}(\psi_h, \tilde{\psi}_h) = \left(\rho_0 \psi_h, \tilde{\psi}_h\right)_{L^2(\Omega)}, \quad (5.13) $$

$$ a_{\varphi}(\varphi_h, \varphi_h) = (\nabla \varphi_h - N^2 \varphi_h \mathbf{e}_z, \nabla \varphi_h - N^2 \varphi_h \mathbf{e}_z)_{\mathcal{H}}, \quad a_{\psi}(\psi_h, \psi_h) = (\rho_0 N^2 \psi_h, \tilde{\psi}_h)_{L^2(\Omega)}, \quad (5.14) $$

$$ c(\psi_h, \tilde{\psi}_h) = (-\rho_0 N \psi_h \mathbf{e}_z, \nabla \tilde{\psi}_h - N \tilde{\psi}_h \mathbf{e}_z)_{\mathcal{H}}. \quad (5.15) $$

The semi-discrete variational formulation reads then: find

$$ (\varphi_h, \psi_h) \in C^2([0, T]; V_h \times L_h) \quad (5.16) $$

solution to

$$ \frac{d^2}{dt^2} m_{\varphi}(\varphi_h, \tilde{\varphi}_h) + a_{\varphi}(\varphi_h, \tilde{\varphi}_h) + c(\psi_h, \tilde{\psi}_h) = -(\rho_0 u_b / g, \tilde{\varphi}_h)_{L^2(\Gamma_b)}, \quad \forall \tilde{\varphi}_h \in V_h, \quad \forall t \in [0, T], \quad (5.17) $$

$$ \frac{d^2}{dt^2} m_{\psi}(\psi_h, \tilde{\psi}_h) + a_{\psi}(\psi_h, \tilde{\psi}_h) + c(\tilde{\psi}_h, \varphi_h) = 0, \quad \forall \tilde{\psi}_h \in L_h, \quad \forall t \in [0, T], \quad (5.18) $$

An approximated velocity $\mathbf{U}_h$ is computed from the variables $(\varphi_h, \psi_h)$ using formally the relation $\mathbf{U} = G^*\Phi$. At the discrete level we solve

$$ (\mathbf{U}_h, \tilde{\mathbf{U}}_h)_{\mathcal{H}} = b_{\varphi}(\varphi_h, \tilde{\mathbf{U}}_h) + b_{\psi}(\psi_h, \tilde{\mathbf{U}}_h), \quad \forall V_h^d \in V_h, \quad (5.19) $$

where $b_{\varphi}: V_h \times V_h^d \to \mathbb{R}$ and $b_{\psi,h}: L_h \times V_h^d \to \mathbb{R}$ are the following bilinear forms,

$$ b_{\varphi}(\varphi_h, \tilde{\mathbf{U}}_h) = -(\nabla \varphi_h - \frac{N^2}{g} \mathbf{e}_z \varphi_h, \tilde{\mathbf{U}}_h)_{\mathcal{H}}, \quad b_{\psi,h}(\psi_h, \tilde{\mathbf{U}}_h) = (N \psi_h \mathbf{e}_z, \tilde{\mathbf{U}}_h)_{\mathcal{H}}. \quad (5.20) $$

In practice the spectral element method is used to construct $V_h$ and $L_h$ and all integrals are approximated using an adequate quadrature formula. For the potential-based formulation, the decomposition on each basis yields the global vectors $\Phi_h, \Psi_h$ and $u_{b,h}$, as well as the mass matrices $M_{\varphi}, M_{\psi}$ and $M_b$, the stiffness matrices $K_{\varphi}, K_{\psi}$, and the interaction matrix $C_{\psi,\varphi}$. We also obtain the matrices $B_{\varphi}, B_{\psi}$ used to deduce the velocity from the potential, see Equation (5.19). The semi-discrete system corresponding to the finite-element approximation (5.17)-(5.19) reads

$$ \frac{d^2}{dt^2} M_{\varphi} \Phi_h + K_{\varphi} \Phi_h + C_{\psi,\varphi} \Psi_h = -M_b U_{b,h}, \quad (5.21) $$

$$ \frac{d^2}{dt^2} M_{\psi} \Psi_h + K_{\psi} \Psi_h + C_{\psi,\varphi}^{T} \Phi_h = 0, \quad (5.22) $$

$$ M_{\mathcal{U}} U_h = B_{\varphi} \Phi_h + B_{\psi} \Psi_h. \quad (5.23) $$
hence $\Phi(0) = \Phi_{\alpha}(0) = 0$. Similarly, thanks to the weak convergence in $H^2(0, T; (\mathcal{D}(G^*) \times \mathcal{H})')$ we conclude that

$$ \frac{d}{dt} \Phi(0) = \frac{d}{dt} \Phi_{\alpha}(0) = 0. $$

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

We have shown that $\Phi_{\alpha}$ converges weakly to $\Phi$, that the limit $\Phi$ is solution to (2.8)-(2.9), and that $G_{\alpha}^* \Phi_{\alpha}$ converges weakly to $G^* \Phi$. Hence, $U = G^* \Phi$ is solution to (3.5)-(3.7); this proves Theorem 4.1.

## 5. NUMERICAL ILLUSTRATIONS

We now use the two formulations (3.5)-(3.7) and (2.6) to give several numerical illustrations of tsunamis and hydro-acoustic waves generated by an earthquake or landslide source. Three sets of simulations are presented in this section. The first one reproduces a classical scenario of a submarine earthquake generating a tsunami and acoustic waves. This scenario serves several purposes: convergence analysis, model validation by comparison with the literature, and illustration of the equivalence between the potential-based and velocity-field formulations. In the second set of simulations, the equation $U = G^* \Phi$ is used to investigate the validity of the classical irrotational assumption for the velocity field. The third simulation is a preliminary work towards the study of acoustic waves generated by submarine landslides. We compute the spectrogram of pressure recorded at a given sensor for a source located on the seabed which emits a range of frequencies typical of those observed in the field. This novel simulation highlights the interference pattern caused by the reflection of the acoustic waves at the surface.

We first describe the choice of parameters. The discretization of the velocity-field and the potential-based simulations are then presented. Finally, for each set of simulations, we describe the scenario and provide some illustrations.

As stated in the introduction, the background functions $\rho_0, c_0$ must satisfy the positivity property

$$ -\frac{g}{\rho_0} \frac{d\rho_0}{dz} - \frac{g^2}{c_0^2} = N^2 \geq 0. \qquad (5.1) $$

Unless stated otherwise, we use the simplified case where $N$ and $c_0$ are chosen constant. This gives an ordinary differential equation satisfied by $\rho_0$,

$$ \frac{d\rho_0}{dz} + \rho_0 \left( \frac{N^2}{g} + \frac{g}{c_0^2} \right) = 0, \qquad (5.2) $$

hence the density profile has the form

$$ \rho_0(z) = \rho_0(0) \exp(-n^2 z), \quad n^2 = \frac{N^2}{g} + \frac{g}{c_0^2} > 0. \qquad (5.3) $$

The numerical values used in the simulations are $\rho_0(0) = 1000 \text{ kg m}^{-3}$, $c_0 = 1500 \text{ ms}^{-1}$, $N = 0.001 \text{ s}^{-1}$ and $g = 9.81 \text{ ms}^{-2}$. Note that the model can also handle more complex stratification, where the background functions are computed from a given temperature profile. This yields depth-dependent $N$ and $c_0$. For conciseness, the case with a given temperature profile is presented in Appendix 7.2.

### 5.1. Discretization

In this part, we introduce the variational formulations, write the space discretization obtained with a high-order spectral element method, and describe the fully-discrete scheme.
<figure><img src="image_8.png" /><figcaption>Figure 2: The function $u_b(x, t)$ describing the earthquake used in the first simulation to generate acoustic waves and tsunamis.</figcaption></figure>

By using the Gauss-Lobatto quadrature rule for the integral approximations, we obtain that the mass matrices $M_\varphi$, $M_\psi$ and $M_b$ are diagonal. Again, the time interval $[0, T]$ is partitioned into $M$ equal intervals of length $\Delta t = T/M$ and for the potential-based formulation, we consider the sequences $\{\Phi_h^n \in V_h\}_{n \in \{1, \dots, M\}}$ and $\{\Psi_h^n \in L_h\}_{n \in \{1, \dots, M\}}$ solution to

$$M_\varphi \frac{\Phi_h^{n+1} - 2\Phi_h^n + \Phi_h^{n-1}}{\Delta t^2} + K_\varphi \Phi_h^n + C_{\varphi\psi} \Psi_h^n = -M_b U_{b,h}(n\Delta t), \quad (5.24)$$

$$M_\psi \frac{\psi_h^{n+1} - 2\psi_h^n + \psi_h^{n-1}}{\Delta t^2} + K_\psi \psi_h^n + C_{\psi\psi}^t \Phi_h^n = 0, \quad (5.25)$$

and the velocity is then computed with

$$M_b U_h^n = B_\varphi \Phi_h^n + B_\psi \psi_h^n. \quad (5.26)$$

In the following, the displacement will be sometimes used. The displacement is computed by integrating in time the velocity.

## 5.2 First simulation: waves generation from a submarine earthquake

In order to validate and compare the formulations (5.9) and (5.24)-(5.26), we reproduce the test case presented in [27]. A tsunami and hydro-acoustic waves are generated by a seabed elevation due to a submarine earthquake, in a 2D domain on a flat topography. In [27], the flow is assumed irrotational with a velocity potential $\phi$. The distance between the seabed and the mean water level is denoted $h(x, y, t)$. The model, solved with a finite elements method, reads

$$\frac{\partial^2 \phi}{\partial t^2} - c^2 \Delta \phi = 0,$$

where $c = 1500 \text{ ms}^{-1}$ is constant, and with boundary conditions (with $g = 9.81 \text{ ms}^{-2}$)

$$\frac{\partial^2 \phi}{\partial t^2} + g \frac{\partial \phi}{\partial z} = 0, \quad \text{at } z = h(x, y, t), \quad \nabla \phi \cdot \mathbf{n} = \frac{\partial h}{\partial t}, \quad \text{at } z = 0.$$

In the simulations, we use the following source: $u_b(x, t) = f(x)g(t)$, where $f$ and $g$ are smoothed rectangular functions,

$$f(x) = \frac{1}{1 + e^{-s_x(x-r_x/2)}} - \frac{1}{1 + e^{-s_x(x+r_x/2)}}, \quad g(t) = \frac{1}{1 + e^{-s_t(t-t_0)}} - \frac{1}{1 + e^{-s_t(t-t_0-r_t)}}. \quad (5.27)$$

see Figure 2. Since no conditions are imposed on the lateral boundaries, the computational domain is wider to avoid reflections inside the domain of interest. The numerical values are given in Table 1.
### 5.1.1. Velocity-field formulation

For the continuous problem (3.5)-(3.7), the velocity is sought in a subspace of $H(\text{div}, \Omega)$. The discretization of $H(\text{div}, \Omega)$ can be done using e.g. Raviart-Thomas elements [4]. Here, for the simplicity of implementation, the problem is discretized in $V_h \subset H^1(\Omega)$. Such a strategy is certainly not adequate for harmonic problems [3], or for transient problems with a mean flow [2]. However, for the current problem we have not observed poor behaviour of the solution with this choice. The essential boundary condition is imposed with a Lagrange multiplier, and we denote by $M_h \subset L^2(\Gamma_b)$ its finite-dimensional space. The following discrete variational formulation is considered: Find

$$ \mathbf{U}_h \in C^2([0, T]; V_h^d) \quad \text{and} \quad \lambda_h \in C^0([0, T]; M_h) \qquad (5.4) $$

solution to

$$ \frac{d^2}{dt^2}(\mathbf{U}_h(t), \tilde{\mathbf{U}}_h)_H + (\tilde{G}\mathbf{U}_h(t), \tilde{G}\tilde{\mathbf{U}}_h)_G + (\lambda_h, \tilde{\mathbf{U}}_h \cdot \mathbf{n}_b)_{L^2(\Gamma_b)} = 0, \quad \forall \tilde{\mathbf{U}}_h \in V_h^d, \quad \forall t \in [0, T], \qquad (5.5) $$

$$ (\mathbf{U}_h(t) \cdot \mathbf{n}_b - u_b, \mu_h)_{L^2(\Gamma_b)} = 0 \quad \forall \mu_h \in M_h, \quad \forall t \in [0, T]. \qquad (5.6) $$

The space $M_h$ is constructed using traces of functions in $V_h$. For the spaces $M_h$, $V_h$ and for the discretization of the variational formulation, we use the spectral element method [23]. The spectral element method it a type of finite element method that uses Gauss-Lobatto quadrature rule to compute the integrals of the variational formulation [8]. The main benefit of this method is the mass-lumping: the resulting mass matrices are diagonal. After space discretization, the semi-discrete algebraic problem reads

$$ \frac{d^2}{dt^2} \mathbb{M} U_h + \mathbb{K} U_h + \mathbb{C} \Lambda_h = 0, \qquad (5.7) $$

$$ \mathbb{C}^{\top} U_h - \mathbb{M}_b U_{b,h} = 0, \qquad (5.8) $$

where $U_h$ and $\Lambda_h$ denote the time-dependent vectors of degrees of freedom and $U_{b,h}$ corresponds to the degree of freedoms of the decomposition of $u_b$ in the space $M_h$. The finite-element approximation (5.7) is then discretized in time with a leapfrog scheme. The time interval $[0, T]$ is partitioned into equal intervals of length $\Delta t$. The time step $\Delta t$ will depend on a CFL condition, which is described later. We consider the sequences $\{U_h^n\}_{n \in \{1, \dots, M\}} \subset V_h^d$ and $\{\Lambda_h^n\}_{n \in \{1, \dots, M\}} \subset M_h$ solution to

$$ \mathbb{M} \frac{U_h^{n+1} - 2U_h^n + U_h^{n-1}}{\Delta t^2} + \mathbb{K} U_h^n + \mathbb{C} \Lambda_h^n = 0, \qquad (5.9) $$

$$ \mathbb{C}^{\top} U_h^n - U_{b,h}^n = 0. \qquad (5.10) $$

The system is then solved as follows: given $\mathbf{U}_h^n, \Lambda_h^n$ from the previous time step, the solution $\mathbf{U}_h^{n+1}, \Lambda_h^{n+1}$ is obtained with

$$ U_h^{n+1} = 2U_h^n - U_h^{n-1} - \Delta t^2 \mathbb{M}^{-1} (\mathbb{K} U_h^n + \mathbb{C} \Lambda_h^n), \qquad (5.11) $$

$$ \Lambda_h^{n+1} = -(\mathbb{C}^{\top}\mathbb{M}^{-1}\mathbb{C})^{-1} \left( \mathbb{M}_b \frac{U_{b,h}^{n+2} - 2U_{b,h}^{n+1} + U_{b,h}^n}{\Delta t^2} + \mathbb{C}^{\top}\mathbb{M}^{-1}\mathbb{K} U_h^{n+1} \right). \qquad (5.12) $$

The mass matrix $\mathbb{M}$ is diagonal, hence easily inverted. At each time step, the matrix $\mathbb{C}^{\top}\mathbb{M}^{-1}\mathbb{C}$ must also be inverted. It is also diagonal when adequate quadrature formulae are used thanks to our choice of having $M_h$ as a space of traces of functions of $V_h$. Moreover, using that for all $\mu_h \in M_h$,

$$ (\tilde{\mathbf{U}}_h \cdot \mathbf{n}_b, \mu_h)_{L^2(\Gamma_b)} = 0 \quad \forall \tilde{\mathbf{U}}_h \in V_h^d \Rightarrow \mu_h = 0, $$
<table><caption>Table 1: Parameter values for the first simulation: hydro-acoustic waves and tsunami generation by an earthquake.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain height</td><td>1.5 km</td></tr><tr><td>s<sub>x</sub></td><td>150 m<sup>-1</sup></td></tr><tr><td>r<sub>x</sub></td><td>30 km</td></tr><tr><td>t<sub>0</sub></td><td>2 s</td></tr><tr><td>s<sub>t</sub></td><td>4 s<sup>-2</sup></td></tr><tr><td>r<sub>t</sub></td><td>1 s</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 3: Surface vertical displacement at $x = 50$ km for the velocity formulation at different element orders.</figcaption></figure>

### 5.2.1 Convergence analysis

For each formulation, several element orders are tested. The aim is to show that the simulations converge and to select parameters that are precise enough and not too costly. The tsunami is simulated for 50 s and the computational domain is 101 km long. To evaluate whether convergence is reached, we plot the time evolution of the free surface displacement at a point $x = 50$ km away from the source center. The convergence is estimated visually and is assumed reached when increasing the order does not change the curves. The finite element orders are denoted $P_x$, $P_z$ and the number of mesh subdivision in the $x$ and in the $z$ coordinates are respectively denoted $N_x$, $N_z$. Several combinations of $P_x$ and $P_z$ are tested while keeping $N_x$ and $N_z$ fixed. The total number of Degrees of Freedom (NDoFs) for each case is indicated in each surface displacement plot.

For the velocity-field scheme, the surface displacement is shown in Figure 3. The convergence is reached for $P_x = 4$, $P_z = 5$. For the potential-based scheme, the surface displacement is shown in Figure 4a and we see that the orders with $P_x \ge 3$ converge. However, when taking $P_x = 3$, $P_z = 3$, we have observed increasing oscillations for longer times (not shown here). These oscillations come from the seabed displacement at $x = 0$ km, as shown in Figure 4b. The convergence is actually reached for $P_x = 3$, $P_z = 5$. Such oscillations are not observed for the velocity-field formulation.

### 5.2.2 Comparison and snapshots

We now run the simulation for a longer time and compare it with the results from [27]. The simulation obtained with the velocity-field formulation and the potential-based formulations are first compared with [27], then compared to each other. Finally we present some snapshots of the simulation.

The tsunami is simulated for 1000 s, and the computational domain is 1576.5 km long. The element orders, the DoFs and the computational time are given in Table 2 and 3 respectively for the velocity-field and potential-based formulations. The potential-based formulation requires more computational time than the velocity-field formulation. This additional time is due to the displacement calculation: the displacement is obtained by in-
and using that this property also holds with adequate quadrature formulae, we are able to show that $\text{Ker } \mathbb{C} = \{0\}$. It follows that $\mathbb{C}^\top \mathbb{M}^{-1} \mathbb{C}$ is a positive definite matrix.

### 5.1.2. Potential-based formulation

As presented in Section 2, we use the variational formulation (2.6). The problem is written as a coupled system for the variables $(\varphi(t), \psi(t)) \in H^1(\Omega) \times L^2(\Omega)$ for all time $t \in [0, T]$. By taking test function $\tilde{\psi} = 0$ in (2.6), the variational formulation reads

$$ \frac{d^2}{dt^2} \left( \int_{\Omega} \frac{\rho_0}{c_0^2} \varphi \tilde{\varphi} \, dx + \int_{\Gamma_s} \frac{\rho_0}{g} \varphi \tilde{\varphi} \, dx \right) + \int_{\Omega} \rho_0 \left( \nabla \varphi - \frac{N^2}{g} \varphi \mathbf{e}_z \right) \cdot \left( \nabla \tilde{\varphi} - \frac{N^2}{g} \tilde{\varphi} \mathbf{e}_z \right) \, dx \\ - \int_{\Omega} \rho_0 N \psi \mathbf{e}_z \cdot \left( \nabla \tilde{\varphi} - \frac{N^2}{g} \tilde{\varphi} \mathbf{e}_z \right) \, dx = - \int_{\Gamma_b} \frac{\rho_0}{g} u_b \tilde{\varphi}, \quad (5.13) $$

and by taking test function $\tilde{\varphi} = 0$ in (2.6), the variational formulation reads

$$ \frac{d^2}{dt^2} \int_{\Omega} \rho_0 \psi \tilde{\psi} \, dx - \int_{\Omega} \rho_0 \left( \nabla \varphi - \frac{N^2}{g} \varphi \mathbf{e}_z \right) \cdot \left( N \tilde{\psi} \mathbf{e}_z \right) \, dx + \int_{\Omega} \rho_0 N^2 \psi \tilde{\psi} \, dx = 0. \quad (5.14) $$

Moreover, the velocity $\mathbf{U}_h$ is computed from the variables $(\varphi_h, \psi_h)$ using the relation $\mathbf{U} = G^* \Phi$. The weak formulation is

$$ \int_{\Omega} \mathbf{U} \cdot \tilde{\mathbf{U}} \, dx = \int_{\Omega} \left( \nabla \varphi - \frac{N^2}{g} \mathbf{e}_z \varphi \right) \cdot \tilde{\mathbf{U}} \, dx + \int_{\Omega} N \psi \mathbf{e}_z \cdot \tilde{\mathbf{U}} \, dx. \quad (5.15) $$

We introduce the discrete space $L_h \subset L^2(\Omega)$. The space discretisation of the system (5.13)-(5.15) with the spectral element method yields the semi-discrete system

$$ \frac{d^2}{dt^2} \mathbb{M}_{\varphi} \Phi_h + \mathbb{K}_{\varphi} \Phi_h + \mathbb{C}_{\psi \varphi} \Psi_h = - \mathbb{M}_b U_{b,h}, \quad (5.16) $$

$$ \frac{d^2}{dt^2} \mathbb{M}_{\psi} \Psi_h + \mathbb{K}_{\psi} \Psi_h + \mathbb{C}_{\psi \varphi}^{\top} \Phi_h = 0, \quad (5.17) $$

$$ \mathbb{M}_{\mathbf{U}} U_h = \mathbb{G}_{\varphi} \Phi_h + \mathbb{G}_{\psi} \Psi_h, \quad (5.18) $$

Where $\Phi_h \in V_h$, $\Psi_h \in L_h$ and $U_h \in V_h^d$ denote the vectors of degrees of freedom for $\varphi, \psi$ and $\mathbf{U}$, respectively, and where $U_b$ is the discretization of the right-hand side of (5.13). As mentioned in the previous section, the Gauss-Lobatto quadrature rule ensures that the mass matrices $\mathbb{M}_{\varphi}, \mathbb{M}_{\psi}$ and $\mathbb{M}_{\mathbf{U}}$ are diagonal. For the time discretization, the time interval $[0, T]$ is partitioned into equal intervals of length $\Delta t$. We consider the sequences $\{\Phi_h^n \in V_h\}_{n \in \{1, \dots, M\}}$ and $\{\Psi_h^n \in L_h\}_{n \in \{1, \dots, M\}}$ solution to

$$ \mathbb{M}_{\varphi} \frac{\Phi_h^{n+1} - 2\Phi_h^n + \Phi_h^{n-1}}{\Delta t^2} + \mathbb{K}_{\varphi} \Phi_h^n + \mathbb{C}_{\psi \varphi} \Psi_h^n = -U_{b,h}^n, \quad (5.19) $$

$$ \mathbb{M}_{\psi} \frac{\Psi_h^{n+1} - 2\Psi_h^n + \Psi_h^{n-1}}{\Delta t^2} + \mathbb{K}_{\psi} \Psi_h^n + \mathbb{C}_{\psi \varphi}^{\top} \Phi_h^n = 0, \quad (5.20) $$

and the velocity is then computed with

$$ \mathbb{M}_{\mathbf{U}} U_h^n = \mathbb{G}_{\varphi} \Phi_h^n + \mathbb{G}_{\psi} \Psi_h^n. \quad (5.21) $$

Finally, the discrete displacement, denoted by $d_h^n$, is computed by integrating the velocity in time:

$$ d_h^{n+1} = d_h^n + \Delta t U_h^n. \quad (5.22) $$
<figure><img src="image_10.png" /><figcaption>Figure 4: Displacement of (a) the water free surface at $x = 50$ km, (b) the seabed at $x = 0$ km, obtained from the potential-based formulation for several element orders.</figcaption></figure>

<table><caption>Table 2: Parameters values for the velocity-field simulation.</caption><thead><tr><td>Nx</td><td>Nz</td><td>Px</td><td>Pz</td><td>NDoFs</td><td>Computational time</td></tr></thead><tbody><tr><td>1051</td><td>10</td><td>4</td><td>5</td><td>428 910</td><td>771 s</td></tr></tbody></table>

<table><caption>Table 3: Parameters values for the potential-based simulation.</caption><thead><tr><td>Nx</td><td>Nz</td><td>Px</td><td>Pz</td><td>DoFs</td><td>Computational time</td></tr></thead><tbody><tr><td>1051</td><td>10</td><td>3</td><td>5</td><td>413 094</td><td>1197 s</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIGURE 2. The function $u_b(x, t)$ describing the earthquake used in the first simulation to generate acoustic waves and tsunamis.</figcaption></figure>

### 5.1.3. Mesh and time discretization

The 2D domain $\Omega$ is subdivided into quadrilaterals, with $Nx$ subdivisions in the $x$ direction and $Nz$ subdivisions in the $z$ direction. We use continuous quadrilateral $Q_{P_x,P_z}$ spectral finite elements [8, 23] for constructing $V_h$ ($L_h$ is constructed similarly but with discontinuous elements). These elements are the tensor product of 1D $Q_{P_x}$-elements with 1D $Q_{P_z}$-elements, respectively in the $x$ and $z$ direction. For the time discretization, since an explicit leap-frog scheme is used, a CFL condition must be satisfied [22]. It involves the spectral radius of the matrix $M^{-1}K$ (denoted $\rho(M^{-1}K)$) as follows,

$$ \Delta t \le \frac{2}{\sqrt{\rho(M^{-1}K)}} $$

Hence, the time step $\Delta t$ depends on the mesh size and on the finite element order. Its numerical value will be given for each simulation.

## 5.2. First simulation: waves generation from a submarine earthquake

In order to validate and compare the formulations (5.9) and (5.19)-(5.21), we reproduce the test case presented in [27]. Hydro-acoustic waves and a tsunami are generated by a submarine earthquake in a 2D domain with a flat seabed. In [27], the flow is assumed irrotational with a velocity potential $\phi$. The distance between the seabed and the mean water level is denoted by $h(x, y, t)$. The model, solved with a finite element method, reads

$$ \frac{\partial^2 \phi}{\partial t^2} - c_0^2 \Delta \phi = 0, $$

where $c_0 = 1500 \text{ ms}^{-1}$ is constant, and with boundary conditions (with $g = 9.81 \text{ ms}^{-2}$)

$$ \frac{\partial^2 \phi}{\partial t^2} + g \frac{\partial \phi}{\partial z} = 0, \quad \text{at } z = h(x, t), \quad \nabla \phi \cdot \mathbf{n} = \frac{\partial h}{\partial t}, \quad \text{at } z = 0. $$

It should be noted that in [27], only the width and velocity of the source are given, but not its exact expression.
<figure><img src="image_11.png" /><figcaption>Figure 5: Vertical displacement at x = 50 km, computed from the velocity-field (blue) superimposed to the result in [27].</figcaption></figure>

tegrating the velocity in time, and to obtain a precise integration, the velocity is computed at each time step. Hence the system (5.26) must be solved at each time step, leading to an additional computational cost. The code could be accelerated by using an efficient time integration rule which would not require to compute the velocity at each time step, or by directly solving the potential-based problem associated to the displacement.

For comparison with the literature, we consider the surface vertical displacement at a point x = 50 km evolving with time. The vertical displacement computed with the velocity-field model is superimposed to the result from [27] in Figure 5. Both curves have a high-frequency component starting at time t ~ 20 s, namely the acoustic waves, and a low-frequency wave starting at t ~ 300 s, which is the tsunami. The tsunami arrival times are in good agreement between the two curves. It should be noted that since we do not have the exact description of the boundary condition used in [27], the earthquake description slightly differs between the two simulations. This could explain the phase difference between both displacements.

To compare the potential-based and the velocity-field formulations, we plot in Figure 6a the same vertical displacement obtained from both formulations. The orange curve is the displacement $d_\phi$ computed from the potential formulation, and the blue curve is the displacement $d_u$ computed from the velocity formulation. The curves are almost superimposed, which illustrates the fact that both formulations are equivalent. Figure 6b shows the error $|d_\phi - d_u|$ between the potential-based and the velocity-field formulations. The error has a magnitude of less than 5 percents, and it oscillates with a period of approximately 4 seconds. The error averaged over a period is plotted in the same figure. We see that the averaged error slowly increases, which could be caused by the time integration: small errors add up and become non-negligible for large times.

Finally, to give some new illustration of the wave propagation, we show in Figure 7 snapshots of the Eulerian domain near the earthquake at several times. The Eulerian domain is obtained by deforming the domain at rest with the displacement vector field. The colors correspond to the values of the vertical displacement. On each snapshot, we see the permanent deformation of the seabed over [-10; 10] km caused by the earthquake. The surface deformation is due to the hydro-acoustic waves propagating for early times (t = 49 s), then in the next two snapshots (t = 147 s and t = 343 s), we see the tsunami propagating in both directions. In the last snapshot (t = 833 s), the tsunami is away from the considered domain and it remains only acoustic waves.
<table><caption>TABLE 1. Parameter values for the first simulation: hydro-acoustic waves and tsunami generation by an earthquake.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain depth</td><td>1.5 km</td><td>s<sub>x</sub></td><td>150 m</td></tr><tr><td>T (convergence analysis)</td><td>50 s</td><td>r<sub>x</sub></td><td>15 km</td></tr><tr><td>Domain length (convergence analysis)</td><td>151.5 km</td><td>t<sub>0</sub></td><td>1 s</td></tr><tr><td>T</td><td>1000 s</td><td>s<sub>t</sub></td><td>0.05 s</td></tr><tr><td>Domain length</td><td>1576.5 km</td><td>r<sub>t</sub></td><td>1 s</td></tr><tr><td>g</td><td>9.81 m<sup>-1</sup></td><td>c<sub>0</sub></td><td>1500 m<sup>-1</sup></td></tr></tbody></table>

<table><caption>TABLE 2. Comparison of velocity-field and potential-based formulations: Mesh size, element order values, timestep, final time and computational time for the velocity-field simulation (first line) and the potential-based simulation (second line).</caption><thead><tr><th>Nx</th><th>Nz</th><th>Px</th><th>Pz</th><th>NDoFs</th><th>Δt</th><th>T</th><th>Computational time</th></tr></thead><tbody><tr><td>1051</td><td>10</td><td>4</td><td>5</td><td>428 910</td><td>9.6 10<sup>-3</sup> s</td><td>1000 s</td><td>771 s</td></tr><tr><td>1051</td><td>10</td><td>3</td><td>5</td><td>413 094</td><td>9.6 10<sup>-3</sup> s</td><td>1000 s</td><td>1197 s</td></tr></tbody></table>

For our simulation, the source has the form $u_b(x,t) = f(x)g(t)$, where $f$ and $g$ are smoothed rectangular functions such that the width and velocity of $u_b$ correspond to the source used in [27],

$$f(x) = \frac{1}{1 + e^{-(x-r_x)/s_x}} - \frac{1}{1 + e^{-(x+r_x)/s_x}}, \quad g(t) = \frac{1}{1 + e^{-(t-t_0)/s_t}} - \frac{1}{1 + e^{-(t-t_0-r_t)/s_t}}. \quad (5.23)$$

see Figure 2. Note that $x = 0$ corresponds to the middle of the domain. Natural boundary conditions are used on the lateral boundaries, and the computational domain is wider to avoid non-physical reflections inside the domain of interest. The numerical values are given in Table 1.

### 5.2.1. Convergence analysis

For each formulation, several elements orders are tested. The aim is to show that the simulations converge and to select parameters that are precise enough with a number of degree of freedoms as low as possible. The tsunami is simulated for 50 s and the computational domain is 151.5 km long. To evaluate whether convergence is reached, we plot the time evolution of the free surface displacement at a point $x = 50$ km away from the source center. The convergence is estimated visually and is assumed reached when increasing the order does not change the curves. The finite element orders are denoted by $P_x$, $P_z$ and the number of mesh subdivision in the $x$ and in the $z$ coordinates are respectively denoted by $N_x$, $N_z$. Several combinations of $P_x$ and $P_z$ are tested while keeping $N_x$ and $N_z$ fixed. The total number of Degrees of Freedom (NDoFs) for each case is indicated in each surface displacement plot.

For the velocity-field scheme, the surface displacement is shown in Figure 3. The convergence is reached for $P_x = 4$, $P_z = 5$. For the potential-based scheme, the surface displacement is shown in Figure 4a and we see that the orders with $P_x \ge 3$ converge. However, when taking $P_x = 3$, $P_z = 3$, we have observed increasing oscillations for longer times (not shown here). These oscillations come from the seabed displacement at $x = 0$ km, as shown in Figure 4b. The convergence is actually reached for $P_x = 3$, $P_z = 5$. Such oscillations are not observed for the velocity-field formulation.

### 5.2.2. Comparison and snapshots

We now run the simulation for a longer time and compare it with the results from [27]. The simulation obtained with the velocity-field formulation and the potential-based formulations are first compared with [27], then compared to each other. Finally we present some snapshots of the simulation.
<figure><img src="image_12.png" /><figcaption>(a) Vertical displacement at x =50 km obtained with the velocity-field (blue) and the potential-based (orange) formulations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Absolute difference (blue) and averaged absolute difference (orange) from Figure (a).</figcaption></figure>

Figure 6: Comparison of the velocity-field and the potential-based formulations: (a) superposition of the displacements, (b) absolute error.
<figure><img src="image_6.png" /><figcaption>FIGURE 3. Convergence analysis: Surface vertical displacement at $x = 50$ km for the velocity formulation at different elements orders.</figcaption></figure>

The tsunami is simulated for 1000 s, and the computational domain is 1576.5 km long. Since the final time is larger than for the convergence analysis, the domain size is also larger to avoid non-physical reflections at the lateral boundary and NDoFs increases. The element orders, the DoFs and the computational time are given in Table 2. We note that the potential-based formulation requires more computational time than the velocity-field formulation. The additional time is due to the displacement calculation: the displacement is obtained by integrating the velocity in time, and to obtain a precise integration, the velocity is computed at each time step. Hence the system (5.21) must be solved at each time step, leading to an additional computational cost. The code could be accelerated by using an efficient time integration rule which would not require to compute the velocity at each time step, or by directly solving the potential-based problem associated to the displacement.

For comparison with the literature, we consider the surface vertical displacement at a point $x = 50$ km evolving with time. The vertical displacement computed with the velocity-field model is superimposed to the result from [27] in Figure 5. Both curves have a high-frequency component starting at time $t \sim 20$ s, namely the acoustic waves, and a low-frequency wave starting at $t \sim 300$ s, which is the tsunami. The tsunami arrival times are in good agreement between the two curves. Since we do not have the exact description of the boundary condition used in [27], the earthquake description may slightly differs between the two simulations. This could explain the phase difference between both displacements.

To compare the potential-based and the velocity-field formulations, we plot in Figure 6a the same vertical displacement obtained from both formulations. The orange curve is the displacement $d_\Phi$ computed from the potential formulation, and the blue curve is the displacement $d_U$ computed from the velocity formulation. The curves are almost superimposed, which illustrates the fact that both formulations are equivalent. Figure 6b shows the difference $|d_\Phi - d_U|$ between the potential-based and the velocity-field formulations. The difference has a magnitude of less than 5 percents, and it oscillates with a period of approximately 4 seconds. The difference averaged over a period is plotted in the same figure. We see that the averaged difference slowly increases, which could be caused by the time integration: small errors add up and become non-negligible for large times.

Finally, to give some new illustrations of the wave propagation, we show in Figure 7 snapshots of the Eulerian domain near the earthquake at several times. The Eulerian domain is obtained by deforming the domain at rest with the displacement vector field. The colors correspond to the values of the vertical displacement. On each snapshot, we see the permanent deformation of the seabed over [-10; 10] km caused by the earthquake. The surface deformation is due to the hydro-acoustic waves propagating for early times ($t = 49$ s), then in the
<figure><img src="image_14.png" /><figcaption>Figure 7: Snapshots of the water displacement obtained with the potential formulation in the first simulation. The acoustic waves and the tsunami generated by a submarine earthquake induce a water displacement. The snapshots are taken at the following times, from top to bottom: $t = 49.0$ s, $t = 147.1$ s, $t = 343.3$ s, $t = 833.8$ s.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Convergence analysis: Displacement of (a) the water free surface at $x = 50$ km, (b) the seabed at $x = 0$ km, obtained from the potential-based formulation for several element orders.</figcaption></figure>

next two snapshots ($t = 147$ s and $t = 343$ s), we see the tsunami propagating in both directions. In the last snapshot ($t = 833$ s), the tsunami is away from the considered domain and it remains only acoustic waves.

## 5.3. Second simulation: rotational component of the fluid velocity

In this part, we study the validity of the traditional assumption that the velocity flow is irrotational. This assumption is made for example in [14, 25, 27]. The generalized potential yields a decomposition of the fluid velocity into a main irrotational component $\nabla\varphi$ and a remainder $U_r$,

$$ \mathbf{U} = -\nabla\varphi + \mathbf{U}_r, \quad \mathbf{U}_r = U_r \mathbf{e}_z = N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z. \qquad (5.24) $$

We recall that $N = \sqrt{-g\rho_0'(z)/\rho_0(z) - g^2/c_0^2(z)}$ is the buoyancy frequency [11, 15]. In the case $N \equiv 0$, the decomposition (5.24) shows that the flow is irrotational and the potential-based formulation (2.1)-(2.4) simplifies to the classical model (0.14)-(0.15). Hence the potential-based formulation can be seen as a generalization of
<table><caption>Table 4: Parameter values used in the second set of simulations for various earthquake scenarios: with and without topography, and using two different profiles for $N$.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain width</td><td>15 km</td><td>a</td><td>150 m</td></tr><tr><td>Domain height</td><td>1.5 km</td><td>x<sub>0</sub></td><td>7.5 km</td></tr><tr><td>Px</td><td>8</td><td>s<sub>x</sub></td><td>4 10<sup>-5</sup> m<sup>-2</sup></td></tr><tr><td>Pz</td><td>8</td><td>t<sub>0</sub></td><td>2 s</td></tr><tr><td>Nx</td><td>100</td><td>s<sub>t</sub></td><td>4 s<sup>-2</sup></td></tr><tr><td>Nz</td><td>10</td><td>f<sub>x</sub></td><td>0.07 m<sup>-1</sup></td></tr><tr><td>PML thickness</td><td>1.5 km</td><td>k<sub>x</sub></td><td>0.03 m<sup>-1</sup></td></tr><tr><td>N<sub>const</sub></td><td>0.001 s<sup>-1</sup></td><td>r<sub>x</sub></td><td>1.5 km</td></tr><tr><td></td><td></td><td>b</td><td>300 m</td></tr></tbody></table>

## 5.3 Second simulation: rotational component of the fluid velocity

In this part, we study the validity of the traditional assumption that the velocity flow is irrotational. This assumption is made for example in [25, 27, 13]. The generalized potential yields a decomposition of the fluid velocity into a main irrotational component $\nabla\varphi$ and a remainder $\mathbf{U}_r$,

$$ \mathbf{U} = -\nabla\varphi + \mathbf{U}_r, \quad \mathbf{U}_r = U_r \mathbf{e}_z = N \left( \psi + \frac{N}{g} \varphi \right) \mathbf{e}_z. \qquad (5.28) $$

We recall that $N = \sqrt{-g\rho_0'(z)/\rho_0(z) - g^2/c_0^2(z)}$ is the buoyancy frequency [10, 15]. In the case $N \equiv 0$, the decomposition (5.28) shows that the flow is irrotational and the potential-based formulation (2.1)-(2.4) simplifies to the classical model (0.14)-(0.15). Hence the potential-based formulation can be seen as a generalization of the system (0.14)-(0.15) to the case of a non-irrotational velocity. In this section, we use the potential-based formulation to quantify the contribution of the remainder $\mathbf{U}_r$ when the velocity is not strictly irrotational.

**Remark 5.1.** The remainder $\mathbf{U}_r$ defined in Equation (5.28) is not exactly the rotational part of the velocity as it cannot be written as $\mathbf{U}_r = \nabla \times \psi_r$.

We consider a scenario where one portion of the seafloor goes up and another portion goes down. This scenario is obtained with a source of the form $u_b(x, t) = f(x)g(t)$ where $f(x)$ is the derivative of a Gaussian and $g(t)$ is a Ricker function,

$$ f(x) = -2a s_x(x - x_0)e^{-(x-x_0)^2 s_x}, \quad g(t) = (4s_t^2(t - t_0)^2 - 2s_t)e^{-(t-t_0)^2 s_t}, \qquad (5.29) $$

see Figure 8. The potential-based simulation is used, and to avoid non-physical reflections, Perfectly Matched Layers (PML) were implemented. Two cases are tested: (1) with an initially flat seabed and with a constant buoyancy frequency $N_{const}$, (2) with an initial topography $z_b(x)$ and a depth-dependant $N(z)$. The topography $z_b$ consists of several bumps and flattens in the PML domain,

$$ z_b(x) = b(1 + \sin(k_x x)) \left( \frac{1}{1 + e^{-f_x(x-r_x/2)}} - \frac{1}{1 + e^{-f_x(x+r_x/2)}} \right). \qquad (5.30) $$

The profile $N(z)$ is obtained for profiles $c_0(z)$, $\rho_0(z)$ typically found in the ocean and described in Appendix 7.1. The simulation parameters are described in Table 4. For each case, we show a snapshot at time $t = 8$ s of several quantities: the irrotational component magnitude $|\nabla\varphi|$, the remainder $U_r$ and the ratio $|U_r|/|\mathbf{U}|$. Since all quantities fluctuate in time, we also show the time averaging of $|U_r|/|\mathbf{U}|$ on each DoFs over the whole simulated time.

The case $N \equiv N_{const}$ is shown in Figure 9. From the time averaging (Figure 9d), we see that the remainder is at least four orders of magnitude smaller than the main irrotational component, but that the ratio is not homogeneous in space. Its extremal values are on the seabed, on the upper half of the domain and away from