# Three-dimensional modelling of internal tide generation over isolated seamounts in a rotating ocean

Cécile Le Dizes, Nicolas Grisouard, Olivier Thual, Matthieu Mercier

► To cite this version:

Cécile Le Dizes, Nicolas Grisouard, Olivier Thual, Matthieu Mercier. Three-dimensional modelling of internal tide generation over isolated seamounts in a rotating ocean. Journal of Fluid Mechanics, 2025, 1022, pp.A5. 10.1017/jfm.2025.10647. hal-05440264

HAL Id: hal-05440264

https://hal.science/hal-05440264v1

Submitted on 17 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Three-dimensional modelling of internal tide generation over isolated seamounts in a rotating ocean

Cécile Le Dizes¹ ID, Nicolas Grisouard² ID, Olivier Thual¹ ID and Matthieu J. Mercier¹ ID

¹Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, INPT, Toulouse, France

²Department of Physics, University of Toronto, Toronto, ON, Canada

**Corresponding authors:** Nicolas Grisouard, nicolas.grisouard@utoronto.ca; Cécile Le Dizes, cecile.ledizes@imft.fr

(Received 28 March 2025; revised 20 August 2025; accepted 24 August 2025)

This study introduces a boundary element method to solve the three-dimensional problem of internal tide generation over arbitrary isolated seamounts in a uniformly stratified finite-depth fluid with background rotation, without assumptions on the size or slope of the topography. Focusing on linearly propagating waves with small tidal excursions, the approach employs a vertical mode decomposition to describe the wavefield and the wave energy flux. We apply the model to the generation of internal tides by a unidirectional barotropic tide interacting with an axisymmetric Gaussian seamount. We study the conversion rate and flow field for various topographic configurations. We qualitatively recover some of the two-dimensional results of Papoutsellis et al. (2023 *J. Fluid Mech.* **964**, A20), and find topographies with weak conversion rates, as discussed by Maas (2011 *J. Fluid Mech.* **684**, 5–24). Furthermore, our results reveal the previously underestimated influence of the Coriolis frequency on the wavefield and on the spatial distribution of radiated energy flux. Due to Coriolis effects, the energy fluxes are shifted slightly counter-clockwise in the northern hemisphere. We explain in detail how this shift increases with the magnitude of the Coriolis frequency and the topographic features and why such effects are absent in models based on the weak topography assumption.

**Key words:** internal waves, waves in rotating fluids, boundary integral methods

## 1. Introduction

Internal tides are inertia-gravity waves that propagate within the stratified ocean, generated by the interaction of barotropic tidal currents with underwater topography (Garrett &

© The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Kunze 2007; Musgrave et al. 2022). These waves are associated with significant vertical displacements and represent a major source of baroclinic energy that contributes to the mixing of the stratified fluid when breaking and dissipating (Munk & Wunsch 1998; Egbert & Ray 2000, 2001; Whalen et al. 2020). This mixing, in turn, affects global ocean circulation (Wunsch & Ferrari 2004) and ecosystems through nutrient transport (Sandstrom & Elliott 1984). However, accurately modelling internal tide dynamics is challenging due to the complex interactions across a wide range of scales. Global models typically lack the spatial resolution required to solve these processes explicitly and must rely on parametric approaches to represent the mixing induced by internal-wave breaking (MacKinnon et al. 2017). These choices are critical for climate scenarios, as it can significantly impact the global-scale dynamics (Melet et al. 2013, 2016). Parametrisations traditionally rely on local energy conversion rates, which correspond to the energy transmitted from the barotropic tide to the internal waves (St. Laurent, Simmons & Jayne 2002), to represent the contribution of internal tides to mixing near the generation sites. More recent models describe also the horizontal propagation of internal-wave energy and various energy sinks (Eden & Olbers 2014; de Lavergne et al. 2019), allowing for the estimation of the far-field energy dissipation. In this context, the geometry of energy dissipation is strongly related to local generation maps, highlighting the need for approaches that accurately capture the variability of internal tide generation.

One widely used framework to estimate the energy conversion rate is referred to as the weak topography assumption (WTA). This approach assumes that the topographic height is much smaller than the total fluid depth and that internal waves are subcritical, i.e. propagate at an angle to the horizontal larger than the topographic slope. These assumptions allow for the linearisation of the non-penetration condition at the seafloor, making the problem analytically tractable. Bell (1975a,b) developed the first WTA model for a semi-infinite ocean, which was later extended to a bounded ocean by Llewellyn Smith & Young (2002), using the decomposition into vertical modes. Their work provided estimates of the conversion rate for any three-dimensional (3-D) seamount, which were used for global calculations of tidal energy conversion, such as those by Nycander (2005), Green & Nycander (2013) and Falahat et al. (2014).

To extend these results to more realistic cases, several other approaches have been developed to address the generation of internal tides at supercritical topographies, namely, when the maximum topographic slope is larger than the angle of internal-wave beams to the horizontal. These include analytical methods based on the assumption of infinitely steep topography like the knifeedge or the top-hat barrier (St. Laurent et al. 2003; Llewellyn Smith & Young 2003; Baines 2007), the body-forcing approach (Baines 1973) later revised and corrected by Garrett & Gerkema (2007), asymptotic expansions (Balmforth, Ierley & Young 2002) or coupled-mode approaches (Griffiths & Grimshaw 2007; Kelly 2016; Papoutsellis, Mercier & Grisouard 2023); see Garrett & Kunze (2007) for additional references. Among these, Pétrélis et al. (2006), Echeverri & Peacock (2010) and Mathur, Carter & Peacock (2016) adapted the boundary element method (BEM), commonly used in elasticity, acoustics and electromagnetism, for internal tide generation over topography of various steepness. Their method offers solutions for internal tide generation over both subcritical and supercritical topographies, but it is restricted to 2-D ridges.

More recently, efforts have been made to account not only for the conversion rate, but also for the directional dependence of the baroclinic energy radiated by 3-D topographies. The work of Llewellyn Smith & Young (2002) first showed that the conversion rate differs for symmetric and asymmetric 3-D seamounts depending on the orientation of the barotropic tide, but it did not provide information on the direction of the energy flux.
This missing information is now discussed in Pollmann et al. (2019), using the same WTA approach to describe the energy flux direction, emphasising the dependence on the shape of the topography and the orientation of the barotropic tidal ellipse. This modelling approach has been applied to the global ocean by Pollmann & Nycander (2023), allowing us to decompose, in terms of directions, the barotropic-to-baroclinic energy conversion. However, since it is restricted to the WTA, it leads to invalid results at steep continental slopes and large seamounts, especially for supercritical cases (Geoffroy, Pollmann & Nycander 2024). The 3-D modelling at supercritical topographies has first been addressed using step-like axisymmetric features, the pillbox-shaped seamounts (Baines 2007). The internal waves radiated resemble the beams radiated at 2-D step topographies, but there is variation in amplitude and phase with azimuthal angle. Although it is not the focus of this paper, the author indicates that the direction of maximum flux is not always aligned with that of the barotropic tide, due to the influence of background rotation. The boundary integral representation has also been used for 3-D supercritical bodies by Voisin, Ermanyuk & Flor (2011), Voisin (2021) and Voisin (2024). In this last paper, the author derives the internal waves generated by translational oscillations of an ellipsoidal seamount in a semi-infinite non-rotating and uniformly stratified ocean. In particular, he finds that the conversion rate for the spheroid of length *a* and height *c* can be orders of magnitude lower than its WTA equivalent for large values of the ratio $γ = µc/a$, where $µ^{-1}$ is the slope of the wave beams. He also studies the effects of the anisotropy of the seamount, which are found to be similar to the ones observed by Llewellyn Smith & Young (2002), as well as the influence of the orientation of the tidal ellipse with the major axis of the topography. The studies of Baines (2007) and Voisin (2024) are the only analytical studies that tackle internal-wave generation for 3-D supercritical topographies, but both are limited to the infinitely steep slope limit.

The validation of these analytical approaches with observations is often difficult due to the complex nature of oceanic bathymetry and tidal forcing. Nevertheless, the use of analytical inputs remain crucial to investigate the internal tide lifecycle in specific regions (Vic et al. 2018; Waterhouse et al. 2018; Savage, Waterhouse & Kelly 2020). Better comparisons can be obtained with direct numerical simulations of idealised scenarios, which indeed confirmed that non-axisymmetric topographies convert more energy into internal tides (Holloway & Merrifield 1999; Munroe & Lamb 2005), as discussed earlier. Another important result is that axisymmetric seamounts can generate waves with asymmetric patterns of energy distribution (Munroe & Lamb 2005), induced by the background rotation (Coriolis frequency). However, no clear explanation has been reported. Some specific studies with no background rotation have also focused on the radiated 3-D internal waves, focusing on weakly nonlinear aspects (King, Zhang & Swinney 2010) or on the interactions of multiple isolated seamounts (Zhang et al. 2017).

In this study, we address the limitations of past analytical approaches by developing a boundary element method to solve the 3-D problem of internal tide generation over an isolated seamount in an ocean with finite depth. In particular, our approach allows us to account for arbitrary seamount geometries and slopes, by taking into account the full bottom boundary condition, similarly to the approach developed in 2-D studies (Echeverri & Peacock 2010). The use of boundary integral equations enables an accurate representation of the flow field at a reasonable computational cost. We set the problem in § 2 before providing the analytical details on how to solve the internal tide generation problem in three dimensions, based on the Green's function approach, in § 3. The properties of the internal tide generated by a linear barotropic forcing over a Gaussian seamount are presented in detail in § 4, focusing on the influence of the topographic features on the conversion rate as well as the orientation of the radiated tide (§ 5).
C. Le Dizes, N. Grisouard, O. Thual and M.J. Mercier

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic of the problem with representations of the first three vertical modes $a_1, a_2, a_3$ for constant stratification.</figcaption></figure>

We explain in § 6 the origin of the peculiar azimuthal dependence of the energy flux radiated, influenced by the background rotation and criticality, before concluding in § 7.

## 2. Governing equations

### 2.1. Wave equation

We investigate the internal tides generated by a spatially uniform barotropic tide, interacting with an isolated seamount of horizontal and vertical scales $L$ and $\Lambda$ (see figure 1). For the sake of simplicity, we will consider a rectilinear tide $\tilde{U}_0 = U_0 \cos(\omega t)e_x$, where $\omega$ is the tidal frequency. However, the tide can also be chosen to be elliptical. The vertical coordinate $z$ is aligned with gravity and increases upward, while $e_x$ and $e_y$ represent the horizontal coordinate vectors.

The ocean is modelled as a rotating, incompressible and inviscid layer of finite depth $H_0$. The seafloor bathymetry is defined by its depth $h(x, y)$. We focus on isolated seamounts, meaning that $h$ tends to $-H_0$ in all directions. The ocean is stably stratified with mean density $\bar{\rho}_0$ and background buoyancy frequency $N$. For the results presented hereafter, the buoyancy frequency is chosen uniform. However, the method can be extended to more generic cases with vertically varying stratification $N(z)$ defined between $z = -H_0$ and $z = 0$, as done in two dimensions by Mathur et al. (2016). Under the $f$-plane approximation, the Coriolis frequency $f = 2\Omega \sin(\phi)$, with $\Omega$ the Earth's angular velocity is assumed spatially uniform for a given latitude $\phi$.

Velocity, pressure and buoyancy fluctuations ($\tilde{u}, \tilde{p}, \tilde{b}$) away from the background state verify the linear Boussinesq inviscid equations

$$ \tilde{u}_{,t} - f \tilde{v} = -\tilde{p}_{,x}, \quad \tilde{v}_{,t} + f \tilde{u} = -\tilde{p}_{,y}, \quad \tilde{w}_{,t} - \tilde{b} = -\tilde{p}_{,z}, \quad (2.1a,b,c) $$

$$ \tilde{b}_{,t} + N^2 \tilde{w} = 0 \quad \text{and} \quad \tilde{u}_{,x} + \tilde{v}_{,y} + \tilde{w}_{,z} = 0, \quad (2.1d,e) $$

where $\tilde{u} = (\tilde{u}, \tilde{v}, \tilde{w})$ and $\tilde{p}$ refers to the pressure scaled with the reference density $\bar{\rho}_0$ (with dimensions of velocity squared). Bold characters represent vectors. Commas followed by subscripts denote partial derivatives (e.g. $A_{,t} = \partial A / \partial t$).
We further assume that the waves are monochromatic with frequency $\omega$ and we introduce the complex value $A$, associated with the real quantity $\tilde{A}$ such that

$$ \tilde{A}(x, y, z, t) = \frac{1}{2}[A(x, y, z)e^{-i\omega t} + A^*(x, y, z)e^{i\omega t}] = \mathrm{Re}[A(x, y, z)e^{-i\omega t}], \quad (2.2) $$

where $A^*$ denotes complex conjugation. This assumption leads us to the internal-wave equation for the complex vertical velocity $w(x, y, z)$, namely,

$$ \left( \frac{N^2 - \omega^2}{\omega^2 - f^2} \right) \nabla_H^2 w - w_{,zz} = 0, \quad (2.3) $$

where $\nabla_H^2 = \partial_{xx} + \partial_{yy}$ is the horizontal Laplacian. We define

