<figure><img src="image_8.png" /><figcaption>FIGURE 5. Vertical displacement at $x = 50$ km, computed from the velocity-field (blue) superimposed to the result in [27].</figcaption></figure>

<table><caption>TABLE 3. 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 depth</td><td>1.5 km</td><td>a</td><td>150 m s<sup>-1</sup></td></tr><tr><td>Domain length</td><td>15 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>106 m</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>0.1 s</td></tr><tr><td>Nz</td><td>10</td><td>f<sub>x</sub></td><td>15 m</td></tr><tr><td>T</td><td>120 s</td><td>k<sub>x</sub></td><td>0.03 m<sup>-1</sup></td></tr><tr><td>PML thickness</td><td>1.5 km</td><td>r<sub>x</sub></td><td>750 m</td></tr><tr><td>PML damping coefficient</td><td>10</td><td>b</td><td>300 m</td></tr><tr><td>N<sub>const</sub></td><td>0.001 s<sup>-1</sup></td><td></td><td></td></tr></tbody></table>

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 $U_r$ when the velocity is not strictly irrotational.

**Remark 5.1.** The remainder $U_r$ defined in Equation (5.24) is not exactly the rotational part of the velocity as it cannot be written as $U_r = \nabla \times \psi_r$.

We consider a scenario where a 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) = a f(x) g(t)$ where $f(x)$ is the derivative of a Gaussian and $g(t)$ is a Ricker function (second derivative of a Gaussian),

$$f(x) = -\frac{1}{s_x^2}(x - x_0) \exp\left(-\frac{(x - x_0)^2}{2s_x^2}\right) \quad g(t) = -\frac{1}{s_t^4}(t - t_0)^2 - \frac{1}{s_t^2} \exp\left(-\frac{(t - t_0)^2}{2s_t^2}\right) \quad (5.25)$$

see Figure 8. The potential-based simulation is used, and to avoid non-physical reflections, Perfectly Matched Layers (PML) were implemented in the spirit of [9]. A constant damping function was used in the PML. Two
<figure><img src="image_15.png" /><figcaption>Figure 8: The function $u_b(x, t)$ describing the seabed displacement for the second simulation. The same function $u_b$ is used for every earthquake scenario in this section.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 9: Earthquake on a flat seabed and with a constant $N \equiv 0.001 \text{ s}^{-1}$. Snapshots at $t = 8$s of: (A) the irrotational component magnitude $|\nabla\varphi|$, (B) the remainder $U_r$ and (C) the relative magnitude between the remainder and the total velocity $|U_r|/|\mathbf{U}|$. The time averaging over 120 s of $|U_r|/|\mathbf{U}|$ is shown in (D).</figcaption></figure>
<figure><img src="image_9.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_10.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) super-position of the displacements, (b) absolute difference.

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-dependent $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^{-(x-r_x)/f_x}} - \frac{1}{1 + e^{-(x+r_x)/f_x}} \right). \quad (5.26)$$

The profile $N(z)$ is obtained for profiles $c_0(z)$, $\rho_0(z)$ typically found in the ocean and described in Appendix 7.2. The simulation parameters are described in Table 3. Note that the time step $\Delta t$ depends on $\rho_0$ and $N$: for the case (1) we have $\Delta t = 2.9 \times 10^{-3}$ s and for the case (2) we have $\Delta t = 1.9 \times 10^{-3}$ s. For each case, we show a snapshot at time $t = 8$ s of several quantities: the irrotational component magnitude $|\nabla\phi|$, the remainder $U_r$ and the ratio $|U_r|/|U|$. Since all quantities fluctuate in time, we also show the time averaging of $|U_r|/|U|$ on
<figure><img src="image_17.png" /><figcaption>Figure 10: Earthquake on a topography and with a depth-dependent $N(z)$. Snapshots at $t = 8$s of: (A) the irrotational component magnitude $|∇φ|$, (B) the remainder $U_r$ and (C) the relative magnitude between the remainder and the total velocity $|U_r|/|U|$. The time averaging over 120 s of $|U_r|/|U|$ is shown in (D). The profile $N(z)$ (black curve) is shown on (B).</figcaption></figure>
<figure><img src="image_11.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_18.png" /><figcaption>Figure 11: Earthquake on a flat seabed and with a constant $N \equiv 0.01 \text{ s}^{-1}$: snapshots at $t = 8\text{s}$ of the remainder $U_r$</figcaption></figure>

the source location.

The case $N(z)$ is shown in Figure 10. Note that in this case, the profiles $c_0(z)$ and $\rho_0(z)$ also vary with depth, see Appendix 7.1. For reference, the profile $N(z)$ is plotted in Figure 10b. As before, the rotational component is much smaller than the irrotational component. The more complex profile $N(z)$ and the topography have an impact on the distribution of $U_r$. Indeed, the extreme values of $U_r$ are concentrated around the location of maximal $N(z)$ and on the off-peak parts of the seabed. We also notice that for both cases (Figures 9 and 10) the rotational component vanishes near the surface. Hence the irrotationality assumption seems verified for flows near the surface even when the inner flow is not strictly irrotational.

For the simulations presented here, the source time function $g(t)$ is a Ricker function, so that the time average of the source is zero. Other simulations made with a source having a non-zero time average show that the remainder $U_r$ and the ratio $|U_r|/|U|$ have their maximal values near the source location (see Appendix 7.2). Moreover, both quantities increase with time so that for long times ($t > 100 \text{ s}$) and close to the source, the ratio $|U_r|/|U|$ is of the same magnitude order as the irrotational component $\nabla\varphi$, even for small values of $N$. Hence the approximation $U = \nabla\varphi$ is not uniform in space and time and can become invalid near the source for particular choices of source functions.

Finally, to evaluate how the values of the buoyancy frequency affect the remainder $U_r$, we run the same scenario with $N \equiv 10 N_{\text{const}}$. A snapshot at time $t = 8 \text{ s}$ of the remainder is shown in Figure 11. We see that $U_r$ is distributed very similarly to the case $N \equiv N_{\text{const}}$, but is 100 times larger (compare to Figure 9b). It seems that $U_r$ scales approximately as $N^2$. From the equation (5.28) defining $U_r$, this observation could indicate that the component $N^2/g\varphi$ is dominant over the component $N\psi$.

## 5.4 Third simulation: interference patterns for submarine landslides