$$ \mu^{-1} = \sqrt{\frac{\omega^2 - f^2}{N^2 - \omega^2}}, \quad (2.4) $$

as the slope of the characteristics of the hyperbolic equation (2.3), which corresponds to the slope of the wave rays.

The problem has a rich parameter space, including length scales $H_0$, $L$, $\Lambda$; frequencies $\omega$, $f$, $N$; and the velocity $U_0$. Garrett & Kunze (2007) introduce five relevant non-dimensional numbers. The first three are taken as the frequency ratios

$$ \beta = f/\omega \quad \text{and} \quad \zeta = N/\omega, \quad (2.5) $$

and the height ratio

$$ \delta = \frac{\Lambda}{H_0}. \quad (2.6) $$

Of particular interest is the relative steepness $\epsilon$, defined as the ratio of the maximum slope of the topography in the direction of the barotropic tide to the internal-wave slope, namely,

$$ \epsilon = \mu \max (\partial_x h). \quad (2.7) $$

This parameter measures the criticality of the topography, with $\epsilon < 1$ (respectively $> 1$) corresponding to the subcritical (respectively supercritical) regime. Finally, the last non-dimensional number is defined as the ratio of the tidal excursion $U_0/\omega$ to the horizontal length of the topography $L$, namely,

$$ \alpha = \frac{U_0}{\omega L}. \quad (2.8) $$

The assumption of small excursion parameter $\alpha \ll 1$ allows us to neglect the advective term $U_0 \cdot \nabla$ in the Boussinesq equations (see e.g. Garrett & Kunze 2007 and Bell 1975a,b for the complete solution). Furthermore, if we assume that $u$ is smaller than, or similar to, $U_0$, the nonlinear terms $u \cdot \nabla$ can similarly be neglected. These conditions are verified when the topography is flat enough. As $\epsilon$ increases, beams with increasingly sharp gradients appear, reaching a singularity when $\epsilon = 1$. This results in the breakdown of the linear regime, which can lead to two types of processes: the generation of smaller-scale internal gravity waves (e.g. higher harmonics, overturning or shear instabilities, parametric subharmonic instability) with small group velocities that carry a small portion of the energy flux; and turbulent motions near critical slopes (Gáyen & Sarkar 2011). However, larger scales remain mostly unaffected by this breakdown. Therefore, it is still reasonable to explore solutions for steep topography while assuming linearity.
## 2.2. Boundary conditions

We consider an ocean with upper and lower horizontal boundaries. We disregard the motion of the free surface and treat it as a solid surface. Therefore, the velocity field must satisfy the non-penetration condition at both the bottom and the top boundaries, that is,

$$w|_{z=0} = 0 \quad \text{and} \quad w|_{z=h} = (U_0 + u_H) \cdot \nabla_H h, \qquad (2.9a,b)$$

where $u_H = (u, v)$ is the horizontal velocity, $r = (x, y)$ and $\nabla_H = (\partial_x, \partial_y)$. The horizontal boundary conditions are chosen so that there are only waves radiated outwards as $\|r\| \to \infty$.

In the WTA, the topography is assumed to be small and flat, in the sense that $\delta \ll 1$ and $\epsilon \ll 1$. The term $u_H \cdot \nabla_H h$, of order $\epsilon w$, can thus be neglected in the non-penetration condition at the seamount (2.9b) and the vertical velocity can be evaluated at the bottom ($z = -H_0$). Equation (2.9b) then becomes $w|_{z=-H_0} = U_0 \cdot \nabla_H h$. However, for steeper topographies, the horizontal velocity $u_H$ can be locally larger than, or of the same order as, the tidal velocity $U_0$, and this approximation is no longer valid. Here, we use the full boundary condition, in order to more accurately represent larger topographies whatever the values of the steepness.

## 2.3. Dimensionless modal equations for constant stratification

In the case of constant stratification $N$, the slope of the wave beams $\mu^{-1}$ is uniform along the depth. We introduce the following dimensionless length scales, velocities, pressure and buoyancy:

$$[\hat{x}, \hat{y}] = \frac{\pi [x, y]}{\mu H_0}, \quad \hat{z} = \frac{\pi z}{H_0}, \quad [\hat{u}, \hat{v}] = \frac{[u, v]}{U_0}, \quad \hat{w} = \frac{\mu w}{U_0}, \qquad (2.10)$$

$$\hat{p} = \frac{\pi p}{\omega U_0 H_0 \mu (1 - \beta^2)}, \quad \hat{b} = \frac{\mu \omega b}{N^2 U_0}. \qquad (2.11)$$

The asymmetric scaling between the horizontal and vertical directions allows us to rescale the problem in order to have wave beams with slope 1. In the following sections, all horizontal lengths, such as the topographic width and the horizontal mesh size, are scaled by $\mu H_0 / \pi$, while vertical lengths are scaled using $H_0 / \pi$.

We project the solution onto the vertical modes associated with the stratification $N$, namely,

$$[\hat{w}, \hat{b}](\hat{x}, \hat{y}, \hat{z}) = \sum [\hat{w}_m, \hat{b}_m](\hat{x}, \hat{y}) a_m(\hat{z}), \qquad (2.12a)$$

$$[\hat{u}, \hat{v}, \hat{p}](\hat{x}, \hat{y}, \hat{z}) = \sum \hat{k}_m^{-1} [\hat{u}_m, \hat{v}_m, \hat{p}_m](\hat{x}, \hat{y}) a'_m(\hat{z}). \qquad (2.12b)$$

With the scaling (2.10), the modal functions $a_m$ are determined by the following eigenproblem:

$$\frac{d^2 a_m}{dz^2} + m^2 a_m = 0 \quad \text{with} \quad a_m(0) = a_m(-\pi) = 0, \qquad (2.13)$$

where $m = \hat{k}_m = k_m(\mu H_0 / \pi)$ is the dimensionless horizontal wavenumber of mode $m$. This leads to: $a_m(\hat{z}) = \sin(m\hat{z})$. The first three vertical modes are represented in figure 1. The orthogonality conditions are

$$\int_{-\pi}^{0} \sin(n\hat{z}) \sin(m\hat{z}) \, d\hat{z} = \int_{-\pi}^{0} \cos(n\hat{z}) \cos(m\hat{z}) \, d\hat{z} = \begin{cases} \frac{\pi}{2} & \text{if } m = n, \\ 0 & \text{otherwise.} \end{cases} \qquad (2.14)$$
The Boussinesq equations (2.1) can be rewritten in their dimensionless form, dropping the hats for simplicity

$$u - i\beta v + i(1 - \beta^2)p_{,x} = 0, \quad v + i\beta u + i(1 - \beta^2)p_{,y} = 0, \quad w - ip_{,z} = 0, \qquad (2.15a,b,c)$$

$$-ib + w = 0 \quad \text{and} \quad u_{,x} + v_{,y} + w_{,z} = 0. \qquad (2.15d,e)$$

By projecting onto the vertical modes, i.e. multiplying (2.15a,b,e) by $\cos(mz)$ and (2.15c,d) by $\sin(mz)$ and integrating vertically between $-\pi$ and 0, we obtain the following modal equations:

$$u_m = m^{-1}(w_{m,x} + i\beta w_{m,y}), \quad v_m = m^{-1}(w_{m,y} - i\beta w_{m,x}), \qquad (2.16a,b)$$

$$p_m = im^{-1}w_m \quad \text{and} \quad b_m = -iw_m. \qquad (2.16c,d)$$

In the following sections, we use dimensionless quantities, but we remove the hats for simplicity, except in the figures' axes, legends and captions, where the scaling is fully detailed.

### 3. Solution based on Green's function approach

#### 3.1. Boundary integral equation

We derive an analytical expression for the solution of the problem using the boundary element method, which represents the topography as a layer of sources with unknown distribution $S(x, y)$. This method was first used by Pétrélis et al. (2006) to compute internal tide generation over 2-D topographies. However, their method, using vertical discretisation, was restricted to symmetric seamounts. Echeverri & Peacock (2010) extended it to arbitrary 2-D topographies by implementing horizontal discretisation. The method presented here is an adaptation of their method to 3-D seamounts.

The source term $S(x, y)$ is introduced as a forcing term for the vertical velocity wave equation

$$\nabla_H^2 w - w_{,zz} = S. \qquad (3.1)$$

We look for the Green's function $G$ or fundamental solution of problem (3.1) for a single source point located at $(\mathbf{r}', h(\mathbf{r}')) = (x', y', h(x', y'))$. This function is a solution of equation

$$\nabla_H^2 G - G_{,zz} = \delta^{(2)}(\mathbf{r} - \mathbf{r}')\delta^{(1)}(z - h(\mathbf{r}')), \qquad (3.2)$$

where $\delta^{(i)}$ is the Dirac delta function in $i$-dimensional Cartesian coordinates. We decompose $G$ into vertical modes, namely,

$$G(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) = \sum_{m=1}^{\infty} g_m(\mathbf{r}, \mathbf{r}') \sin(mz). \qquad (3.3)$$

With this decomposition, the solution verifies the non-penetration condition at the top and bottom boundaries ($z = -\pi$ and $z = 0$). Formally, the velocity, pressure and buoyancy fields are defined for the whole domain, both above and below the topography. The method thus implies that such a wave solution exists below the topography, which is not guaranteed, as discussed in Voisin (2021) (see § 6.1), especially in the case of supercritical seamounts. This is a consequence of the use of the indirect formulation of the boundary integral equation. However, we find that the direct formulation, as derived in Martin &
$$
\begin{align*}
\mathcal{G}(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) &= \sum_{m=1}^{\infty} \frac{\sin(mh(\mathbf{r}')) \sin(mz)}{2i\pi} \mathcal{H}_0^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|), \\
\mathcal{G}_u(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) &= -\sum_{m=1}^{\infty} \frac{(x-x') + i\beta(y-y')}{2i\pi \|\mathbf{r} - \mathbf{r}'\|} \sin(mh(\mathbf{r}')) \cos(mz) \mathcal{H}_1^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|), \\
\mathcal{G}_v(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) &= -\sum_{m=1}^{\infty} \frac{(y-y') - i\beta(x-x')}{2i\pi \|\mathbf{r} - \mathbf{r}'\|} \sin(mh(\mathbf{r}')) \cos(mz) \mathcal{H}_1^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|), \\
\mathcal{G}_p(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) &= \sum_{m=1}^{\infty} \frac{\sin(mh(\mathbf{r}')) \cos(mz)}{2m\pi} \mathcal{H}_0^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|), \\
\mathcal{G}_b(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) &= -\sum_{m=1}^{\infty} \frac{\sin(mh(\mathbf{r}')) \sin(mz)}{2\pi} \mathcal{H}_0^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|).
\end{align*}
$$

<caption>Table 1. Expression of the Green’s functions for the velocities w, u, v, the pressure p and the buoyancy b.</caption>

Llewellyn Smith (2012) and Voisin (2021) and which would avoid this assumption, cannot be easily applied in this case.

By projecting (3.2) on mode m, we obtain the following Helmholtz equation for the mode g_m:

$$
\nabla_H^2 g_m + m^2 g_m = \frac{2 \sin(mh(\mathbf{r}'))}{\pi} \delta^{(2)}(\mathbf{r} - \mathbf{r}'). \quad (3.4)
$$

The solution of this equation with outward radiation as $r = \|\mathbf{r}\| \to \infty$ is given by

$$
g_m(\mathbf{r}; \mathbf{r}') = \frac{\sin(mh(\mathbf{r}'))}{2i\pi} \mathcal{H}_0^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|), \quad (3.5)
$$

with $\mathcal{H}_0^{(1)}$ the Hankel function of the first type and order zero (Couto 2013). Thus, the 3-D Green’s function can be expressed as

$$
\mathcal{G}(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) = \sum_{m=1}^{\infty} \frac{\sin(mh(\mathbf{r}')) \sin(mz)}{2i\pi} \mathcal{H}_{0}^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|). \quad (3.6)
$$

The vertical velocity field can then be calculated at any point using the convolution

$$
w(x, y, z) = \iint_{\mathcal{R}} \mathcal{G}(\mathbf{r}, z; \mathbf{r}', h(\mathbf{r}')) S(\mathbf{r}') \, dx' dy' = \mathcal{G} * S, \quad (3.7)
$$

with $\mathcal{R}$ being the domain where the topography is defined, i.e. where the source $S$ is not zero. Moreover, knowing that the modal horizontal velocity and pressure can be expressed as function of $w_n$ (see (2.16)) and that $w_n = g_n * S$, the internal wavefields can be expressed as convolutions with $S$, that is,

$$
u = \mathcal{G}_u * S, \quad v = \mathcal{G}_v * S, \quad p = \mathcal{G}_p * S, \quad b = \mathcal{G}_b * S. \tag{3.8}
$$

The expressions of the functions $\mathcal{G}_u$, $\mathcal{G}_v$ and $\mathcal{G}_p$ are detailed in table 1.

Hence, the source distribution $S$ is the only unknown of the problem and has to verify
the condition at the seafloor, which involves both the vertical velocity field $w$ and the
horizontal velocity fields $u$ and $v$. Using (3.7) and (3.8), the non-penetration condition at
the topography

$$
U_0 \cdot \nabla_H h(\mathbf{r}) = w[\mathbf{r}, h(\mathbf{r})] - h_{,x}(\mathbf{r})u[\mathbf{r}, h(\mathbf{r})] - h_{,y}(\mathbf{r})v[\mathbf{r}, h(\mathbf{r})], \quad (3.9)
$$
can be rewritten as an integral equation problem

$$ U_0 \cdot \nabla_H h(\mathbf{r}) = \mathcal{F} * S = \iint_{\mathcal{R}} \mathcal{F}(\mathbf{r}, h(\mathbf{r}); \mathbf{r}', h(\mathbf{r}')) S(\mathbf{r}') \, dx'dy', \quad (3.10) $$

where

$$ \mathcal{F} = \mathcal{G} - h_{,x}\mathcal{G}_u - h_{,y}\mathcal{G}_v. \quad (3.11) $$

Equations (3.7) and (3.8) are valid anywhere outside the seamount surface. To write (3.10), it is necessary to ensure that the convolution integrals are continuous as the evaluation point $(x, y, z)$ approaches the topography, as discussed in Martin & Llewellyn Smith (2012). While a complete analytical proof of this continuity – such as the one outlined in Brebbia et al. (2012, Chapter 2.3) – eluded us in our configuration, we conducted a numerical analysis to assess the behaviour of the Green’s function $\mathcal{G}$ and its gradient components $\mathcal{G}_u$ and $\mathcal{G}_v$ near the surface. In subcritical cases, this analysis showed that all three functions can be bounded by integrable functions whose integrals over vanishingly small neighbourhoods around $(\mathbf{r}', z')$ tend to zero. This ensures the continuity of the integrals and validates the boundary integral formulation. However, the same method cannot be directly applied in supercritical cases, where a singular behaviour is expected along the critical wave beam (Le Dizès 2024). To the best of our knowledge, a rigorous proof of continuity in this regime remains an open problem. Despite this, we extended the method into the supercritical regime. However, we find that the BEM struggles to solve the problem for large values of $\epsilon$. We attribute these difficulties to numerical limitations due to the singular structure of the wave beams, which introduces large numerical errors, rather than to the continuity issue discussed here. However, we acknowledge that a formal proof of continuity is still needed and that the results of the BEM for supercritical cases need to be considered with caution.

To solve the integral problem in (3.10), we discretise the topography into square cells as shown in figure 1. These cells are denoted as $C_{ij} = [x_i - \Delta/2, x_i + \Delta/2] \times [y_j - \Delta/2, y_j + \Delta/2]$ in the following, where $\Delta$ is the cell side length (see figure 1). We then assume that $S$ is constant over a cell: $\forall \mathbf{r} \in C_{ij}, S(\mathbf{r}) = S_{ij}$. This allows us to write

$$ U_0 \cdot \nabla_H h(\mathbf{r}) = \sum_{ij} S_{ij} \iint_{C_{ij}} \mathcal{F}(\mathbf{r}, h(\mathbf{r}); \mathbf{r}', h(\mathbf{r}')) \, dx'dy'. \quad (3.12) $$

Integrating (3.12) over the cell $C_{kl}$, we obtain the following matrix equation:

$$ B = AS, \quad \text{with} \quad \left\{ \begin{aligned} B_{kl} &= \iint_{C_{kl}} U_0 \cdot \nabla_H h(x, y) \, dx \, dy, \\ A_{klij} &= \iint_{C_{kl}} \iint_{C_{ij}} \mathcal{F}(\mathbf{r}, h(\mathbf{r}); \mathbf{r}', h(\mathbf{r}')) \, dx \, dy \, dx' \, dy'. \end{aligned} \right. \quad (3.13) $$

The function $\mathcal{F}$ involves an infinite sum over the vertical modes. In practice, we truncate the sum at a maximum mode number $M_s$ to compute the matrix $A$ and solve the boundary integral equation (3.13). Moreover, the function $\mathcal{F}$ is a linear combination of the Hankel functions of the first type $\mathcal{H}_0^{(1)}$ and $\mathcal{H}_1^{(1)}$, which are both singular at $r = 0$. To handle this singularity, we analytically integrate their asymptotic forms.

Details on the computation of the matrix coefficients are given in Appendix A. Equation (3.13) is then solved numerically using the Generalized minimal residual solver from the PETSc library.
## 3.2. Wavefield and energetics

Once the boundary integral equation (3.10) is solved and the source distribution $S$ is computed, the convolutions (3.7)–(3.8) can be evaluated to reconstruct the velocity components $u, v, w$, the buoyancy $b$ and pressure $p$ at any point $(x, y, z)$ in the domain.

To discuss the solution in terms of energetic considerations, we consider the depth-integrated time-averaged energy density flux defined as

$$ \tilde{\mathbf{J}} = \int_{-\pi}^{0} \langle \tilde{p} \tilde{u} \rangle dz = \frac{1}{2} \int_{-\pi}^{0} \mathrm{Re}(p^* \mathbf{u}) dz = \sum_{m=1}^{\infty} \tilde{\mathbf{J}}_m, \quad (3.14) $$

where $\langle \cdot \rangle$ is the average over a time period and the modal flux $\tilde{\mathbf{J}}_m$ is given by

$$ \tilde{\mathbf{J}}_m = \mathrm{Re} \left[ \frac{\pi}{4} p_m^*(\mathbf{r}) \mathbf{u}_m(\mathbf{r}) \right]. \quad (3.15) $$

Another quantity of interest for internal tide generation is the energy conversion rate $\tilde{C}$, defined as the integral of the energy flux across a closed cylindrical boundary $\partial \Omega$ of height $\pi$ around the topography, namely,

$$ \tilde{C} = \iint_{\partial \Omega} \langle \tilde{p} \tilde{u} \rangle \cdot \mathbf{e}_r \, ds = \int_0^{2\pi} \tilde{\mathbf{J}} \cdot \mathbf{e}_r r \, d\theta = \sum_{m=1}^{\infty} \tilde{C}_m, \quad (3.16) $$

where $\mathbf{e}_r$ is the outward normal vector, $ds$ the dimensionless horizontal differential surface element, $(r, \theta)$ represent the horizontal cylindrical coordinates and

$$ \tilde{C}_m = \int_0^{2\pi} \tilde{\mathbf{J}}_m \cdot \mathbf{e}_r r \, d\theta \quad (3.17) $$

is the modal conversion rate. Far away from all sources ($r \gg 1$), the energy flux decays as $1/r$, as discussed hereafter, which means that the conversion rates will be independent of the distance. Here, the flux $\tilde{\mathbf{J}}$ and the conversion rate $\tilde{C}$ are scaled respectively by the quantities

$$ J_0 = \bar{\rho}_0 U_0^2 \mu \omega (1 - \beta^2) \frac{H_0^2}{\pi^2} \quad \text{and} \quad C_0 = \bar{\rho}_0 U_0^2 \mu^2 \omega (1 - \beta^2) \frac{H_0^3}{\pi^3}. \quad (3.18) $$

All of the convolutions (3.7)–(3.8), the energy flux (3.14) and the conversion rate (3.16) involve infinite sums over the vertical modes. In practice, these series are truncated at a finite mode number, $M_f$. To compute energy fluxes and for the analysis presented in § 6, we use a cutoff $M_f = M_f^*$, defined as the smallest mode such that the contribution $C_{M_f^*}$ represents less than 1% of the cumulative energy of all the lower modes, i.e. $C_{M_f^*} \le 0.01 \sum_{m=1}^{M_f^*} C_m$. The values of $M_f^*$, which were found to be independent of $\beta$, are listed in table 2 for the different topographies. They range from 2 to 12, increasing with $\epsilon$ and decreasing with $\delta$. However, it was observed that using $M_f = M_f^*$ is insufficient to accurately reconstruct the velocity field or clearly visualise the wave beams. This can be explained by the fact that modal amplitudes decay more rapidly with mode number for the energy flux than for the velocity components. As a result, for visualising the velocity field, we use an arbitrary value of $M_f = 20$ modes, which ensures that the wave beams are well defined.

Finally, as done by Pollmann et al. (2019), we are also interested in an asymptotic expression of the modal energy flux. This asymptotic expression is used in § 6 to compare
<table><caption>Table 2. Axisymmetric Gaussian cases, varying the criticality parameter $\epsilon$ for a fixed height ratio $\delta = 0.3$ (left) and varying $\delta$ for a fixed $\epsilon = 0.7$ (right). All cases were tested with five different values of $\beta = f/\omega = \{0; 0.1; 0.5; 0.7; 0.9\}$.</caption><thead><tr><th>ε</th><th>δ</th><th>M<sub>s</sub></th><th>Δ</th><th>M<sup>*</sup><sub>f</sub></th><th>ε</th><th>δ</th><th>M<sub>s</sub></th><th>Δ</th><th>M<sup>*</sup><sub>f</sub></th></tr></thead><tbody><tr><td>0.3</td><td>0.3</td><td>62</td><td>0.0508</td><td>2</td><td>0.7</td><td>0.1</td><td>300</td><td>0.0082</td><td>12</td></tr><tr><td>0.4</td><td>0.3</td><td>150</td><td>0.0264</td><td>2</td><td>0.7</td><td>0.2</td><td>240</td><td>0.0131</td><td>7</td></tr><tr><td>0.5</td><td>0.3</td><td>137</td><td>0.0229</td><td>3</td><td>0.7</td><td>0.3</td><td>160</td><td>0.0196</td><td>5</td></tr><tr><td>0.6</td><td>0.3</td><td>151</td><td>0.0208</td><td>4</td><td>0.7</td><td>0.4</td><td>132</td><td>0.0238</td><td>4</td></tr><tr><td>0.7</td><td>0.3</td><td>160</td><td>0.0196</td><td>5</td><td>0.7</td><td>0.5</td><td>115</td><td>0.0272</td><td>3</td></tr><tr><td>0.8</td><td>0.3</td><td>200</td><td>0.0286</td><td>6</td><td>0.7</td><td>0.6</td><td>96</td><td>0.0327</td><td>2</td></tr><tr><td>0.9</td><td>0.3</td><td>200</td><td>0.0254</td><td>7</td><td>0.7</td><td>0.7</td><td>96</td><td>0.0327</td><td>2</td></tr><tr><td>1.1</td><td>0.3</td><td>750</td><td>0.0208</td><td>9</td><td>0.7</td><td>0.8</td><td>200</td><td>0.0697</td><td>2</td></tr></tbody></table>

the results obtained by our method and the WTA. First, we write the modal horizontal velocity in the polar coordinates $(r, \theta)$

$$ \mathbf{u}_m = -i \begin{pmatrix} p_{m,x} + i\beta p_{m,y} \\ p_{m,y} - i\beta p_{m,x} \end{pmatrix}_{(\mathbf{e}_x, \mathbf{e}_y)} = -i \begin{pmatrix} p_{m,r} + i\beta p_{m,\theta}/r \\ p_{m,\theta}/r - i\beta p_{m,r} \end{pmatrix}_{(\mathbf{e}_r, \mathbf{e}_\theta)} . \quad (3.19) $$

Given that $p_m = im^{-1}g_m * S$, we can express the modal pressure as

$$ p_m = \frac{1}{2\pi m} \iint_{\mathbb{R}} \sigma_m(\mathbf{r}') \mathcal{H}_0^{(1)}(m \|\mathbf{r} - \mathbf{r}'\|) \, dx'dy', \quad (3.20) $$

with $\sigma_m$ a modal forcing term defined as

$$ \sigma_m(\mathbf{r}) = \sin(mh(\mathbf{r}))S(\mathbf{r}). \quad (3.21) $$

Using the asymptotic expansion of the Hankel functions, similar to what is done by Pollmann et al. (2019) (see their equation (2.23)) but at a higher order of approximation, we can express the modal pressure at a point $(r, \theta)$ far away from all sources as

$$ p_m(\mathbf{r}) \sim \frac{1}{2\pi m} \sqrt{\frac{2}{\pi m r}} e^{i(mr - \pi/4)} \left[ \left(1 - \frac{i}{8mr}\right) f_m(\theta) + i \frac{g_m(\theta)}{2r} \right], \quad (3.22) $$

where

$$ f_m(\theta) = \iint \sigma_m(\mathbf{r}') e^{-ime_r \cdot \mathbf{r}'} dx'dy' = \mathfrak{F}[\sigma_m](me_r) = \mathfrak{F}[\sigma_m](m \cos \theta, m \sin \theta), \quad (3.23) $$

$$ g_m(\theta) = \iint (-ie_r) \cdot \mathbf{r}' \sigma_m(\mathbf{r}') e^{-ime_r \cdot \mathbf{r}'} dx'dy' = (\partial_k \mathfrak{F}[\sigma_m])(me_r), \quad (3.24) $$

where $\mathfrak{F}[\sigma_m]$ is the Fourier transform of the forcing term $\sigma_m$ and $\partial_k \mathfrak{F}[\sigma_m]$ is the derivative of the Fourier transform with respect to the wavevector norm $k$. By definition, $f_m$ and $g_m$ are solely functions of the polar coordinate $\theta$ via the unit vector $e_r$ and do not depend on the radius $r$, but vary depending on the nature of the topography and the ratio of Coriolis and tidal frequencies $\beta = f/\omega$.

We use a similar approach to express the derivatives of the modal pressure. Details of the calculation are given in Appendix B. This leads to an approximation at $O(r^{-2})$ of the
<figure><img src="image_4.png" /><figcaption>Figure 2. Source distribution for an axisymmetric Gaussian of parameters $\delta = 0.3$ and $\epsilon = 0.7$, solved for $M_s = 160$ modes: magnitude (a) and phase (b).</figcaption></figure>

modal energy flux

$$ \tilde{J}_m = \frac{1}{8\pi^2 m^2 r} \left\{ \left( |f_m|^2(\theta) + \beta \frac{\text{Re}[f_m^*(\theta) f_m'(\theta)]}{mr} + \frac{\text{Im}[f_m(\theta) g_m^*(\theta)]}{r} \right) e_r + \left( \beta \frac{|f_m|^2(\theta)}{2mr} + \frac{\text{Im}[f_m^*(\theta) f_m'(\theta)]}{mr} \right) e_\theta \right\}. \quad (3.25) $$

This expression shows the leading-order $1/r$ decay of the energy flux in the far field. The azimuthal component decreases as $1/r^2$ and, therefore, far from the generation site, the flux is mainly radial.

## 4. Internal tide radiated by Gaussian seamounts

### 4.1. Examples of subcritical and supercritical seamounts

To validate our boundary integral method and illustrate how topographic criticality influences both the internal tide field and the numerical performance, we investigate the internal tide generated by a unidirectional $M_2$ tide along $e_x$ over two references topographies, one subcritical and one supercritical. In both cases, the seamount is modelled as an axisymmetric Gaussian topography, defined as:

$$ h(x, y) = \begin{cases} -\pi + \pi \delta e^{-\frac{|r|^2}{2L^2}} & \text{if } |x|, |y| \le 3L, \\ -\pi & \text{otherwise.} \end{cases} \quad (4.1) $$

We focus first on the subcritical topography with height ratio $\delta = 0.3$ and criticality $\epsilon = 0.7$. This yields a dimensionless topographic width of $L = \pi \delta / (\epsilon \sqrt{e}) \approx 0.82$. Rotational effects are initially neglected by setting $\beta = 0$. Using dimensional parameters relevant to the oceanic context, for example an ocean depth of $H_0 = 3$ km, a tidal frequency $\omega = 1.4 \times 10^{-4}$ s$^{-1}$ and a buoyancy frequency $N = 2 \times 10^{-3}$ s$^{-1}$, this would correspond to a dimensional seamount of height 900 m and width 11 113 m.

The topography spans a mesh of 250 × 250 points, resulting in a horizontal mesh size of $\Delta \approx 0.020$. Figure 2 displays the source distribution, computed with $M_s = 160$, with magnitude and phase represented, respectively, in panels (a) and (b). The magnitude of
<figure><img src="image_5.png" /><figcaption>Figure 3. Isocontours of internal-wave vertical velocity field $\mu w/U_0 = \pm 0.37$ generated by a unidirectional tide over a subcritical Gaussian topography with $\delta = 0.3$ and $\epsilon = 0.7$, reconstructed with $M_f = 20$ modes. The blue and yellow panels correspond, respectively, to the vertical and horizontal slices shown in figure 4.</figcaption></figure>

the source distribution varies continuously between values of the order [0.1 – 10] except for a specific line at $x = 0$ and along specific circular ones, where very weak amplitudes are observed ($\|S\| < 10^{-3}$). The regions separated by these lines display changes in sign of the phase distribution. Increasing the maximum number of modes to compute $S$ leads to the appearance of finer length scales in the source distribution but does not cause significant differences in the wavefield presented below.

Figure 3 shows the isocontours $w = \pm 0.37$ of the wave’s vertical velocity, obtained from (3.7) using the previously shown source distribution. As discussed previously, the wavefield is reconstructed with the first 20 modes. Notably, we observe the conical structure of the wave beams, as expected for 3-D internal-wave generation. The beams then reflect at the top and bottom boundaries as they propagate away from the topography. Figure 4(a) shows horizontal and vertical slices of the vertical velocity field, at $z = -1$ and $y = 0$. We note the symmetry about the x-axis, consistent with the unidirectional nature of the tide and the axisymmetry of the topography. The amplitude of the vertical velocity is maximised along the direction of the tide. Moving azimuthally away from this axis, the amplitude progressively decays, reaching zero on the y-axis, perpendicular to the generating tide.

The method also allows us to compute a solution for the internal tide field for supercritical cases. However, as mentioned in the previous section, the solver fails to converge for large values of $\epsilon$. Because of this, we restrict our analysis to $\epsilon < 1.1$. Larger values may require numerical adaptations to better handle singularities at critical points.

Figure 4(b) shows the horizontal and vertical slices of the vertical velocity field for the case $\delta = 0.3$ and $\epsilon = 1.1$, reconstructed using the first 20 modes. The horizontal slices shown in the top panels appear qualitatively similar for both subcritical and supercritical topographies. In both cases, we observe radially radiating rings with azimuthal amplitude modulation. However, the vertical slices in the bottom panels reveal that the wave beams are narrower in the supercritical case and reach a maximum amplitude of 3.9 (saturated in the colour bar) roughly 3.5 times larger than in the subcritical case. This large amplitude comes from the singular nature of the wave beam in supercritical cases. Although, here, the use of a finite number of modes ($M_f$) regularises the solution and smooths out
<figure><img src="image_6.png" /></figure>

Figure 4. Internal-wave vertical velocity field $\mu w/U_0$ reconstructed with $M_f = 20$ for (a) a subcritical topography ($\delta = 0.3, \epsilon = 0.7$) and (b) a supercritical topography ($\delta = 0.3, \epsilon = 1.1$). The top orange panels represents horizontal slices at $\pi z/H_0 = -1.0$, and the bottom blue panels vertical slices at $\pi y/(\mu H_0) = 0$. For the case (b), the colour bar is saturated at $\pm 2.0$, while the true extrema are $\pm 3.9$.

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

Figure 5. Energy conversion rate for the first 20 modes for a subcritical (black, $\epsilon = 0.7$) and supercritical topography (grey, $\epsilon = 1.1$): (a) fraction $C_m$ going into each mode $m$, normalised by $C_0$, with decay rates represented by the dashed lines: 3 m$^{-6}$ in black and 500 m$^{-3}$ in grey; (b) cumulative sum up to mode $M_f$ scaled by model predictions: $C_{\text{model}} = C^{\text{WTA}}$ (subcritical) and $C_{\text{model}} = C^{\text{Voisin}}$ (supercritical).

the singularity. Additionally, we observe in figure 4(b) an oscillating signal with small wavelength near the centre of the topography in (x, y) = (0, 0), which we interpret as a Gibbs phenomenon.

For further analysis, we compute for both topographies the modal energy conversion rate $C_m$ going into mode $m$, as shown in figure 5(a). In both regimes, the energy is concentrated in the first 10 modes and decreases with mode number from $m = 2$ onward.
In the subcritical case, mode 1 dominates. However, for the supercritical case, this is no longer the case and the energy content appears to decrease more slowly with mode number than for the subcritical case. The black and grey dashed lines are indication for the rates of decay, which are proportional to $m^{-6}$ for the subcritical case and $m^{-3}$ for the supercritical case. Similar results for subcritical and supercritical cases were previously observed in two dimensions by Papoutsellis et al. (2023) (see their figure 3a). Additionally, we compute the global conversion rate defined as the cumulative sum of the modal components $C_m$ up to mode $m = M_f$. Figure 5(b) displays this quantity scaled by a reference value $C_{model}$ as a function of the maximum mode number $M_f$. For the subcritical case, we use the WTA estimate as a reference for the axisymmetric Gaussian in a finite-depth ocean from Llewellyn Smith & Young (2002), computed over the first 20 modes

$$C^{WTA} = \frac{\pi^3}{2} \delta^2 L^4 \sum_{m=1}^{20} C_m^{WTA}, \quad \text{with} \quad C_m^{WTA} = m^2 \exp(-m^2 \hat{L}^2), \quad (4.2)$$

where $C^{WTA}$ and $C_m^{WTA}$ are scaled by $C_0$ and $L$ by $\mu H_0 / \pi$. Using $\delta = 0.3$ and $\epsilon = 0.7$ (i.e. $L \approx 0.82$), this yields $C^{WTA} \approx 0.505$. For the supercritical seamount, we compare our estimate with the value from Voisin (2024) for an equivalent spheroid in an infinite ocean. As discussed by this author, the equivalent spheroid has the same height $\Gamma$ and a radius given by $d = 3/(2\sqrt{\pi}) \times L \approx 0.846 L$. We compute the conversion rate using (5.6) and (5.7) from that reference with $\gamma'_v = \mu \Lambda / d \approx 2.143$, which yields $C^{Voisin} \approx 0.104$. In both cases, the conversion increases with $M_f$ and reaches a plateau at a certain threshold: $M_f^* = 5$ (subcritical) and $M_f^* = 9$ (supercritical). For the subcritical case, the WTA provides a good estimate, with a relative error around 5% that can be attributed to convergence errors (see Appendix C). In the supercritical case, however, Voisin's estimate for an equivalent spheroid underpredicts the conversion by a factor of 2.5. Finite-depth effects are unlikely to explain this discrepancy, as the topography is small enough for the reflected beams not to re-interact with the topography. Instead, we believe that the difference in shape between the Gaussian and the ellipsoid accounts for the higher conversion rate observed with our model.

We verified that the results presented above are converged with respect to mesh size, mode truncation and domain extent, with a relative error lower than 5%. A detailed convergence analysis is provided in Appendix C. We find that the convergence behaviour with mesh size is similar in both regimes. However, the number of modes $M_s$ required to reach convergence for the supercritical case is significantly higher. Specifically, the convergence rate scales approximately as $M_s^{-3/2}$ in the subcritical case and as $M_s^{-3/4}$ in the supercritical case, indicating slower convergence in the latter. This increased modal requirement might partially explain the numerical difficulties encountered for larger values of $\epsilon$, as resolving the singular beam structure necessitates a larger number of vertical modes.

## 4.2. Energy conversion rate

We now evaluate the energy conversion rate $C$, computed using the first 20 modes, for a series of Gaussian topographies with different values of height ratio $\delta$, relative steepness $\epsilon$ and ratio $\beta = f/\omega$.

For each pair $(\delta, \epsilon)$, we examine multiple values of $\beta = f/\omega$ ranging from 0 to 1. Negative values of $\beta$, corresponding to latitudes in the southern hemisphere, were also tested but display a north/south symmetry with the corresponding positive $\beta$ and are therefore not shown here. To maintain a fixed criticality $\epsilon$ for each $\beta$, the characteristic
<figure><img src="image_8.png" /><figcaption>Figure 6. Evolution of the dimensionless energy conversion rate $C/C_0$ (computed using the first 20 modes) with the criticality parameter $\epsilon$ for $\delta = 0.3$ (a) and with the height ratio $\delta$ for $\epsilon = 0.7$ (b). The black dashed lines correspond to the value given by the WTA prediction for a Gaussian topography.</figcaption></figure>

width $L$ of the Gaussian is adapted according to $L = \mu H_0 \delta / (\epsilon \sqrt{\epsilon})$. All tested cases are summarised in table 2. The grid cell size $\Delta$ and maximum mode number $M_s$ were selected to balance computational efficiency with numerical accuracy, following the guidelines discussed in Appendix C. For each configuration, we systematically test multiple values $M_s$ and retain the lowest value that ensures convergence of the total energy conversion rate to within 3% relative error. The mesh resolution, however, is limited by numerical constraints. The associated error is estimated by extrapolation from the detailed convergence studies in Appendix C. Across all cases considered, the coarsest mesh corresponds to $\Delta = 0.0697$, which yields a maximum relative error on the conversion rate of approximately 5%.

Figure 6(a) shows the variation of $C$ with $\epsilon$ for a Gaussian seamount with a given height $\delta = 0.3$, the colours of the symbol reflecting the values of the ratio $\beta = f/\omega$. The black solid line corresponds to the WTA estimate, defined in (4.2). As discussed in Papoutsellis et al. (2023) (see their figure 6), energy conversion rates are in good agreement with the WTA when $0.5 \le \epsilon \le 0.9$ for a Gaussian topography with $\delta = 0.3$. When the topography becomes supercritical $\epsilon \ge 1$, both our model and theirs deviate smoothly from the WTA. However, a stronger and more abrupt deviation is observed for $\epsilon \le 0.4$. This corresponds to configurations where a beam emanates from one side of the seamount, propagates up and towards the centre, reflects at the top and interacts with the other side of the seamount.

Similarly, figure 6(b) shows the evolution of $C$ with $\delta$ for a Gaussian with $\epsilon = 0.7$ and varying values of $\beta$. As predicted by Papoutsellis et al. (2023), we observe a strong deviation from the WTA estimate for $\delta \ge 0.6$ for a fixed $\epsilon = 0.7$. The deviation is related to instances of weakly radiating topographies, where the radiated flux is several orders of magnitude smaller than the WTA estimate, for example the cases ($\delta = 0.7, \epsilon = 0.7$) and ($\delta = 0.8, \epsilon = 0.7$). This occurs when opposite sides of the seamount emit waves that interact with each other in a destructive manner (see figure 7). This phenomenon was previously described in two dimensions (Pétrélis et al. 2006; Maas 2011; Papoutsellis et al. 2023; Geoffroy et al. 2024), but to our knowledge, it has not been previously documented in 3-D cases.

The Coriolis frequency appears to have only a minor effect on the conversion rate for topographies well described by the WTA. However, this effect is non-monotonic. The conversion rate increases with $\beta = f/\omega$ for small $\epsilon$ and large $\delta$, while it decreases
<figure><img src="image_9.png" /><figcaption>Figure 7. Weakly radiating topography ($\delta = 0.7$, $\epsilon = 0.7$): (a) vertical slice of vertical velocity (with $M_f = 20$) and (b) horizontal mapping of the vertically integrated energy flux $\tilde{J}$ reconstructed with $M_f^* = 2$ modes.</figcaption></figure>

with $\beta$ for larger $\epsilon$ and smaller $\delta$. This non-trivial effect of $\beta$ on the conversion rate was observed by Baines (2007, see their figure 7) for a cylindrical pillbox of height $h_m$ and radius $a$, in the limit of high $a\omega/NH_0$. However, this limit cannot be directly accessed in our case, as the pillbox has an infinite slope, which gives $\epsilon = \infty$. For a Gaussian, taking the limit $\epsilon \to \infty$ for a given $\delta$ leads to $a = L = 0$. Moreover, it is interesting to note that the weak radiation phenomenon seems robust to variations of the Coriolis frequency.

In the rest of this study, we focus on the impact of $\beta$ on the wavefield patterns and on the directivity of the energy flux. To ensure reliable analysis, we restrict our study to radiative topographies, as weakly radiating cases exhibit energy flux levels comparable to numerical noise, leading to significant uncertainties.

## 5. Influence of the Coriolis frequency

As discussed in the previous section, the Coriolis frequency only has a limited effect on the overall energy conversion rate. However, its impact on the velocity field patterns as well as on the spatial distribution and direction of the energy flux, is significant.

### 5.1. Velocity patterns

Figure 8 presents horizontal slices of the velocity components $w$ (a), $u$ (b) and $v$ (c) at depth $z = -1$, computed for an axisymmetric Gaussian seamount with ($\delta = 0.3$, $\epsilon = 0.7$). These fields are reconstructed by summing the contributions of the first 20 vertical modes. The panels from left to right show the effect of increasing the non-dimensional Coriolis parameter $\beta$ from 0 to 0.9. In the non-rotating case ($\beta = 0$), the vertical velocity $w$ in the upper-left panel exhibits a pattern of concentric waves radiating outward from the topography, with an anti-symmetry about the x-axis, as described in the previous section. However, as $\beta$ increases, a distinct spiralling pattern emerges in the wavefield. This effect is most pronounced in the horizontal velocity components $u$ and $v$. Similar spiral structures were previously identified by Munroe & Lamb (2005, see their figure 3) as a consequence of the rotation, for a supercritical asymmetric Gaussian with horizontal aspect ratio $L_x/L_y = 0.7$, where $L_x$ and $L_y$ are the characteristic length scales in the x- and y-directions.
<figure><img src="image_10.png" /><figcaption>Figure 8. Horizontal slices at $\pi z/H_0 = -1.0$ of the dimensionless vertical (a) and horizontal (b,c) velocity field generated by a uni-directional tide on a axisymmetric Gaussian with $\delta = 0.3$ and $\epsilon = 0.7$ for different values of $\beta = 0$; 0.5; 0.9, increasing from left to right. The wavefield is reconstructed using $M_f = 20$ modes.</figcaption></figure>

## 5.2. Energy fluxes

Figure 9(a) displays the depth-integrated energy fluxes generated for the same Gaussian topography ($\delta = 0.3$, $\epsilon = 0.7$) and increasing values of $\beta = 0$, 0.5 and 0.9. The flux is computed using the first five modes, which together account for approximately 99% of the total flux (see table 2). When $\beta = 0$, the energy flux is symmetric about the x-axis and propagates radially outward in two dominant lobes on each side of the topography for positive and negative x. The radiation is strongest along the x-axis, while the y-axis exhibits minimal flux. As $\beta$ increases, the flux field becomes increasingly asymmetric. An azimuthal component emerges, strongest near the seamount and decaying with distance. This component also increases with $\beta$, resulting in a flux that is no longer purely radial. Moreover, the dominant lobes are progressively rotated counter-clockwise with increasing $\beta$. This means that the most intense fluxes are no longer aligned with the direction of the barotropic tide, i.e. the x-axis, but appear slightly deviated.

To quantify this behaviour, we introduce the deviation angle $\theta_{max}(r)$, defined as the azimuthal angle at which the radial component of the energy flux reaches its maximum for a given distance $r$ from the seamount. The locations $(r, \theta_{max}(r))$, tracked for each lobe,
<figure><img src="image_11.png" /><figcaption>Figure 9. Energy density flux generated by a uni-directional tide on an axisymmetric Gaussian with $\delta = 0.3$ and $\epsilon = 0.7$ using our model (a) and the WTA (b) for different values of $\beta$, increasing from left to right. The hatched region corresponds to the region closed to the seamount, for which the asymptotic expression is not valid. The flux is calculated by summation of the first modes up until $M_f^* = 5$. The dashed lines represent the locations of the maximum radial flux for a given radius.</figcaption></figure>

are represented by the dashed lines in figure 9. In the following, we only focus on the right side of the topography ($x > 0$), for which $\theta_{max} \in [-\pi/2, \pi/2]$. The energy fluxes maps and specific deviations obtained with our model are compared with the WTA prediction in the following section, using the asymptotic expression of the energy flux developed in § 3.2.

## 5.3. Revisiting the WTA case

Previous studies by Llewellyn Smith & Young (2002) and later Pollmann et al. (2019) derived analytical solutions for internal tide generation in the hydrostatic limit under the WTA (see equation 2.19 of Pollmann et al. 2019). By extending their results to the non-hydrostatic case and using our dimensionless values, we can apply their approach to express the far-field energy flux as done in § 3.2, using the modified forcing term (see (3.21)) defined by

$$ \sigma_m^{WTA} = (-1)^m m U_0 \cdot \nabla_H h(\mathbf{r}). \qquad (5.1) $$

Figure 10 compares their mode-1 forcing terms $\sigma_1^{WTA}$ with $\sigma_1$ from our model, obtained for a unidirectional tide along $e_x$ over an axisymmetric Gaussian topography with $\delta = 0.3$, $\epsilon = 0.7$ and increasing values of $\beta = 0, 0.5$ and $0.9$. As is visible in the upper panels, the two models exhibit very different forcing terms with larger values of the amplitude in our case, and a different spatial distribution. In the WTA, the amplitude peaks along the flanks of the seamount where the slopes are the steepest. In contrast, our model predicts stronger amplitude near the summit, with concentric modulations and rings of zero forcing. The lower panels represent the phase distribution. The WTA produces two phase values, 0 and $\pi$, representing an anti-symmetry between the left and right sides of the topography. Our model instead exhibits a rich phase structure displaying concentric
<figure><img src="image_12.png" /><figcaption>Figure 10. Forcing term $\sigma_1$ for the WTA (a) and our model for an axisymmetric Gaussian topography with $\delta = 0.3$, $\epsilon = 0.7$ and different values of $\beta$: (b) $\beta = 0$; (c) $\beta = 0.5$ and (d) $\beta = 0.9$. The top panels show the amplitudes and the bottom panels the phase distribution.</figcaption></figure>

rings that are correlated to the modulations observed for the amplitude. For $\beta = 0$, each concentric ring displays as well a $\pi$-shift across the y-axis. Increasing $\beta$ disrupts this anti-symmetry. In particular, the phase distribution features partial azimuthal shifts that depend on the distance from the centre of the topography. This is not accounted for by the WTA, for which the source term does not depend on $\beta$.

From the WTA forcing term in (5.1), we derive analytical expressions for the functions

$$f_m^{WTA}(\theta) = \zeta_m \cos(\theta) \quad \text{and} \quad g_m^{WTA}(\theta) = \xi_m \cos(\theta), \qquad (5.2)$$

with $\zeta_m = -i(-1)^m m^2 \mathfrak{F}[h](me_r)$ and $\xi_m = -i(-1)^m m (1 - m^2 L^2) \mathfrak{F}[h](me_r)$. These coefficients are imaginary constants for a given mode $m$ in the case of an axisymmetric topography. In practice, the topography is defined as a Gaussian seamount truncated using a square window of side length $6L$ (see (4.1)). The topography is thus not rigorously axisymmetric. However, we expect the effects of the truncation to induce an angular dependence of the coefficients only for the higher modes. In this context, the WTA modal flux is given by

$$\tilde{J}_m^{WTA} = \frac{|\zeta_m|^2}{8\pi^2 m^2 r} \cos^2(\theta) \left\{ \left(1 - \frac{\beta \tan(\theta)}{mr}\right) e_r + \frac{\beta}{2mr} e_\theta \right\}. \qquad (5.3)$$

An interesting observation is that, at $O(1/r^2)$, the flux is not purely radial but has an azimuthal component proportional to $\beta$.

Figure 9(b) displays the total energy flux, calculated using (5.3) for the first five modes. The expression is applied only far enough away from the seamount centre, i.e. outside of the hatched region. While the amplitudes appear similar to that of our model, the direction and spatial distribution differ significantly. In particular, the direction of maximum radial flux, materialised by the dashed lines, occurs at negative angles. From (5.3), we derive the corresponding deviation angle

$$\theta_{max}^{WTA}(r) \sim -\frac{\beta}{2\chi r}, \quad \text{with} \quad \chi = \frac{\sum m^{-2} |\zeta_m|^2}{\sum m^{-3} |\zeta_m|^2}. \qquad (5.4)$$
<figure><img src="image_13.png" /><figcaption>Figure 11. (a) Radial flux along the azimuthal coordinate $\theta$ at different values of $r$, and different values of $\beta$. The solid lines represent (6.1), summing up to $M_r^* = 5$. (b) Direction of the maximum radial flux $\theta_{max}$ as a function of the inverse radii $1/r$ (symbols '+'), with the least-squares linear fit materialised by the dashed lines. The symbols '•' are the deviations $\theta_{\infty}$ given by the asymptotic expression.</figcaption></figure>

This implies that the angular deviation in the WTA decays like $1/r$ and vanishes in the far field ($r \to \infty$).

These effects were not accounted for in the WTA framework used by Pollmann et al. (2019), who considered only the leading-order $O(1/r)$ terms for which the flux remains radial and symmetric about the x-axis.

# 6. Deviation of the energy fluxes

## 6.1. Direction of the maximum radial flux

To characterise the deviation of the energy fluxes predicted by our model, we investigate how the angles $\theta_{max}$ at which the radial flux is maximum, evolve with the distance $r$ and the Coriolis parameter $\beta$.

To do so, we plot in figure 11(a) the radial component of the energy flux as a function of the azimuthal angle $\theta$ for different $r$ and $\beta$. The fluxes are here computed using the first five modes via (3.14). For each configuration, we identify the angle $\theta_{max}$ of maximum radial flux and plot the values in figure 11(b) as a function of the inverse radius $1/r$. The near linear trend, materialised by the dashed lines using a least squares, indicates that the deviations increase with distance. However, contrary to the WTA model, the angles $\theta_{max}$ in our model are positive and increase with $\beta$. Moreover, extrapolation shows that the deviation does not vanish as $r \to \infty$, but instead converges to a finite positive value $\theta_{\infty}$.

To estimate this far-field limit $\theta_{\infty}$ accurately, we use the asymptotic expression from (3.25). At $O(1/r)$, the energy flux is radial and its radial component can be approximated by

$$ (r \tilde{J}_r)_\infty = \frac{1}{8\pi^2} \left( \sum_m \frac{|f_m|^2(\theta)}{m^2} \right). \qquad (6.1) $$

The azimuthal dependence is given by the functions $|f_m|^2(\theta) = |\mathfrak{F}[\sigma_m](n \cos \theta, n \sin \theta)|^2$. We compute $(r \tilde{J}_r)_\infty$ numerically up to $M_r^*$ modes. This ensures that we accurately represent the energy content of the internal waves, but avoid introducing spurious noise due to higher modes. Moreover, we extend the source terms $\sigma_m$ outside of the domain where the topography is defined, using zero padding up to $60L$, before computing the
<figure><img src="image_14.png" /><figcaption>Figure 12. Evolution of the deviation angles $\theta_{\infty}$ for which the radial flux at infinity is maximal, with $\beta$, $\epsilon$ and $\delta$. (a) The height ratio of the topography is kept at $\delta = 0.3$ and the criticality parameter $\epsilon$ is varied; (b) $\delta$ is varied, while $\epsilon$ is kept constant at $\epsilon = 0.7$.</figcaption></figure>

Fourier transform $\mathfrak{F}[\sigma_m](me_r)$. We then interpolate the function $(r\tilde{J}_r)_{\infty}(\theta)$ between $-\pi/2$ and $\pi/2$ to locate its maximum. This function is plotted for a few values of $\beta$ in figure 11(a) using solid lines. Finally, we measure the far-field deviation angles $\theta_{\infty}$ for which the maximum is reached and plot the values as the ‘•’ symbols along the y-axis in figure 11(b).

The results confirm that our model predicts persistent positive deviation angles in the far field, which increase with $\beta$, in contrast to the WTA for which the deviation vanished in the far field. This discrepancy is a direct consequence of the introduction of nonlinear terms in the boundary condition at the seafloor (2.9b). Adding the horizontal velocity components $u$ and $v$ breaks the symmetry and leads to the observed asymmetric spatial distribution in the flux field. As the Coriolis force increases, the radiated energy pattern rotates counter-clockwise, shifting the direction of maximum radiation away from the axis of the barotropic tide. This effect has also been observed by Baines (2007) (see figure 17b) for the pillbox topography. This author found that, when $\beta = 0$, the direction of maximum flux aligns with that of the barotropic tide ($\theta = 0$), with minima at $\theta = \pm \pi/2$. However, as $\beta$ increases, this direction moves to positive angles, and for mode-1 waves, it may reach values in excess of 1 rad close to the limit $\beta = 1$.

## 6.2. Parametric dependence of the far-field deviation angle

Figure 12 explores how the far-field deviation angle $\theta_{\infty}$ evolves with $\beta$ and the topography’s geometric properties ($\delta$, $\epsilon$). Panel (a) presents $\theta_{\infty}$ as a function of $\beta$ and $\epsilon$ for a fixed $\delta = 0.3$. The deviation increases linearly with $\beta$, with a stronger effect for seamounts with larger values of steepness $\epsilon$. For smaller values of $\epsilon$, the deviation converges towards zero, consistent with the WTA prediction. In the inset of figure 12(a), the data collapse on a single linear trend

$$ \theta_{\infty} \sim \frac{\beta \epsilon}{2}. \qquad (6.2) $$

As shown in figure 12(b), there is no visible influence of the height ratio $\delta$ on the deviation angle for values smaller than 0.5, and the results can be fitted with the linear trend given by (6.2). For higher values of $\delta$, namely, 0.6, 0.7 and 0.8, the left and right sides of the topography interact with each other, giving a more complex spatial distribution and very small values of the conversion rate, as discussed in § 4. We leave these cases to future investigations.
It is important to note that this result is specific to axisymmetric topographies. For elongated topographies (not shown here), we verified that the deviation of the energy flux with the Coriolis frequency is significantly reduced and the orientation between the tidal forcing and the topography's major axis governs the angular dependence of the energy flux radiated, as discussed in § 1.

## 7. Conclusion and perspectives

This study presents a novel 3-D analytical model for internal tide generation over isolated topographies, extending previous theoretical approaches based on the WTA (Llewellyn Smith & Young 2002; Pollmann et al. 2019) by removing restrictions on the height and the slope of the topography. Our method generalises the 2-D boundary element method of Echeverri & Peacock (2010) to 3-D cases. It assumes linearly propagating inviscid waves that verify the Boussinesq approximation and uses the decomposition into vertical modes to calculate the wavefield. The solution of the problem is obtained by solving the fully nonlinear boundary condition at the seafloor (recall 2.9b), which is a noticeable difference from past studies. The model is implemented using a discrete boundary element method on a regular Cartesian grid. It provides a full 3-D solution to the internal tide generation problem for arbitrary seamounts. In particular, the method offers a description of the various components of the wavefield, as well as energy conversion for both subcritical and supercritical seamounts. However, the supercritical regime presents significant numerical challenges. In this case, the internal-wave beams are singular and a large number of vertical modes are required to accurately capture the velocity fields, thus increasing computational cost. Moreover, for large values of the criticality parameter $\epsilon$, the solver fails to converge, and no solution is obtained. As a result, our analysis is restricted to cases with $\epsilon \le 1.1$. We attribute these difficulties to numerical limitations in the computation of the matrix coefficients. To be able to solve the problem for seamounts with larger $\epsilon$, the method would need to be adapted, in order to better handle the singularities at critical points.

Focusing on the internal tide generated by a rectilinear barotropic flow on axisymmetric Gaussian topographies, we compute the energy conversion rates obtained for different values of the topographic height and steepness. In the range of values considered for the steepness ratio $0.3 \le \epsilon \le 1.1$ and the height ratio $0.1 \le \delta \le 0.8$, we find that the WTA model overestimates the conversion rates for large topographies ($\epsilon = 0.7$ and $\delta \ge 0.5$), and more surprisingly, also for some relatively small and flat topographies ($\delta = 0.3$; $\epsilon = 0.3$). For both cases, this is due to destructive interferences of waves generated from the opposite flanks of the seamounts, in the direction of the tide. Similar observations regarding the range of validity of the WTA model in two dimensions have been obtained by Papoutsellis et al. (2023) using a coupled-mode approach. This phenomenon is thought to be relevant for a range of oceanic seamounts. In particular, prominent features with relatively gentle slopes – such as the Batavia Rise (25.67° S, 100.5° E), the Magellan Rise (7° N, 177° W) or the Ninety-East Ridge (90° E) – are potentially subject to such interference effects. The Ninety-East Ridge in the Indian Ocean has been shown to emit relatively weak internal tide energy through satellite altimetry (Zhao et al. 2016), which may be partially attributed to this mechanism. We also provide a crude estimate of the conversion rate for the Batavia Rise. By fitting an axisymmetric Gaussian to this seamount, we find $\epsilon \approx 0.4$ and $\delta \approx 0.6$. These values were obtained by taking a Brunt–Vaisälä frequency of $N^2 = 1.25 \times 10^{-6} \text{ s}^{-2}$, which corresponds to an average over the vertical extent of the seamount, based on the salinity and temperature profiles from NOAA World Ocean Atlas 2023 (Reagan et al. 2024). The $M_2$ barotropic tide is assumed to be rectilinear, with components derived from the inverse model of Egbert & Erofeeva (2002).
(version TPXO10), i.e. $U_0 = (U_0, V_0)e^{-i\omega t}$ with $U_0 = 0.21 \text{ cm s}^{-1}$, $V_0 = -0.96 \text{ cm s}^{-1}$ and $\omega = 1.4 \times 10^{-4} \text{ s}^{-1}$ for the M₂ tide. Under the WTA, the predicted conversion rate is $3 \times 10^5 \text{ W}$, while our model yields $1.8 \times 10^4 \text{ W}$, indicating a reduction of approximately 95% relative to the WTA estimate.

More surprisingly, our computations have shown for the first time the possible existence of non-radiating 3-D topographies, analogous to the 2-D configuration described by Maas (2011). While our model does not yield strictly zero radiation in the cases presented here, it would be interesting to compute the conversion rate for the 3-D equivalent of the non-radiating ridge from Maas (2011), or to investigate in more details the values of the parameter space $(\epsilon, \delta)$. In particular, an open question is whether specific 3-D topographies can lead to minimal conversion rates, as reported in two dimensions by Papoutsellis et al. (2023). If so, a natural extension would be to develop predictive criteria for these minima.

In addition, we investigate the influence of the Coriolis frequency, via the parameter $\beta$ on the internal tide generation. While rotation has little effect on the global conversion rate, we find that it significantly alters the spatial structure of the wavefield and energy flux, a result that had been previously noticed (Holloway & Merrifield 1999; Munroe & Lamb 2005; Baines 2007) but not described in detail. In particular, we observe an asymmetric distribution of the energy fluxes. In the Northern Hemisphere (where $f > 0$), the energy is radiated preferentially in a direction that is not aligned with the tidal forcing direction, but deviated counter-clockwise. We discuss how this deviation increases with both $\beta$ and the criticality parameter $\epsilon$, and find that the far-field deviation angle $\theta_\infty$ can be described by $\theta_\infty = \beta\epsilon/2$. The height ratio $\delta$ has a weaker influence, except in regimes where the topography becomes tall enough for the waves radiated from the left and right flanks to interact. The direction of the energy fluxes in these cases was not investigated here, as the radiated energy is significantly reduced. The deviation of the energy propagation observed with our model originates from taking into consideration the fully nonlinear boundary condition against the topography (recall 2.9b). This effect is not captured by the previous WTA model of Pollmann et al. (2019), which derives an asymptotic expression of the energy flux at $O(1/r)$. By expanding this asymptotic derivation at $O(1/r^2)$, we showed that the WTA does predict a small deviation, but only close enough to the topography. We believe that this effect of the Coriolis frequency on the energy flux direction could have an importance in the ocean, especially given the increasing interest of the oceanography community on the direction of internal tide energy fluxes in order to predict wave-wave interactions (Eden & Olbers 2014; Pollmann & Nycander 2023). However, the discussion provided here, as well as the empirical law for $\theta_\infty$, does not account for the influence of non-axisymmetric topographies nor the ellipticity of the barotropic tide, whose effect on the direction remains unknown. This could be the subject of a specific study.

Overall, our model offers an alternative approach to direct computational methods (Arbic et al. 2018) for solving the full 3-D problem for internal tide radiation. It opens the door to exploring more complex scenarios that go beyond the scope of this study, such as more realistic tidal forcing and bathymetry, non-radiative topographies or 3-D attractors. Even though the current work is limited to uniform stratification, it can be extended to vertically varying stratification which will allow comparisons with realistic oceanic scenarios. The case of horizontally varying stratification cannot, however, be treated using this method. Finally, upon abandoning the WTA approximation, we find a significant influence on the Coriolis frequency on both the near- and far-field wavefields and energy flux distributions. Considering the nonlinear boundary condition at the seafloor could be of similar importance for other problems dealing with 3-D inertial/internal-wave problems involving topographies, such as lee-wave generation or topographic wave scattering.
**Acknowledgements.** This work was granted access to the HPC resources of CALMIP supercomputing centre under the allocation P23045. The authors acknowledge the IMFT's scientific computing service COSINUS, and in particular R. Allis for his help with the parallelisation of the boundary element method. We also thank B. Voisin for the insightful discussions on the formulation of the integral equation. The present document is licensed under CC-BY-4.0. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/

**Funding.** This study has been partially supported through the grant EUR TESS N°ANR-18-EURE-0018 in the framework of the Programme des Investissements d'Avenir. N.G. acknowledges the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) (funding reference number RGPIN-2022-04560) and of the Canadian Space Agency (14SUSWOTTO).

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

**Data availability statement.** The data that support the findings of this study are openly available on Recherche Data Gouv at https://doi.org/10.57745/VHJRUV.

## Appendix A. Matrix coefficients

The topography is discretised using square cells $C_{ij} = [x_i - \Delta/2, x_i + \Delta/2] \times [y_j - \Delta/2, y_j + \Delta/2]$, where $\Delta$ is the cell side length. The boundary element method amounts to solving the matrix equation

$$ B = AS, \quad \text{with} \quad \left\{ \begin{array}{l} B_{kl} = \iint_{C_{kl}} U_0 \cdot \nabla_H h(x, y) \, dx \, dy, \\ \mathbf{A}_{klij} = \iint_{C_{kl}} \iint_{C_{ij}} \mathcal{F}(\mathbf{r}, h(\mathbf{r}); \mathbf{r}', h(\mathbf{r}')) \, dx \, dy \, dx' \, dy'. \end{array} \right. \quad (A1) $$

Here, the kernel function $\mathcal{F}$ involves an infinite sum on the vertical modes. In practice, we truncate this sum at a maximum mode number $M_s$. The full expression of $\mathbf{A}_{klij}$ with truncation at the mode $M_s$ is given below

$$ \mathbf{A}_{klij} = \sum_{m=1}^{M_s} \mathbf{A}_{klij}^{m} = \sum_{m=1}^{M_s} \iint_{C_{kl}} \iint_{C_{ij}} \frac{\sin(mh(\mathbf{r}'))}{2i\pi} \left( \sin(mh(\mathbf{r}))\mathcal{H}_{0}^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|) \right. \\ \left. + h_{,x}(\mathbf{r}) \frac{(x-x') + i\beta(y-y')}{\|\mathbf{r} - \mathbf{r}'\|} \cos(mh(\mathbf{r}))\mathcal{H}_{1}^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|) \right. \\ \left. + h_{,y}(\mathbf{r}) \frac{(y-y') - i\beta(x-x')}{\|\mathbf{r} - \mathbf{r}'\|} \cos(mh(\mathbf{r}))\mathcal{H}_{1}^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|) \right) \, dxdydx'dy'. \quad (A2) $$

Each coefficient involves the Hankel functions of the first type $\mathcal{H}_{0}^{(1)}$ and $\mathcal{H}_{1}^{(1)}$. These functions have an oscillatory behaviour at infinity and exhibit singularities at 0 in their imaginary part. Their expansion close to 0 is given by

$$ \mathcal{H}_{\nu}^{(1)}(z) = \begin{cases} 1 + \frac{2i}{\pi} \left( \ln\left(\frac{z}{2}\right) + \gamma \right) + o(1) & \text{if } \nu = 0, \\ -\frac{2i}{\pi z} + o(1) & \text{if } \nu = 1. \end{cases} \quad (A3) $$

In our case, there is integration of the singularity when $m\|\mathbf{r} - \mathbf{r}'\| = 0$, which occurs when the cells $C_{kl}$ and $C_{ij}$ are overlapping, i.e. $C_{kl} = C_{ij}$. This corresponds to the diagonal terms $(i, j) = (k, l)$. To deal with this singularity, as well as the oscillatory behaviour, we develop computation methods that are specific to diagonal and extra-diagonal terms.
A.1. Different cells

This case corresponds to all extra-diagonal terms of the matrix **A**. We define **r**<sub>kl</sub> = (x<sub>k</sub>, y<sub>l</sub>) and **r**<sub>ij</sub> = (x<sub>i</sub>, y<sub>j</sub>) as the centres of the cells C<sub>kl</sub> and C<sub>ij</sub>, d<sub>klij</sub> = ||**r**<sub>kl</sub> - **r**<sub>ij</sub>|| as the distance between those centres and Δ as the cell size of the regular grid.

When the mode *m* is small enough (*m*Δ ≪ 1), we can approximate the value of one mode in the sum by its value at the centres of the cells. In practice, we define a small τ such that if *m*Δ < τ, we use the values at the centres of the cells.

However, when *m* becomes large (*m*Δ > τ), the Hankel kernels $\mathcal{H}_0^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|)$ and $\mathcal{H}_1^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|)$, as well as the eigenfunction $a_m$, can become highly oscillatory. To deal with these oscillations, we adapt the numerical method from Wu & Sun (2022). In particular, we transform the Hankel function $\mathcal{H}_\nu^{(1)}$ into the following improper integral (see their (2.3) and (2.4)).

For $m > 0, -(\pi/2) < \arg x < (3\pi/2)$ and $\mathrm{Re}(\nu) > -1/2$

$$ \mathcal{H}_{\nu}^{(1)}(mx) = \frac{2 e^{-i \left(\frac{\pi}{2} \nu + \frac{\pi}{4}\right)}}{\sqrt{\pi} \Gamma(\nu + 0.5)} \left(\frac{1}{2mx}\right)^{\nu} e^{imx} I_{\nu}(mx), \quad (A4) $$

where

$$ I_{\nu}(mx) = \int_{0}^{\infty} (2mx + iy)^{\nu-0.5} y^{\nu-0.5} e^{-y} dy. \quad (A5) $$

The advantage of this method is to separate the Hankel function into an oscillatory part with a slowly varying part $I_\nu$. We then replace in the expression of the matrix coefficient for the mode *m* and we approximate the non-oscillatory parts by their value at the centres of the cells. For the oscillatory functions, we assume that the slope of the topography is constant over a cell

$$ h(\mathbf{r}) = h_{kl} + h_{kl,x}x + h_{kl,y}y, \quad (A6) $$

with $h_{kl}$, $h_{kl,x}$ and $h_{kl,y}$ are the values of the topography and its derivatives at the centre of the cell and $(x, y) = \mathbf{r} - \mathbf{r}_{kl}$. Moreover, we write

$$ \|\mathbf{r} - \mathbf{r}'\| \sim d_{klij} + \frac{(x_{kl} - x_{ij})}{d_{klij}}(x - x') + \frac{(y_{kl} - y_{ij})}{d_{klij}}(y - y'). \quad (A7) $$

This leads to

$$ \mathbf{A}_{klij}^{m} = \frac{e^{-i(md_{klij} + \frac{\pi}{4})}}{i\pi^2} \left[ I_0(md_{klij}) U_{klij}(m) U_{ijkl}(m) - i \frac{I_1(md_{klij})}{md_{klij}} V_{klij}(m) U_{ijkl}(m) \times \left( \frac{x_{kl} - x_{ij}}{d_{klij}}(h_{kl,x} - i\beta h_{kl,y}) + \frac{y_{kl} - y_{ij}}{d_{klij}}(h_{kl,y} + i\beta h_{kl,x}) \right) \right], \quad (A8) $$

with

$$ U_{klij}(m) = \iint_{C_{kl}} \sin(mh(\mathbf{r})) e^{i \frac{(x_{kl}-x_{ij})}{d_{klij}} x} e^{\frac{(y_{kl}-y_{ij})}{d_{klij}} y} dxdy, \quad (A9) $$

and

$$ V_{klij}(m) = \iint_{C_{kl}} \cos(mh(\mathbf{r})) e^{i \frac{(x_{kl}-x_{ij})}{d_{klij}} x} e^{\frac{(y_{kl}-y_{ij})}{d_{klij}} y} dxdy. \quad (A10) $$

The values of $U_{klij}(m)$ and $V_{klij}(m)$ are calculated analytically. As the integrand of $I_\nu$ is non-oscillatory and decays exponentially, we calculate $I_0(md_{klij})$ and $I_1(md_{klij})$
using a generalised Gauss–Laguerre quadrature rule with respect to the weight functions $y^{\nu-0.5}e^{-y}$.

## A.2. Same cells

This case corresponds to all diagonal terms of the matrix **A**. To deal with the singularity of the Hankel functions, we replace them by their limiting form at zero in (A2) and the expression is integrated analytically.

The diagonal terms for a given mode $m$ can then be expressed by

$$ 
\begin{aligned}
\mathbf{A}_{klkl}^m = & \frac{\sin(mh_{kl})}{2i\pi} \underbrace{\iint_{C_{kl}} \iint_{C_{kl}} \sin(mh_{kl}) \mathcal{H}_0^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|) \, dxdydx'dy'}_{I_{kl}(m)} \\
& + h_{kl,x} \underbrace{\iint_{C_{kl}} \iint_{C_{kl}} \frac{\cos(mh_{kl}) \frac{(x-x') + i\beta(y-y')}{\|\mathbf{r}-\mathbf{r}'\|} \mathcal{H}_1^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|) \, dxdydx'dy'}_{0} \\
& + h_{kl,y} \underbrace{\iint_{C_{kl}} \iint_{C_{kl}} \frac{\cos(mh_{kl}) \frac{(y-y') - i\beta(x-x')}{\|\mathbf{r}-\mathbf{r}'\|} \mathcal{H}_1^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|) \, dxdydx'dy'}_{0}.
\end{aligned}
\quad (A11) $$

The second and third integrals on the right-hand side evaluate to 0 since their integrands are odd functions of $\mathbf{r}' - \mathbf{r}$.

Moreover, using the expansion close to 0 and changing variable, we can write

$$ 
\begin{aligned}
I_{kl}(m) &= 4 \int_{x=0}^{\Delta} \int_{y=0}^{\Delta} \int_{x'=0}^{x} \int_{y'=0}^{y} \mathcal{H}_0^{(1)}(m \|\mathbf{r}'\|) \, dx'dy'dxdy \\
&= 4 \int_{x=0}^{\Delta} \int_{y=0}^{\Delta} \int_{x'=0}^{x} \int_{y'=0}^{y} \left\{ 1 + \frac{2i}{\pi} \left[ \ln\left(\frac{m}{2}\right) + \gamma \right] + \frac{i}{\pi} \ln(x'^2 + y'^2) \right\} dx'dy'dxdy \\
&= \Delta^4 \left\{ 1 + \frac{2i}{\pi} \left[ \ln\left(\frac{m\Delta}{2}\right) + \gamma - \frac{25}{12} + \frac{\ln(2)}{3} + \frac{\pi}{3} \right] \right\}.
\end{aligned}
\quad (A12) $$

This leads to

$$ \Rightarrow \mathbf{A}_{klkl}^m = \frac{\Delta^4}{2i\pi} \sin^2(mh_{kl}) \left\{ 1 + \frac{2i}{\pi} \left[ \ln\left(\frac{m\Delta}{2}\right) + \gamma - \frac{25}{12} + \frac{\ln(2)}{3} + \frac{\pi}{3} \right] \right\}. \quad (A13) $$

This expression is only true for a mode $m$ such that $m\Delta \ll 1$. In practice, we use the same parameter $\tau$ as for the non-superposed case, and use the expression in (A13) when $m\Delta < \tau$.

When the mode $m$ is larger, we divide the cell into $k$ smaller cells such that the new cell size $\Delta_{new}$ verifies $m\Delta_{new} < \tau$ and apply the same procedure as for the bigger cells. This method is expensive numerically. However, since it is only performed for the diagonal of the matrix, the impact on the full construction of the matrix is still limited.

## Appendix B. Expression of modal energy flux

Given that $p_m = im^{-1}g_m * S$, we express the modal pressure as

$$ p_m = \frac{1}{2\pi m} \iint_{\mathbb{R}} \sigma_m(\mathbf{r}') \mathcal{H}_0^{(1)}(m \|\mathbf{r} - \mathbf{r}'\|) \, dx'dy' \quad \text{with} \quad \sigma_m(\mathbf{r}) = \sin(mh(\mathbf{r}))S(\mathbf{r}). \quad (B1) $$
At a point $\mathbf{r}$ far from all sources ($\|\mathbf{r}'\| \ll \|\mathbf{r}\|$), the Hankel functions $\mathcal{H}_0^{(1)}$ and $\mathcal{H}_1^{(1)}$ can be expressed using asymptotic expansions, namely,

$$ \mathcal{H}_0^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|) = \sqrt{\frac{2}{\pi m r}} e^{i(mr-me_r \cdot \mathbf{r}' - \frac{\pi}{4})} \left(1 + \frac{e_r \cdot \mathbf{r}'}{2r} - \frac{i}{8mr} + \mathcal{O}(r^{-2})\right), \quad (B2) $$

$$ \mathcal{H}_1^{(1)}(m\|\mathbf{r}-\mathbf{r}'\|) = \sqrt{\frac{2}{\pi m r}} e^{i(mr-me_r \cdot \mathbf{r}' - \frac{3\pi}{4})} \left(1 + \frac{e_r \cdot \mathbf{r}'}{2r} + \frac{3i}{8mr} + \mathcal{O}(r^{-2})\right). \quad (B3) $$

Replacing $\mathcal{H}_0^{(1)}$ by its expansion in (B2), as was done in Pollmann et al. (2019), we obtain the approximation (3.22) for the modal pressure, a combination of integrals depending only on the mode number $m$ and the polar coordinate $\theta$.

We use the same approach to express the radial and azimuthal components of the modal pressure gradient

$$ \begin{aligned} p_{m,r} &= \cos(\theta)p_{m,x} + \sin(\theta)p_{m,y} \\ &= -\frac{1}{2\pi} \iint_{\mathcal{R}} \frac{e_r \cdot (\mathbf{r} - \mathbf{r}')}{\|\mathbf{r} - \mathbf{r}'\|} \mathcal{H}_1^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|) \sigma_m(\mathbf{r}') \, dx'dy' \\ &= -\frac{1}{2\pi} \sqrt{\frac{2}{\pi m r}} e^{i(mr - \frac{3\pi}{4})} \iint_{\mathcal{R}} \left(1 + \frac{e_r \cdot \mathbf{r}'}{2r} + \frac{3i}{8mr} + \mathcal{O}(r^{-2})\right) \sigma_m(\mathbf{r}') e^{-ime_r \cdot \mathbf{r}'} \, dx'dy' \\ &= \frac{i}{2\pi} \sqrt{\frac{2}{\pi m r}} e^{i(mr - \frac{\pi}{4})} \left(\left(1 + \frac{3i}{8mr}\right) f_m(\theta) + i \frac{g_m(\theta)}{2r} + \mathcal{O}(r^{-2})\right), \end{aligned} \quad (B4) $$

and

$$ \begin{aligned} \frac{p_{m,\theta}}{r} &= -\sin(\theta)p_{m,x} + \cos(\theta)p_{m,y} \\ &= -\frac{1}{2\pi} \iint_{\mathcal{R}} \frac{e_{\theta} \cdot (\mathbf{r} - \mathbf{r}')}{\|\mathbf{r} - \mathbf{r}'\|} \mathcal{H}_1^{(1)}(m\|\mathbf{r} - \mathbf{r}'\|) \sigma_m(\mathbf{r}') \, dx'dy' \\ &= \frac{1}{2\pi} \sqrt{\frac{2}{\pi m r}} e^{i(mr - \frac{3\pi}{4})} \iint_{\mathcal{R}} \left(\frac{e_{\theta} \cdot \mathbf{r}'}{r} + \mathcal{O}(r^{-2})\right) \sigma_m(\mathbf{r}') e^{-ime_r \cdot \mathbf{r}'} \, dx'dy' \\ &= \frac{1}{2\pi} \sqrt{\frac{2}{\pi m r}} e^{i(mr - \frac{\pi}{4})} \left(\frac{f'_m(\theta)}{mr} + \mathcal{O}(r^{-2})\right). \end{aligned} \quad (B5) $$

This leads to an $\mathcal{O}(r^{-2})$ approximation of the modal energy flux (3.15), given in (3.25).

## Appendix C. Model validation

We analyse the internal tide generated by a unidirectional tide $U_0 = e_x$ over an axisymmetric Gaussian topography, defined as

$$ h(x, y) = -\pi + \pi \delta e^{-\frac{|x|^2}{2L^2}}, \quad \text{if } |x|, |y| < x_{\text{max}}, \quad (C1) $$

where $\delta$ is the height ratio and $L$ the characteristic width of the seamount. The topography is confined within the domain $[-x_{\text{max}}, x_{\text{max}}]$, with $h$ set to $-\pi$ and $S$ to zero beyond this region. The numerical implementation relies on solving a boundary integral equation using a square mesh with cell size $\Delta$ in both the $x$ and $y$ directions. The method considers up to $M_s$ modes for the computation of matrix coefficients.

To evaluate the accuracy and convergence of the approach, we examine the energy conversion rate for the two seamounts presented in § 4, both the subcritical ($\epsilon = 0.7$) and the supercritical ($\epsilon = 1.1$) seamounts. The height ratio is kept fixed to $\delta = 0.3$ for both
<figure><img src="image_15.png" /><figcaption>Figure 13. Evolution of the energy conversion rate, computed with 20 modes, with respect to the maximum mode number $M_s$ and the mesh cell size $\Delta$ for (a) $\epsilon = 0.7$ and (b) $\epsilon = 1.1$. The inset figures show the convergence of the relative error, as defined in (C2).</figcaption></figure>

cases. Rotational effects are neglected by setting $\beta = 0$. The convergence of the boundary element method when considering rotation was also investigated and yields similar results. They are thus not presented here.

We focus first on the convergence of the method with respect to the cell size $\Delta$ and mode number $M_s$, keeping the domain size fixed at $x_{max} = 3L$. Figure 13 illustrates the evolution of the energy conversion rate, calculated by summing the contributions of the first 20 modes, with the maximum mode number $M_s$ for different mesh refinements. For the subcritical case (a), at a given mesh resolution $\Delta$, the conversion rate is initially overestimated for small $M_s$ and decreases monotonically toward an asymptotic value. This value increases with mesh resolution and converges to a finite limit. For the supercritical case (b), the conversion rate is similarly overestimated for small $M_s$, then reaches a minimum and increases again toward an asymptotic value. Improving the mesh resolution, however, leads to a decrease of the estimated conversion rate. Both cases require a high number of modes to achieve proper convergence of the source distribution ($M_s > 100$), even though these higher modes do not account for a larger amount of the energy. This is linked to the fact that modal velocities amplitude decrease more slowly with mode number than modal energy flux amplitude. A large number of modes are thus necessary to accurately describe the boundary condition at the topography. This might also be related to the fact that our formulation does not account for an hydrostatic adjustment of the barotropic tide at the seamount, as discussed in Papoutsellis et al. (2023). This adjustment, which corresponds to a mode 0, is thus redistributed in all the other modes in our case.

We define as $C_{best}$ the 'best' calculated conversion rate, i.e. the conversion rate for which the source distribution was computed with the smallest cell size $\Delta$ and largest number of mode $M_s$. The inset figures 13(a) and 13(b) displays the relative error

$$ Error = \frac{|C - C_{best}|}{C_{best}}, \qquad (C2) $$

as a function of $M_s$. For both topographies, the convergence rate with $M_s$ is similar for all mesh sizes. However, it is approximately proportional to $M_s^{-3/2}$ for the subcritical seamount and to $M_s^{-3/4}$ for the supercritical one. The convergence with $M_s$ is thus slower in the supercritical case. This partially explain the numerical difficulties observed in supercritical cases, as the solution requires more vertical modes.
BREBBIA, C.A., TELLES, J.C.F. & WROBEL, L.C. 2012 *Boundary Element Techniques: Theory and Applications in Engineering*. Springer Sci.& Bus. Media.

COUTO, R.T. 2013 Green's functions for the wave, Helmholtz and Poisson equations in a two-dimensional boundless domain. *Rev. Bras. Ens. Fís.* **35**, 1304.

ECHEVERRI, P. & PEACOCK, T. 2010 Internal tide generation by arbitrary two-dimensional topography. *J. Fluid Mech.* **659**, 247–266.

EDEN, C. & OLBERS, D. 2014 An energy compartment model for propagation, nonlinear interaction, and dissipation of internal gravity waves. *J. Phys. Oceanogr.* **44**, 2093–2106.

EGBERT, G.D. & EROFEEVA, S.Y. 2002 Efficient inverse modeling of barotropic ocean tides. *J. Atmos. Ocean. Technol.* **19**, 183–204.

EGBERT, G.D. & RAY, R.D. 2001 Estimates of $m_2$ tidal energy dissipation from TOPEX/Poseidon altimeter data. *J. Geophys. Res.* **106**, 22475–22502.

EGBERT, G.D. & RAY, R.D. 2000 Significant dissipation of tidal energy in the deep ocean inferred from satellite altimeter data. *Nature* **405**, 775–778.

FALAHAT, S., NYCANDER, J., ROQUET, F. & ZARROUG, M. 2014 Global calculation of tidal energy conversion into vertical normal modes. *J. Phys. Oceanogr.* **44**, 3225–3244.

GARRETT, C. & GERKEMA, T. 2007 On the body-force term in internal-tide generation. *J. Phys. Oceanogr.* **37**, 2172–2175.

GARRETT, C. & KUNZE, E. 2007 Internal tide generation in the deep ocean. *Annu. Rev. Fluid Mech.* **39**, 57–87.

GEOFFROY, G., POLLMANN, F. & NYCANDER, J. 2024 Tidal conversion into vertical normal modes by near-critical topography. *J. Phys. Oceanogr.* **54** (9), 1949–1970.

GREEN, J.A.M. & NYCANDER, J. 2013 A comparison of tidal conversion parameterizations for tidal models. *J. Phys. Oceanogr.* **43**, 104–119.

GRIFFITHS, S.D. & GRIMSHAW, R.H.J. 2007 Internal tide generation at the continental shelf modeled using a modal decomposition: two-dimensional results. *J. Phys. Oceanogr.* **37**, 428–451.

GAYEN, B. & SARKAR, S. 2011 Direct and large-eddy simulations of internal tide generation at a near-critical slope. *J. Fluid Mech.* **681**, 48–79.

HOLLOWAY, P.E. & MERRIFIELD, M.A. 1999 Internal tide generation by seamounts, ridges, and islands. *J. Geophys. Res.* **104**, 25937–25951.

KELLY, S.M. 2016 The vertical mode decomposition of surface and internal tides in the presence of a free surface and arbitrary topography. *J. Phys. Oceanogr.* **46**, 3777–3788.

KING, B., ZHANG, H.P. & SWINNEY, H.L. 2010 Tidal flow over three-dimensional topography generates out-of-forcing-plane harmonics. *Geophys. Res. Lett.* **37**, L14606.

DE LAVERGNE, C., FALAHAT, S., MADEC, G., ROQUET, F., NYCANDER, J. & VIC, C. 2019 Toward global maps of internal tide energy sinks. *Ocean Model.* **137**, 52–75.

LE DIZÉS, S. 2024 Critical slope singularities in rotating and stratified fluids. *Phys. Rev. Fluids* **9**, 034803.

LLEWELLYN SMITH, S. & YOUNG, W.R. 2002 Conversion of the barotropic tide. *J. Phys. Oceanogr.* **32**, 1554–1566.

LLEWELLYN SMITH, S.G. & YOUNG, W.R. 2003 Tidal conversion at a very steep ridge. *J. Fluid Mech.* **495**, 175–191.

MAAS, L.R.M. 2011 Topographies lacking tidal conversion. *J. Fluid Mech.* **684**, 5–24.

MACKINNON, J.A., et al. 2017 Climate process team on internal wave-driven ocean mixing. *Bull. Am. Meteorol. Soc.* **98**, 2429–2454.

MARTIN, P.A. & LLEWELLYN SMITH, S.G. 2012 Internal gravity waves, boundary integral equations and radiation conditions. *Wave Motion* **49**, 427–444.

MATHUR, M., CARTER, G.S. & PEACOCK, T. 2016 Internal tide generation using green function analysis: to WKB or not to WKB? *J. Phys. Oceanogr.* **46**, 2157–2168.

MELET, A., HALLBERG, R., LEGG, S. & POLZIN, K. 2013 Sensitivity of the ocean state to the vertical distribution of internal-tide-driven mixing. *J. Phys. Oceanogr.* **43**, 602–615.

MELET, A., LEGG, S. & HALLBERG, R. 2016 Climatic impacts of parameterized local and remote tidal mixing. *J. Clim.* **29**, 3473–3500.

MUNK, W. & WUNSCH, C. 1998 Abyssal recipes II: energetics of tidal and wind mixing. *Deep-Sea Res. I* **45**, 1977–2010.

MUNROE, J.R. & LAMB, K.G. 2005 Topographic amplitude dependence of internal wave generation by tidal forcing over idealized three-dimensional topography. *J. Geophys. Res.* **110**, C02001.

MUSGRAVE, R., POLLMANN, F., KELLY, S. & NIKURASHIN, M. 2022 Chapter 6 - the lifecycle of topographically-generated internal waves. In *Ocean Mixing*, (ed. M. MEREDITH & A. NAVEIRA GARABATO), pp. 117–144. Elsevier.
NYCANDER, J. 2005 Generation of internal waves in the deep ocean by tides. *J. Geophys. Res.* **110**, C10028.

PAPOUTSELLIS, C.E., MERCIER, M.J. & GRISOUARD, N. 2023 Internal tide generation from non-uniform barotropic body forcing. *J. Fluid Mech.* **964**, A20.

POLLMANN, F. & NYCANDER, J. 2023 Resolving the horizontal direction of internal tide generation: global application for the M2 tide's first mode. *J. Phys. Oceanogr.* **53**, 1251–1267.

POLLMANN, F., NYCANDER, J., EDEN, C. & OLBERS, D. 2019 Resolving the horizontal direction of internal tide generation. *J. Fluid Mech.* **864**, 381–407.

PÉTRÉLIS, F., LLEWELLYN SMITH, S. & YOUNG, W.R. 2006 Tidal conversion at a submarine ridge. *J. Phys. Oceanogr.* **36**, 1053–1071.

REAGAN, J.R., et al. 2024 World ocean atlas 2023. Dataset: NCEI Accession 0270533. https://www.ncei.noaa.gov/archive/accession/0270533.

SANDSTROM, H. & ELLIOTT, J.A. 1984 Internal tide and solitons on the Scotian shelf: a nutrient pump at work. *J. Geophys. Res.* **89** (C4), 6415–6426.

SAVAGE, A.C., WATERHOUSE, A.F. & KELLY, S.M. 2020 Internal tide nonstationarity and wave–Mesoscale interactions in the Tasman sea. *J. Phys. Oceanogr.* **50**, 2931–2951.

ST. LAURENT, L., STRINGER, S., GARRETT, C. & PERRAULT-JONCAS, D. 2003 The generation of internal tides at abrupt topography. *Deep-Sea Res. I* **50**, 987–1003.

ST. LAURENT, L.C., SIMMONS, H.L. & JAYNE, S.R. 2002 Estimating tidally driven mixing in the deep ocean. *Geophys. Res. Lett.* **29**, 2106.

VIC, C., NAVEIRA GARABATO, A.C., GREEN, J.A.M., SPINGYS, C., FORRYAN, A., ZHAO, Z. & SHARPLES, J. 2018 The lifecycle of semidiurnal internal tides over the Northern Mid-Atlantic ridge. *J. Phys. Oceanogr.* **48**, 61–80.

VOISIN, B. 2021 Boundary integrals for oscillating bodies in stratified fluids. *J. Fluid Mech.* **927**, A3.

VOISIN, B. 2024 Internal tides and the inviscid dynamics of an oscillating ellipsoid in a stratified fluid. *J. Fluid Mech.* **999**, A59.

VOISIN, B., ERMANYUK, E.V. & FLOR, J.-B. 2011 Internal wave generation by oscillation of a sphere, with application to internal tides. *J. Fluid Mech.* **666**, 308–357.

WATERHOUSE, A.F., KELLY, S.M., ZHAO, Z., MACKINNON, J.A., NASH, J.D., SIMMONS, H., BRAHZNIKOV, D., RAINVILLE, L., ALFORD, M. & PINKEL, R. 2018 Observations of the Tasman sea internal tide beam. *J. Phys. Oceanogr.* **48**, 1283–1297.

WHALEN, C.B., DE LAVERGNE, C., NAVEIRA GARABATO, A.C., KLYMAK, J.M., MACKINNON, J.A. & SHEEN, K.L. 2020 Internal wave-driven mixing: governing processes and consequences for climate. *Nat. Rev. Earth Environ.* **1**, 606–621.

WU, Q. & SUN, M. 2022 Efficient algorithms for integrals with highly oscillatory Hankel kernels. *Results Appl. Maths* **14**, 100254.

WUNSCH, C. & FERRARI, R. 2004 Vertical mixing, energy, and the general circulation of the oceans. *Annu. Rev. Fluid Mech.* **36**, 281–314.

ZHANG, L., BUIJSMAN, M.C., COMINO, E. & SWINNEY, H.L. 2017 Internal wave generation by tidal flow over periodically and randomly distributed seamounts. *J. Geophys. Res. Oceans* **122**, 5063–5074.

ZHAO, Z., ALFORD, M.H., GIRTON, J.B., RAINVILLE, L. & SIMMONS, H.L. 2016 Global observations of open-ocean mode-1 M2 internal tides. *J. Phys Oceanogr.* **46**, 1657–1684.