In this third scenario, we use the model to produce new simulations of waves generated by an submarine landslide. Many studies have been conducted on the generation of hydro-acoustic waves by submarine earthquakes (see [17, 1, 7] and the references therein). However, wave generation from submarine landslides have received less attention. In [5], it is highlighted that the hydro-acoustic signals generated by submarine landslides have a characteristic interference pattern. This pattern could be then used to detect and characterize landslides.

We start by a brief explanation about the interference pattern, then describe a simulation that reproduce the interference in a simplified case.

### 5.4.1 The Lloyd mirror effect

When acoustic waves propagate in a bounded medium, the reflection on the boundary create interference. This effect is also called the Lloyd mirror effect [20]. We consider a harmonic point source emitting rays in all directions, and neglect the reflection on the bottom. Every point in the domain is connected to the source by two rays: the direct ray, and the ray reflected by the surface (see Figure 12). For points far away from the source, the pressure field is approximated by

$$ \rho(x, z, t) \sim -A \frac{2i}{x} \sin(kz_s \sin\theta) e^{ikx} e^{-i\omega t}, \quad \text{for } x \gg z_s, \qquad (5.31) $$
<figure><img src="image_12.png" /><figcaption>FIGURE 8. The function $u_b(x, t)$ describing the seabed displacement for the second simulation. The same function $u_b$ is used for every earthquake scenario in this section.</figcaption></figure>

each DoFs over the whole simulated time. $T = 120$ s.

The case $N \equiv N_{\text{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 the source location.

The case $N(z)$ is shown in Figure 10. Note that in this case, the profiles $c_0(z)$ and $\rho_0(z)$ also vary with depth; see Appendix 7.2. For reference, the profile $N(z)$ is plotted in Figure 10b. As before, the rotational component is much smaller than the irrotational component. The more complex profile $N(z)$ and the topography have an impact on the distribution of $U_r$. Indeed, the extreme values of $U_r$ are concentrated around the location of maximal $N(z)$ and on the off-peak parts of the seabed. We also notice that for both cases (Figures 9 and 10) the rotational component vanishes near the surface. Hence the irrotationality assumption seems verified for flows near the surface even when the inner flow is not strictly irrotational.

For the simulations presented here, the source time function $g(t)$ is a Ricker function, so that the time average of the source is zero. Other simulations made with a source having an almost non-zero time average show that the remainder $U_r$ and the ratio $|U_r|/|\mathbf{U}|$ have their maximal values near the source location (see Appendix 7.3). Moreover, both quantities increase with time so that for long times ($t > 100$ s) and close to the source, the ratio $|U_r|/|\mathbf{U}|$ is of the same magnitude order as the irrotational component $\nabla\varphi$, even for small values of $N$. Hence the approximation $\mathbf{U} = \nabla\varphi$ is not uniform in space and time and can become invalid near the source for particular choices of source functions.

Finally, to evaluate how the values of the buoyancy frequency affect the remainder $U_r$, we run the same scenario with $N \equiv 10 N_{\text{const}}$. A snapshot at time $t = 8$ s of the remainder is shown in Figure 11. We see that $U_r$ is distributed very similarly to the case $N \equiv N_{\text{const}}$, but is 100 times larger (compare to Figure 9b). It seems that $U_r$ scales approximately as $N^2$. From the equation (5.24) defining $U_r$, this observation could indicate that the component $N^2/g\varphi$ is dominant over the component $N\psi$.
<figure><img src="image_19.png" /><figcaption>Figure 12: The direct path (in red) and the reflected path (in blue) connecting the source S to a point P. The location of the reflection (point R) is deduced from the image source S'. The declination angle θ is also indicated.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 13: The domain for the third simulation: acoustic waves and tsunami generation by a submarine landslide. The source is in black and the four receivers are indicated with 'x'.</figcaption></figure>

where A and z_s are respectively the source magnitude and depth, k is the wavenumber, θ is the declination angle, ω is the angular frequency. The pressure minima |p_max| = 0 are reached for

$$ \sin \theta = (m - 1) \frac{\pi}{k z_s}, \quad \text{with } m \in \mathbb{N}^*. \qquad (5.32) $$

Thanks to the relation (5.32), each point in space is associated to a "frequency bandwidth", namely the distance between two frequencies corresponding to minimum pressure. For a point with declination angle θ, the associated bandwidth Δf is given by

$$ \Delta f = \frac{c}{2z_s \sin \theta}. \qquad (5.33) $$

When recording the pressure at a fixed point and for a moving source (e.g. a landslide), the bandwidth Δf should vary with time. In [5], spectrograms computed from several fixed hydrophone show interference patterns. The measured bandwidths evolve with time, and are consistent with the theoretical value predicted by Equation (5.33). These results suggest that the Lloyd mirror effect could be used to detect submarine landslides and to recover their velocity.

### 5.4.2 Simulation description and result

We reproduce here the interference pattern for a static emitter in the 2D case with a flat seabed. Since the source is static, the bandwidth should not change with time. To illustrate the dependency of the interference bandwidth with the emitter-receiver distance, we record the pressure at four different locations. The four receivers are denoted R_i with i ∈ {1, 2, 3, 4}. The domain, with the emitter and receivers, is shown in Figure 13. The receivers locations and the corresponding theoretical bandwidth Δf, computed with Equation (5.33), are indicated in Table 6. The buoyancy frequency N is constant in this scenario.
<figure><img src="image_13.png" /><figcaption>FIGURE 9. Earthquake on a flat seabed and with a constant $N \equiv 0.001 \text{ s}^{-1}$. Snapshots at $t = 8$ s of: (A) the irrotational component magnitude $|\nabla\varphi|$, (B) the remainder $U_r$ and (C) the relative magnitude between the remainder and the total velocity $|U_r|/|U|$. The time averaging over 120 s of $|U_r|/|U|$ is shown in (D).</figcaption></figure>

## 5.4. Third simulation: interference patterns for submarine landslides

In this third scenario, we use the potential-based model to produce new simulations of waves generated by a submarine landslide. Many studies have been conducted on the generation of hydro-acoustic waves by submarine earthquakes (see [1, 7, 17] and the references therein). However, wave generation from submarine landslides have received less attention. In [5], it is highlighted that the hydro-acoustic signals generated by submarine landslides have a characteristic interference pattern. This pattern could be then used to detect and characterize landslides.

We start by a brief explanation about the interference pattern, then describe a simulation that reproduces the interference in a simplified case.

### 5.4.1. *The Lloyd mirror effect*

When acoustic waves propagate in a bounded medium, the reflection on the boundary create interference. This effect is also called the Lloyd mirror effect [21]. We consider a harmonic point source emitting rays in all directions, and neglect the reflection on the bottom. Every point in the domain is connected to the source by two rays: the direct ray, and the ray reflected by the surface (see Figure 12). For points far away from the
<table><caption>Table 5: Parameter values for the third simulation: acoustic waves and tsunami generation by a submarine landslide.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain width</td><td>150 km</td><td>A</td><td>1 ms<sup>-1</sup></td></tr><tr><td>Domain height</td><td>1.5 km</td><td>s<sub>x</sub></td><td>0.07 m<sup>-1</sup></td></tr><tr><td>Px</td><td>6</td><td>x<sub>0</sub></td><td>75 km</td></tr><tr><td>Pz</td><td>6</td><td></td><td></td></tr><tr><td>Nx</td><td>600</td><td></td><td></td></tr><tr><td>Nz</td><td>10</td><td></td><td></td></tr></tbody></table>

<table><caption>Table 6: Interference of acoustic waves generated by a submarine landslide: coordinates of the source center and the receivers, and theoretical Δf obtained with Equation (5.33).</caption><thead><tr><th></th><th>x (km)</th><th>z (km)</th><th>Theoretical Δf (Hz)</th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>5</td><td>0.3</td><td>2</td></tr><tr><td>R<sub>2</sub></td><td>8</td><td>0.3</td><td>3</td></tr><tr><td>R<sub>3</sub></td><td>5</td><td>1.35</td><td>16</td></tr><tr><td>R<sub>4</sub></td><td>8</td><td>1.35</td><td>26</td></tr></tbody></table>

It remains to define the source. Based on in-field data [5], the source is assumed to be a movement of the seabed over 600 m generating a continuous range of frequencies. In [5] the frequencies range up to 400 Hz, but for computational reason, we restrain the frequencies to 20 Hz maximum. The Lloyd-mirror effect is still visible for those lower frequencies. The function $u_b(x, t) = f(x)g(t)$, where f is a Gaussian function with magnitude A,

$$f(x) = A \exp(-s_x^2|x - x_0|^2),$$

and $g(t)$ is a white noise with frequencies ranging from 0 to 20 Hz, see Figure 14. For this simulation, the computational domain is long enough to avoid non-physical reflections. The simulation parameters are described in Table 5.

The spectrograms of the recorded pressures for each receiver are shown in Figure 15. In the four spectrograms, we see that the frequency ranges from 0 to 20 Hz, which is consistent with the source spectrogram. Constructive and destructive interferences are also clearly visible for each receiver. For the receivers $R_1$ and $R_2$, both at $z = 0.3$ km, the measured bandwidth is approximately 2 Hz. The receivers $R_3$ and $R_4$ are both at $z = 1.35$ km and their measured bandwidth is around 5 Hz. For the receivers near the surface ($R_3$ and $R_4$), the measured bandwidth do not correspond to the theoretical bandwidths presented in Table 6. One reason for this difference could be that the theoretical value is computed with the assumption of negligible reflections on the seabed. This assumption is probably not valid for our model, as the seabed is assumed rigid. Even though the numerical values do not correspond, we note that the measured bandwidths share some qualitative properties with the theoretical values: (1) the bandwidth is more sensitive to the depth of the hydro-acoustic sensor than to its distance to the landslide, and (2) the bandwidth is larger for the receivers closer to the surface.

In addition to the acoustic waves forming complex patterns, a tsunami is also generated. We show in Figure 16 the free-surface vertical displacement at a point 1.5 km away from the source center and the spectrogram of the signal. The tsunami is characterized on the spectrogram by the presence of a low frequency (below 5 Hz) component. We note that the low frequency is more prominent in the tsunami spectrogram than in the pressure spectrograms (see Figure 15). As in the previous tsunami simulation (see Figure 5), the curve shows a combination of high-frequency components starting at early times ($t \sim 10$ s) and a lower-frequency component starting at later time ($t \sim 40$ s). The high-frequency components correspond to the acoustic waves and the low-frequency one is the tsunami passing through the point.
<figure><img src="image_14.png" /><figcaption>FIGURE 10. Earthquake on a topography and with a depth-dependent $N(z)$. Snapshots at $t = 8s$ of: (A) the irrotational component magnitude $|∇φ|$, (B) the remainder $U_r$ and (C) the relative magnitude between the remainder and the total velocity $|U_r|/|U|$. The time averaging over 120 s of $|U_r|/|U|$ is shown in (D). The profile $N(z)$ (black curve) is shown on (B).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 11. Earthquake on a flat seabed and with a constant $N \equiv 0.01 \text{ s}^{-1}$: snapshots at $t = 8s$ of the remainder $U_r$</figcaption></figure>

source, the pressure field is approximated by

$$p(x, z, t) \sim -A \frac{2i}{x} \sin(kz_s \sin \theta) e^{ikx} e^{-i\omega t}, \quad \text{for } x \gg z_s, \qquad (5.27)$$
<figure><img src="image_21.png" /><figcaption>Figure 14: The function *g* generating a uniform range of frequencies up to 20 Hz: time serie (top) and spectrogram of *g* (bottom).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 15: Interference of acoustic waves generated by an submarine landslide: each spectrogram corresponds to one receiver. The coordinates for each receiver are, in km: $R_1=(5,0.3)$; $R_2=(8,0.3)$; $R_3=(5,1.35)$ and $R_4=(8,1.35)$. The measured bandwidth is indicated in white.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 12. The direct path (in red) and the reflected path (in blue) connecting the source S to a point P. The location of the reflection (point R) is deduced from the image source S'. The declination angle θ is also indicated.</figcaption></figure>

where A and z_s are respectively the source magnitude and depth, k is the wavenumber, θ is the declination angle, ω is the angular frequency. The pressure minima |p_max| = 0 are reached for

$$ \sin \theta = (m - 1) \frac{\pi}{k z_s}, \quad \text{with } m \in \mathbb{N}^*. \qquad (5.28) $$

Thanks to the relation (5.28), each point in space is associated to a “frequency bandwidth”, namely the distance between two frequencies corresponding to minimum pressure. For a point with declination angle θ, the associated bandwidth Δf is given by

$$ \Delta f = \frac{c}{2z_s \sin \theta}. \qquad (5.29) $$

When recording the pressure at a fixed point and for a moving source (e.g. a landslide), the bandwidth Δf should vary with time. In [5], spectrograms computed from several fixed hydrophone show interference patterns. The measured bandwidths evolve with time, and are consistent with the theoretical value predicted by Equation (5.29). These results suggest that the Lloyd mirror effect could be used to detect submarine landslides and to recover their velocity.

### 5.4.2. Computing the pressure

The pressure p is related to the velocity U by ∂_t p = −ρ₀c₀²∇·U, see [11]. This equation is written in weak form

$$ \frac{d}{dt} \int_{\Omega} p \tilde{p} dx = - \int_{\Omega} \rho_0 c_0^2 \nabla \cdot \mathbf{U} \tilde{p} dx, \quad \forall \tilde{p} \in L^2(\Omega). \qquad (5.30) $$

The weak form is then discretized. Let P^n denote the vector of degree of freedom for the discrete pressure. For simplicity, we choose P^n ∈ V_h. The pressure is computed with

$$ P^n = P^{n-1} + \Delta t \mathbb{M}_p^{-1} \mathbb{K}_p U^n, $$

where M_p is the mass matrix and K_p is the discretization of the right-hand side operation of (5.30).

### 5.4.3. Simulation description and result

We reproduce here the interference pattern for a static emitter in the 2D case with a flat seabed. The domain is assumed to be infinite in the x-direction (in practice we have used a very large domain and Neumann boundary conditions on the vertical boundaries). Since the source is static, the bandwidth should not change with time. To illustrate the dependency of the interference bandwidth with the emitter-receiver distance, we
<figure><img src="image_23.png" /><figcaption>Figure 16: Tsunami generated by a submarine landslide: time series (top) and spectrogram (bottom) of the surface vertical displacement for a point x = 1.5 km away from the source center.</figcaption></figure>

# 6 Conclusion and future work

In this work, we have presented two different formulations for a model describing the propagation of acoustic-gravity waves in a stratified ocean. The novel potential-based formulation is easier to handle from a mathematical point of view since it requires no essential boundary condition. Moreover, it offers several advantages for the numerical approximations: there are only 2 unknowns even for 3D problems, and the formulations include Neumann-type boundary conditions, easier to implement than the boundary condition of the velocity-field formulation – in particular when the seabed is not flat.

We have first shown the well-posedness of both formulations. For the velocity-field formulation, the study was complicated by the presence of a non-homogeneous boundary condition of Dirichlet type. To show the existence of solution, an equivalent dissipative problem was introduced. The equivalence between the velocity-field and potential-based formulations was then proved by using energy estimates. A discretization of both formulations using the spectral element method was presented, and the schemes were validated on a two-dimensional numerical test case. The test case also illustrated the equivalence between both formulation. We then used the relation between the velocity and the generalized potential to study the error made by the hypothesis of an irrotational flow when the fluid is not barotropic. Our numerical experiments suggest that the assumption is well-justified, but that it is not valid uniformly in space and time. Particular types of source can lead to a difference between the irrotational and the general model that increases in time on some very located portions of the domain. Finally, we presented a simulation bringing to light the interference pattern caused by the surface reflection of hydro-acoustic waves created by submarine landslides. Generation of hydro-acoustic waves by landslides is a topic that has not been very much explored yet, even though the analysis of those hydro-acoustic waves may provide a unique tool to detect and characterize submarine landslides and the potential hazard related to generated tsunamis. Notably, numerical experiments are very scarce.

The potential-based formulation could offer a new approach for the study of the Galbrun equation. Indeed, as mentioned in the introduction, the velocity-field problem can be seen as the Galbrun equation with no mean flow. For the Galbrun equation with a mean flow, the choice of the correct functional framework for the analysis is still an open question. It could be interesting to investigate whether a potential-based formulation could be obtained for the Galbrun equation with a mean flow. Moreover, from a numerical point of view, the transient problem with a mean flow and the harmonic problem – even without mean flow – are known to present spurious modes when discretized in a simple way [3, 2]. One could study the discretization of the potential-based formulation in the harmonic regime and check whether spurious modes are present.

Our simulations of a schematic landslide source and of the associated Lloyd mirror effect open new avenues to simulate in the same framework tsunamis and hydro-acoustic waves generated by landslides in the linear
<figure><img src="image_17.png" /><figcaption>FIGURE 13. The domain for the third simulation: acoustic waves and tsunami generation by a submarine landslide. The source is in black and the four receivers are indicated with 'x'.</figcaption></figure>

<table><caption>TABLE 4. Parameter values for the third simulation: acoustic waves and tsunami generation by a submarine landslide.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Domain dept</td><td>1.5 km</td><td>A</td><td>1 m s<sup>-1</sup></td></tr><tr><td>Domain length</td><td>150 km</td><td>s<sub>x</sub></td><td>106 m</td></tr><tr><td>Px</td><td>6</td><td>x<sub>0</sub></td><td>75 km</td></tr><tr><td>Pz</td><td>6</td><td>&Delta;t</td><td>6.2 10<sup>-3</sup> s</td></tr><tr><td>Nx</td><td>600</td><td>T</td><td>50 s</td></tr><tr><td>Nz</td><td>10</td><td></td><td></td></tr></tbody></table>

record the pressure at four different locations.

The four receivers are denoted by $R_i$ with $i \in \{1, 2, 3, 4\}$. The domain, with the emitter and receivers, is shown in Figure 13. The receivers locations and the corresponding theoretical bandwidth $\Delta f$, computed with Equation (5.29), are indicated in Table 5. The buoyancy frequency $N$ is constant in this scenario. It remains to define the source. Based on in-field data [5], the source is assumed to be a movement of the seabed over 400 m generating a continuous range of frequencies. In [5] the frequencies range up to 400 Hz, but for computational reason, we restrain the frequencies to 20 Hz maximum. The Lloyd-mirror effect is still visible for those lower frequencies. The function $u_b$ is given by $u_b(x, t) = f(x)g(t)$, where $f$ is a Gaussian function with magnitude $A$,

$$f(x) = A \exp \left( - \frac{|x - x_0|^2}{2s_x^2} \right)$$

and $g(t)$ is a white noise with frequencies up to 20 Hz (see Figure 14). For this simulation, the computational domain is long enough to avoid non-physical reflections. The simulation parameters are described in Table 4.

The spectrograms of the recorded pressures for each receiver are shown in Figure 15. In the four spectrograms, we see that the frequency ranges from 0 to 20 Hz, which is consistent with the source spectrogram. Constructive and destructive interferences are also clearly visible for each receiver. For the receivers $R_1$ and $R_2$, both at $z = 0.3$ km, the measured bandwidth is approximately 2 Hz. The receivers $R_3$ and $R_4$ are both at $z = 1.35$ km and their measured bandwidth is around 5 Hz. For the receivers near the surface ($R_3$ and $R_4$), the measured bandwidths do not correspond to the theoretical bandwidths presented in Table 5. One reason for this difference could be that the theoretical value is computed with the assumption of negligible reflections on the seabed. This assumption is probably not valid for our model, as the seabed is assumed rigid. Even though the numerical
approximation. Such simulations could ultimately provide insight into the most appropriate hydro-acoustic
sensor configuration in the field, making it possible to detect these events.

References

[1] F. Auclair, L. Debreu, E. Duval, M. Hilt, P. Marchesiello, E. Blayo, F. Dumas, and Y. Morel. Theory and analysis of acoustic-gravity waves in a free-surface compressible and stratified ocean. *Ocean Modelling*, 168:101900, December 2021.

[2] Anne-Sophie Bonnet-Ben Dhia, Kamel Berriri, and Patrick Joly. Régularisation de l'équation de Galbrun pour l'aéroacoustique en régime transitoire. *Revue Africaine de Recherche en Informatique et Mathématiques Appliquées*, Volume 5, Special Issue TAM'05:1855, September 2006.

[3] Anne-Sophie Bonnet-Ben Dhia, Eve-Marie Duclairoir, Guillaume Legendre, and Jean-François Mercier. Time-harmonic acoustic propagation in the presence of a shear flow. *Journal of Computational and Applied Mathematics*, 204(2):428-439, July 2007.

[4] Franco Brezzi and Michel Fortin. *Mixed and hybrid finite element methods*. Number 15 in Springer series in computational mathematics. Springer, New York Berlin Heidelberg, 1991.

[5] J. Caplan-Auerbach, R. P. Dziak, D. R. Bohnenstiehl, W. W. Chadwick, and T.-K. Lau. Hydroacoustic investigation of submarine landslides at West Mata volcano, Lau Basin. *Geophysical Research Letters*, 41(16):5927-5934, August 2014.

[6] Jacqueline Caplan-Auerbach, Christopher G. Fox, and Frederick K. Duennebier. Hydroacoustic detection of submarine landslides on Kilauea Volcano. *Geophysical Research Letters*, 28(9):1811-1813, 2001.

[7] C. Cecioni, G. Bellotti, A. Romano, A. Abdolali, P. Sammarco, and L. Franco. Tsunami Early Warning System based on Real-time Measurements of Hydro-acoustic Waves. *Procedia Engineering*, 70:311-320, 2014.

[8] Gary C Cohen. *Higher-Order Numerical Methods for Transient Wave Equations*. Scientific Computation. Springer Berlin, Heidelberg, 2001.

[9] Robert Dautray and Jacques-Louis Lions. *Mathematical analysis and numerical methods for science and technology. 5: Evolution problems*, volume 5. Springer, Berlin Heidelberg, nachdr. edition, 2000.

[10] Juliette Dubois, Sébastien Imperiale, Anne Mangeney, François Bouchut, and Jacques Sainte-Marie. Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible Euler equation. *Journal of Fluid Mechanics*, 970, 2023.

[11] Lawrence C. Evans. *Partial Differential Equations*. Graduate Studies in Mathematics. American Mathematical Society (AMS), 2004.

[12] Maurice Ewing, Ivan Tolstoy, and Frank Press. Proposed use of the T phase in tsunami warning systems*. *Bulletin of the Seismological Society of America*, 40(1):53-58, January 1950.

[13] Erez Eyov, Assaf Klar, Usama Kadri, and Michael Stiassnie. Progressive waves in a compressible-ocean with an elastic bottom. *Wave Motion*, 50(5):929-939, July 2013.

[14] International Association for the Properties of Water and Steam. IAPWS SR7-09, Supplementary Release on a Computationally Efficient Thermodynamic Formulation for Liquid Water for Oceanographic Use, 2009.

[15] Adrian E. Gill. *Atmosphere-ocean dynamics*. International geophysics series. Academic Press, New York, 1982.
<table><caption>TABLE 5. Interference of acoustic waves generated by a submarine landslide: coordinates of the source center and the receivers, and theoretical Δf obtained with Equation (5.29).</caption><thead><tr><th></th><th>x (km)</th><th>z (km)</th><th>Theoretical Δf (Hz)</th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>5</td><td>0.3</td><td>2</td></tr><tr><td>R<sub>2</sub></td><td>8</td><td>0.3</td><td>3</td></tr><tr><td>R<sub>3</sub></td><td>5</td><td>1.35</td><td>16</td></tr><tr><td>R<sub>4</sub></td><td>8</td><td>1.35</td><td>26</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>FIGURE 14. The function *g* generating a uniform range of frequencies up to 20 Hz: time serie (top) and spectrogram of *g* (bottom).</figcaption></figure>

values do not correspond, we note that the measured bandwidths share some qualitative properties with the theoretical values: (1) the bandwidth is more sensitive to the depth of the hydro-acoustic sensor than to its distance to the landslide, and (2) the bandwidth is larger for the receivers closer to the surface.

## 6. CONCLUSION AND FUTURE WORK

In this work, we have presented two different formulations for a model describing the propagation of acoustic-gravity waves in a stratified ocean. The novel potential-based formulation is easier to handle from a mathematical point of view since it requires no essential boundary condition. Moreover, it offers several advantages for the numerical approximations: there are only two unknowns even for 3D problems, and the formulation includes Neumann-type boundary conditions, easier to implement than the boundary condition of the velocity-field formulation – in particular when the seabed is not flat.

We have first shown the well-posedness of both formulations. For the velocity-field formulation, the study was complicated by the presence of a non-homogeneous boundary condition of Dirichlet type. To show the existence of solution, an equivalent dissipative problem was introduced. The equivalence between the velocity-field and potential-based formulations was then proved by using energy estimates. A discretization of both formulations using the spectral element method was presented, and the schemes were validated on a two-dimensional numerical test case. The test case also illustrated the equivalence between both formulations.
[16] Vivette Girault and Pierre-Arnaud Raviart. *Finite element methods for Navier-Stokes equations: theory and algorithms*. Number 5 in Springer series in computational mathematics. Springer, Berlin Heidelberg, 1. ed. 1986. edition, 1986.

[17] Bernabe Gomez and Usama Kadri. Near real-time calculation of submarine fault properties using an inverse model of acoustic signals. *Applied Ocean Research*, 109:102557, 2021.

[18] Gerd Grubb. *Distributions and operators*. Number 252 in Graduate texts in mathematics. Springer, New York, 2009. OCLC: ocn258078542.

[19] Linus Hägg and Martin Berggren. On the well-posedness of Galbrun's equation. *Journal de Mathématiques Pures et Appliquées*, 150:112–133, 2021.

[20] Finn B. Jensen, William A. Kuperman, Michael B. Porter, and Henrik Schmidt. *Computational Ocean Acoustics*. Springer New York, New York, NY, 2011.

[21] Patrick Joly, C. Tsogka, G. Derveaux, and J. Rodriguez. *Effective Computational Methods for Wave Propagation, Part. 3 Numerical Methods for Elastic Wave Propagation*. Chapman and Hall/CRC, 1 edition, 2008.

[22] Dimitri Komatitsch and Jeroen Tromp. Introduction to the spectral element method for three-dimensional seismic wave propagation. *Geophysical Journal International*, 139(3):806–822, December 1999.

[23] Dimitri Komatitsch and Jean-Pierre Vilotte. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. *Bulletin of the Seismological Society of America*, 88(2):368–392, April 1998.

[24] J. L. Lions and E. Magenes. *Non-Homogeneous Boundary Value Problems and Applications*. Springer Berlin Heidelberg, Berlin, Heidelberg, 1972.

[25] M. S. Longuet-Higgins. A Theory of the Origin of Microseisms. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 243(857):1–35, 1950.

[26] Marcus Maeder, Gwénaël Gabard, and Steffen Marburg. 90 Years of Galbrun's Equation: An Unusual Formulation for Aeroacoustics and Hydroacoustics in Terms of the Lagrangian Displacement. *Journal of Theoretical and Computational Acoustics*, 28(04):2050017, December 2020.

[27] P. Sammarco, C. Cecioni, G. Bellotti, and A. Abdolali. Depth-integrated equation for large-scale modelling of low-frequency hydroacoustic waves. *Journal of Fluid Mechanics*, 722:R6, 2013.

# 7 Appendix

## 7.1 A realistic profile for the buoyancy frequency N

In [10], we have shown how profiles for $\rho(z)$, $c(z)$ can be computed from a given temperature profile. For completeness we recall here the equations. The compressible Euler equations under the assumption of an initial state at equilibrium yield the following system of equations for the pressure $p_0$, temperature $T_0$ and density $\rho_0$,

$$ \nabla p_0 = -\rho_0 g \mathbf{e}_z, \quad \rho_0 = f_\rho(p_0, T_0), \quad p_0(H) = p^a. \qquad (7.1) $$

For the sound speed we use the expression $c_0 = c(p_0, T_0)$ given in [14]. Hence, if $T_0(z)$ and an equation of state $f_\rho$ is given, the system

$$ \frac{dp_0}{dz} = -g f_\rho(p_0, T_0), \quad z \in (0, H), \qquad (7.2) $$

$$ p_0 = p^a, \quad z = H. \qquad (7.3) $$
<figure><img src="image_19.png" /><figcaption>FIGURE 15. Interference of acoustic waves generated by a submarine landslide: each spectrogram corresponds to one receiver. The coordinates for each receiver are, in km: $R_1=(5,0.3)$; $R_2=(8,0.3)$; $R_3=(5,1.35)$ and $R_4=(8,1.35)$. The measured bandwidth is indicated in white.</figcaption></figure>

We then used the relation between the velocity and the generalized potential to study the error made by the hypothesis of an irrotational flow when the fluid is not barotropic. Our numerical experiments suggest that the assumption is well-justified, but that it is not valid uniformly in space and time. Particular types of source can lead to a difference between the irrotational and the general model that increases in time on some very located portions of the domain. Finally, we presented a simulation bringing to light the interference pattern caused by the surface reflection of hydro-acoustic waves created by submarine landslides. Generation of hydro-acoustic waves by landslides is a topic that has not been very much explored yet, even though the analysis of those hydro-acoustic waves may provide a unique tool to detect and characterize submarine landslides and the potential hazard related to generated tsunamis. Notably, numerical experiments are very scarce.

The potential-based formulation could offer a new approach for the study of the Galbrun equation. Indeed, as mentioned in the introduction, the velocity-field problem can be seen as the Galbrun equation with no mean flow. For the Galbrun equation with a mean flow, the choice of the correct functional framework for the analysis is still an open question. It could be interesting to investigate whether a potential-based formulation could be obtained for the Galbrun equation with a mean flow. Moreover, from a numerical point of view, the transient problem with a mean flow and the harmonic problem – even without mean flow – are known to present spurious modes when discretized in a simple way [2,3]. One could study the discretization of the potential-based formulation in the harmonic regime and check whether spurious modes are present.

Our simulations of a schematic landslide source and of the associated Lloyd mirror effect open new avenues to simulate in the same framework tsunamis and hydro-acoustic waves generated by landslides in the linear approximation. Such simulations could ultimately provide insight into the most appropriate hydro-acoustic sensor configuration in the field, making it possible to detect these events.

## REFERENCES

[1] F. Auclair, L. Debreu, E. Duval, M. Hilt, P. Marchesiello, E. Blayo, F. Dumas, and Y. Morel. Theory and analysis of acoustic-gravity waves in a free-surface compressible and stratified ocean. 168:101900.

[2] Anne-Sophie Bonnet-Ben Dhia, Kamel Berriri, and Patrick Joly. Régularisation de l'équation de galbrun pour l'aéroacoustique en régime transitoire. Volume 5, Special Issue TAM'05:1855.

[3] Anne-Sophie Bonnet-Ben Dhia, Eve-Marie Duclairoir, Guillaume Legendre, and Jean-François Mercier. Time-harmonic acoustic propagation in the presence of a shear flow. 204(2):428-439.

[4] Franco Brezzi and Michel Fortin. *Mixed and hybrid finite element methods*. Number 15 in Springer series in computational mathematics. Springer.

[5] J. Caplan-Auerbach, R. P. Dziak, D. R. Bohnenstiehl, W. W. Chadwick, and T.-K. Lau. Hydroacoustic investigation of submarine landslides at west mata volcano, lau basin. 41(16):5927-5934.
<figure><img src="image_24.png" /><figcaption>Figure 17: Profiles used in Section 5.3: (a) Temperature $T(z)$, (b) density $\rho_0(z)$ and (c) sound speed $c_0(z)$ used for the computation of (d) the buoyancy frequency $N(z)$.</figcaption></figure>

is solved to obtain $\rho_0$. Then $\rho_0$ and $c_0$ are computed from (7.1) and [14]. The equation of state $f_\rho$ is given in [14]. The differential equation for the pressure (7.3) is numerically solved for the temperature profile shown in Figure 17a, and we obtain the profiles $\rho_0(z)$ and $c_0(z)$ shown in Figure 17b and 17c. The buoyancy frequency $N$ is then computed from (0.3) and we obtain the profile shown in Figure 17d.

## 7.2 The irrotational flow assumption

We study the same scenario as in Section 5.3 but with a different function $g(t)$ representing the time dependency of the source $u_b$. Here $g(t)$ is a smoothed rectangle,

$$g(t) = \frac{1}{1 + e^{-s_t(t-t_0)}} - \frac{1}{1 + e^{-s_t(t-t_0-r_t)}}$$

with $s_t = 20 \text{ s}^{-1}$, $t_0 = 1 \text{ s}$ and $r_t = 1 \text{ s}$. The difference with Section 5.3 is that the time average of the source $u_b$ is not zero. The function $u_b(x, t)$ is shown in Figure 18. The simulation uses the topography described by Equation (5.30) and the profile $N(z)$ described in Section 7.1. Snapshots of the remainder $|U_r|$ (see Figure 19) show that even though being small, the values near the source increase with time.
[6] Jacqueline Caplan-Auerbach, Christopher G. Fox, and Frederick K. Duennebier. Hydroacoustic detection of submarine land-slides on kilaua'e volcano. 28(9):1811-1813.

[7] C. Cecioni, G. Bellotti, A. Romano, A. Abdolali, P. Sammarco, and L. Franco. Tsunami early warning system based on real-time measurements of hydro-acoustic waves. 70:311-320.

[8] Gary C Cohen. *Higher-Order Numerical Methods for Transient Wave Equations*. Scientific Computation. Springer Berlin, Heidelberg.

[9] Gary C Cohen and Sébastien Imperiale. Perfectly matched layer with mixed spectral elements for the propagation of linearized water waves. *Communications in Computational Physics*, 11(2):285-302, 2012.

[10] Robert Dautray and Jacques-Louis Lions. *Mathematical analysis and numerical methods for science and technology. 5: Evolution problems*, volume 5. Springer, nachdr. edition.

[11] Juliette Dubois, Sébastien Imperiale, Anne Mangeney, François Bouchut, and Jacques Sainte-Marie. Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible euler equation. 970.

[12] Lawrence C. Evans. *Partial Differential Equations*. Graduate Studies in Mathematics. American Mathematical Society (AMS).

[13] Maurice Ewing, Ivan Tolstoy, and Frank Press. Proposed use of the t phase in tsunami warning systems*. 40(1):53-58.

[14] Erez Eyov, Assaf Klar, Usama Kadri, and Michael Stiassnie. Progressive waves in a compressible-ocean with an elastic bottom. 50(5):929-939.

[15] Adrian E. Gill. *Atmosphere-ocean dynamics*. International geophysics series. Academic Press.

[16] Vivette Girault and Pierre-Arnaud Raviart. *Finite element methods for Navier-Stokes equations: theory and algorithms*. Number 5 in Springer series in computational mathematics. Springer, 1. ed. 1986. edition.

[17] Bernabe Gomez and Usama Kadri. Near real-time calculation of submarine fault properties using an inverse model of acoustic signals. 109:102557.

[18] Gerd Grubb. *Distributions and operators*. Number 252 in Graduate texts in mathematics. Springer. OCLC: ocn258078542.

[19] Martin Halla and Thorsten Hohage. On the well-posedness of the damped time-harmonic galbrun equation and the equations of stellar oscillations. 53(4):4068-4095.

[20] Linus Hägg and Martin Berggren. On the well-posedness of galbrun's equation. 150:112-133.

[21] Finn B. Jensen, William A. Kuperman, Michael B. Porter, and Henrik Schmidt. *Computational Ocean Acoustics*. Springer New York.

[22] Patrick Joly, C. Tsogka, G. Derveaux, and J. Rodriguez. *Effective Computational Methods for Wave Propagation, Part. 3 Numerical Methods for Elastic Wave Propagation*. Chapman and Hall/CRC, 1 edition.

[23] Dimitri Komatitsch and Jeroen Tromp. Introduction to the spectral element method for three-dimensional seismic wave propagation. 139(3):806-822.

[24] J. L. Lions and E. Magenes. *Non-Homogeneous Boundary Value Problems and Applications*. Springer Berlin Heidelberg.

[25] M. S. Longuet-Higgins. A theory of the origin of microseisms. 243(857):1-35.

[26] Marcus Maeder, Gwénaël Gabard, and Steffen Marburg. 90 years of galbrun's equation: An unusual formulation for aeroacoustics and hydroacoustics in terms of the lagrangian displacement. 28(4):2050017.

[27] P. Sammarco, C. Cecioni, G. Bellotti, and A. Abdolali. Depth-integrated equation for large-scale modelling of low-frequency hydroacoustic waves. 722:R6.

[28] UNESCO. Tenth report of the joint panel on oceanographic tables and standards.

# 7. APPENDIX

## 7.1. Further details on operators

### 7.1.1. Properties of $G$ and $\tilde{G}$

We recall that $G: \mathcal{D}(G) \subset \mathcal{H} \to \mathcal{G}$ is defined by

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

where $\mathcal{H} = L^2(\Omega)^d$ and $\mathcal{G} = L^2(\Omega) \times L^2(\Omega) \times L^2(\Gamma_s)$, moreover

$$ \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 (7.1) $$
<figure><img src="image_25.png" /><figcaption>Figure 18: The source $u_b(x, t)$ over time and space.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 19: Snapshots of $U_r$ at various times: (a) $t = 8$ s, (b) $t = 60$ s and (c) $t = 120$ s. Note how the color scale changes at each snapshot.</figcaption></figure>
Smooth functions compactly supported in $\Omega$ belong to $\mathcal{D}(G)$ and are dense in $L^2(\Omega)^d$, this implies that $G$ is densely defined. We recall that $\tilde{G}$ is an extension of $G$, its domain of definition is

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

hence it is densely defined since $\mathcal{D}(G) \subset \mathcal{D}(\tilde{G})$. We now show that $\tilde{G}$ is a closed operator. Let $\mathbf{U}_n$ be a sequence of function in $\mathcal{D}(\tilde{G})$ such that $\mathbf{U}_n$ converges towards $\mathbf{U} \in \mathcal{H}$ and $\tilde{G}\mathbf{U}_n$ converges towards $\mathbf{\Psi} \in \mathcal{G}$. Our aim is to show that $\mathbf{U} \in \mathcal{D}(\tilde{G})$ and $\mathbf{\Psi} = \tilde{G}\mathbf{U}$.

Looking at the first component of (1.5) or (7.1), we have the existence of $\psi \in L^2(\Omega)$ such that $\nabla \cdot \mathbf{U}_n$ converges to $\psi$ in $L^2(\Omega)$. Since $\mathbf{U}_n \in H(\mathrm{div}, \Omega)$, one can show – using test functions in $\mathcal{D}(\Omega)$ – that $\mathbf{U} \in H(\mathrm{div}, \Omega)$ and $\psi = \nabla \cdot \mathbf{U}$ (i.e. the divergence operator from $H(\mathrm{div}, \Omega)$ to $L^2(\Omega)$ is a closed operator). Now, since the normal trace operator $\gamma_{1,s}$ is continuous from $H(\mathrm{div}, \Omega)$ to $H^{-1/2}(\Gamma_s)$, we have $\gamma_{1,s}(\mathbf{U}) \in H^{-1/2}(\Gamma_s)$ and

$$ \mathbf{U}_n \xrightarrow{L^2(\Omega)^d} \mathbf{U}, \quad \nabla \cdot \mathbf{U}_n \xrightarrow{L^2(\Omega)} \nabla \cdot \mathbf{U} \Rightarrow \gamma_{1,s}(\mathbf{U}_n) \xrightarrow{H^{-1/2}(\Gamma_s)} \gamma_{1,s}(\mathbf{U}). $$

Moreover, by assumption $\gamma_{1,s}(\mathbf{U}_n)$ converges towards an element in $L^2(\Gamma_s)$, therefore, by identification $\gamma_{1,s}(\mathbf{U})$ belongs to $L^2(\Gamma_s)$. We have shown that $\mathbf{U} \in \mathcal{D}(\tilde{G})$ and $\tilde{G}\mathbf{U}_n$ converges towards $\tilde{G}\mathbf{U}$. The fact that $G$ is closed can be proven similarly.

### 7.1.2. The range of G and $G^*$

We recall and prove the theorem 3.1.

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

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

$$ \|G\mathbf{U}\|_{\tilde{G}}^2 \geq C \|\mathbf{U}\|_{\mathcal{H}}^2 \quad (7.2) $$

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

$$ \mathrm{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 \mathrm{Ker}(G)$, the function $U_x$ is constant in the $x$ direction. Since $U_x \in L^2(\Omega)$ and the domain is unbounded in the $x$ direction, the function $U_x$ is equal to zero. Hence $\mathcal{D}(G^*) \cap \mathrm{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(\mathbf{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)}. $$

We let $\mathbf{U}_n \in \mathcal{D}(G)$ denote the following divergence-free function:

$$ \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)^\top$ and

$$ \|G\mathbf{U}_n\|_{\tilde{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\|_{\mathcal{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. \quad (7.3) $$

For $n$ going to infinity we have $\|G\mathbf{U}_n\|_{\tilde{G}}^2$ goes to zero and $\|\mathbf{U}_n\|_{\mathcal{H}}^2$ converges to $\|\sqrt{\rho_0} \partial_z u\|_{L^2(\Omega)}^2$, hence (7.2) cannot hold. $\square$
<figure><img src="image_20.png" /><figcaption>FIGURE 16. Profiles used in Section 5.3: (a)Temperature T(z), (b) density ρ₀(z) and (c) sound speed c₀(z) used for the computation of (d) the buoyancy frequency N(z).</figcaption></figure>

## 7.2. A realistic profile for the buoyancy frequency N

In [11], we have shown how profiles for $ρ(z)$, $c(z)$ can be computed from a given temperature profile. For completeness we recall here the equations. The compressible Euler equations under the assumption of an initial state at equilibrium yield the following system of equations for the pressure $p_0$, temperature $T_0$ and density $ρ_0$,

$$ \nabla p_0 = -\rho_0 g \mathbf{e}_z, \quad \rho_0 = f_\rho(p_0, T_0), \quad p_0(H) = p^a. \qquad (7.4) $$

For the sound speed we use the expression $c_0 = c(p_0, T_0)$ given in [28]. Hence, if $T_0(z)$ and an equation of state $f_\rho$ is given, the system

$$ \frac{dp_0}{dz} = -g f_\rho(p_0, T_0), \quad z \in (0, H), \qquad (7.5) $$

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

is solved to obtain $p_0$, and then $ρ_0$ and $c_0$ are computed from (7.4) and [28]. The equation of state $f_\rho$ is given in [28]. The differential equation for the pressure (7.6) is numerically solved for the temperature profile shown in Figure 16a, and we obtain the profiles $ρ_0(z)$ and $c_0(z)$ shown in Figure 16b and 16c. The buoyancy frequency N is then computed from (0.3) and we obtain the profile shown in Figure 16d.
<figure><img src="image_21.png" /><figcaption>FIGURE 17. The source $u_b(x, t)$ over time and space.</figcaption></figure>

## 7.3. The irrotational flow assumption

We study the same scenario as in Section 5.3 but with a different function $g(t)$ representing the time dependency of the source $u_b$. Here $g(t)$ is a smoothed rectangle,

$$g(t) = \frac{1}{1 + e^{-(t-t_0)/s_t}} - \frac{1}{1 + e^{-(t-t_0-r_t)/s_t}},$$

with $s_t = 0.05$ s, $t_0 = 1$ s and $r_t = 1$ s. The difference with Section 5.3 is that the time average of the source $u_b$ is not zero. The function $u_b(x, t)$ is shown in Figure 17. The simulation uses the topography described by Equation (5.26) and the profile $N(z)$ described in Section 7.2. Snapshots of the remainder $|U_r|$ (see Figure 18) show that even though being small, the values near the source increase with time.
<figure><img src="image_22.png" /><figcaption>FIGURE 18. Snapshots of $U_r$ at various times: (a) $t = 8$ s, (b) $t = 60$ s and (c) $t = 120$ s.
Note how the color scale changes at each snapshot.</figcaption></figure>