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 $\varLambda$ (see figure 1). For the sake of simplicity, we will consider a rectilinear tide $\boldsymbol{\tilde {U}_0} = U_0 \cos (\omega t) \boldsymbol{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 $\boldsymbol{e}_x$ and $\boldsymbol{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. (Reference Mathur, Carter and Peacock2016). Under the $f$ -plane approximation, the Coriolis frequency $f = 2\varOmega \sin (\phi )$ , with $\Omega$ the Earth’s angular velocity is assumed spatially uniform for a given latitude $\phi$ .

Figure 1. Schematic of the problem with representations of the first three vertical modes $a_1$ , $a_2$ , $a_3$ for constant stratification.

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

(2.1a,b,c) \begin{align}& \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}, \end{align}

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

where $\boldsymbol{\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

(2.2) \begin{equation} \tilde {A}(x,y,z,t) = \frac {1}{2} \big [ A(x,y,z) e^{- i \omega t} + A^*(x,y,z) e^{i \omega t} \big ] = {Re} \big [ A(x,y,z) e^{- i \omega t} \big ], \end{equation}

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

(2.3) \begin{equation} \left ( \frac {N^2 -\omega ^2}{\omega ^2 - f^2} \right ) {\nabla} ^2_H w - w_{,zz} = 0, \end{equation}

where ${\nabla} _H^2 = \partial _{xx} + \partial _{yy}$ is the horizontal Laplacian. We define

(2.4) \begin{equation} \mu ^{-1} = \sqrt {\frac {\omega ^2-f^2}{N^2-\omega ^2}}, \end{equation}

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$ , $\varLambda$ ; frequencies $\omega$ , $f$ , $N$ ; and the velocity $U_0$ . Garrett & Kunze (Reference Garrett and Kunze2007) introduce five relevant non-dimensional numbers. The first three are taken as the frequency ratios

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

and the height ratio

(2.6) \begin{equation} \delta = \frac {\varLambda }{H_0}. \end{equation}

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,

(2.7) \begin{equation} \epsilon = \mu \max \left ( \partial _x h \right )\!. \end{equation}

This parameter measures the criticality of the topography, with $\epsilon \lt 1$ (respectively ${\gt } 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,

(2.8) \begin{equation} \alpha = \frac {U_0}{\omega L}. \end{equation}

The assumption of small excursion parameter $\alpha \ll 1$ allows us to neglect the advective term $\boldsymbol U_0 \boldsymbol{\cdot }\boldsymbol{\nabla}$ in the Boussinesq equations (see e.g. Garrett & Kunze Reference Garrett and Kunze2007 and Bell Reference Bell1975a , Reference Bellb for the complete solution). Furthermore, if we assume that $u$ is smaller than, or similar to, $U_0$ , the nonlinear terms $\boldsymbol u \boldsymbol{\cdot }\boldsymbol{\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 Reference Gáyen and Sarkar2011). 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,

(2.9a,b) \begin{equation} w|_{z=0} = 0\quad \text{and}\quad w|_{z=h} = (\boldsymbol{U}_0 + \boldsymbol{u}_H) \boldsymbol{\cdot }\boldsymbol{\nabla }_H h, \end{equation}

where $\boldsymbol{u}_H = (u,v)$ is the horizontal velocity, $\boldsymbol{r} = (x,y)$ and $\boldsymbol{\nabla }_H = (\partial _x, \partial _y)$ . The horizontal boundary conditions are chosen so that there are only waves radiated outwards as $\| \boldsymbol r \| \rightarrow \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 $\boldsymbol{u}_H\boldsymbol{\cdot }\boldsymbol{\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} = \boldsymbol{U}_0 \boldsymbol{\cdot }\boldsymbol{\nabla }_H h$ . However, for steeper topographies, the horizontal velocity $\boldsymbol{u}_H$ can be locally larger than, or of the same order as, the tidal velocity $\boldsymbol{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:

(2.10) \begin{align}& [\hat {x}, \hat {y}] = \frac {\pi [x,y]}{\mu H_0},\,\,\,\hat {z} = \frac {\pi z}{H_0}\,,\quad [\hat {u}, \hat {v}] = \frac {[u,v]}{U_0},\,\,\, \hat {w} = \frac {\mu w}{U_0}, \end{align}

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

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,

(2.12a) \begin{align}&\qquad \big [ \hat {w}, \hat {b} \big ] (\hat {x},\hat {y},\hat {z}) = \sum \big[ \hat {w}_m , \hat {b}_m \big] (\hat {x},\hat {y}) a_m(\hat {z}), \end{align}

(2.12b) \begin{align}& \left [ \hat {u}, \hat {v}, \hat {p} \right ](\hat {x},\hat {y},\hat {z}) = \sum \hat {k}_m^{-1} \left [ \hat {u}_m, \hat {v}_m, \hat {p}_m \right ] (\hat {x},\hat {y}) a'_m(\hat {z}). \end{align}

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

(2.13) \begin{equation} \frac {{\rm d}^2 a_m}{{\textrm d}z^2} + m^2 a_m = 0 \,\,\,\,\text{with } \,\,a_m(0) = a_m(-\pi ) = 0, \end{equation}

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

(2.14) \begin{equation} {\int _{-\pi }^{0} \sin (n \hat {z}) \sin (m \hat {z})\, {\textrm d}\hat {z}} = \int _{-\pi }^{0} \cos (n \hat {z}) \cos (m \hat {z})\, {\textrm d}\hat {z} = \left \{ \begin{array}{cl} \frac {\pi }{2} &\text{if $m=n$,} \\[4pt] 0 &\text{otherwise.} \end{array} \right . \end{equation}

The Boussinesq equations (2.1) can be rewritten in their dimensionless form, dropping the hats for simplicity

(2.15a,b,c) \begin{align}& u - i \beta v + i \big(1 - \beta ^2 \big)p_{,x} = 0,\quad v + i \beta u + i \big(1 - \beta ^2 \big) p_{,y} = 0,\quad w - i p_{,z} = 0, \end{align}

(2.15d,e) \begin{align}&\qquad\qquad\qquad -i b + w = 0 \quad \text{and}\quad u_{,x} + v_{,y} + w_{,z} = 0. \end{align}

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

(2.16a,b) \begin{align}& 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}), \end{align}

(2.16c,d) \begin{align}&\qquad\qquad\quad p_m = i m^{-1} w_m \quad \text{and}\quad b_m = - i w_m. \end{align}

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.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. (Reference Pétrélis, Llewellyn Smith and Young2006) to compute internal tide generation over 2-D topographies. However, their method, using vertical discretisation, was restricted to symmetric seamounts. Echeverri & Peacock (Reference Echeverri and Peacock2010) 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

(3.1) \begin{equation} {\nabla} ^2_H w - w_{,zz} = S. \end{equation}

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

(3.2) \begin{equation} {\nabla} ^2_H \mathcal{G} - \mathcal{G}_{,zz} = \delta ^{(2)}(\boldsymbol{r}-\boldsymbol{r'}) \delta ^{(1)}(z-h(\boldsymbol{r'})), \end{equation}

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

(3.3) \begin{equation} \mathcal{G}(\boldsymbol{r}, z \,; \boldsymbol{r'}, h(\boldsymbol{r'})) = \sum _{m=1}^{\infty } g_m(\boldsymbol{r}, \boldsymbol{r'}) \sin (m z). \end{equation}

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 (Reference Voisin2021) (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 & Llewellyn Smith (Reference Martin and Llewellyn Smith2012) and Voisin (Reference Voisin2021) 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$ :

(3.4) \begin{equation} {\nabla} ^2_H g_m + m^2 g_m = \frac {2 \sin (m h(\boldsymbol{r'}))}{\pi } \delta ^{(2)} (\boldsymbol{r}-\boldsymbol{r'}). \end{equation}

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

(3.5) \begin{equation} g_m (\boldsymbol{r}; \boldsymbol{r'}) = \frac { \sin (m h(\boldsymbol{r'}))}{2 i \pi } \mathcal{H}_0^{(1)}(m \| \boldsymbol{r}-\boldsymbol{r'} \|), \end{equation}

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

(3.6) \begin{equation} \mathcal{G}(\boldsymbol{r}, z; \boldsymbol{r'}, h(\boldsymbol{r'})) = \sum _{m=1}^{\infty } \frac {\sin (m h(\boldsymbol{r'})) \sin (m z) }{2 i \pi } \mathcal{H}_0^{(1)}(m \| \boldsymbol{r}-\boldsymbol{r'} \|). \end{equation}

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

(3.7) \begin{align} w(x,y,z) &= \int\!\!\!\int _{\mathcal{R}} \mathcal{G}(\boldsymbol{r}, z; \boldsymbol{r'}, h(\boldsymbol{r'})) S(\boldsymbol{r'})\, {\textrm d}x' {\textrm d}y' = \mathcal{G}*S, \end{align}

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,

(3.8) \begin{equation} u = \mathcal{G}_u*S, \,\, v = \mathcal{G}_v*S,\,\,p = \mathcal{G}_p*S,\,\,b = \mathcal{G}_b*S. \end{equation}

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

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

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

(3.9) \begin{equation} \boldsymbol{U}_0\boldsymbol{\cdot }\boldsymbol{\nabla }_H h(\boldsymbol{r}) = w\left [ \boldsymbol{r},h(\boldsymbol{r}) \right ] - h_{,x}(\boldsymbol{r}) u\left [ \boldsymbol{r},h(\boldsymbol{r}) \right ] - h_{,y}(\boldsymbol{r}) v\left [ \boldsymbol{r},h(\boldsymbol{r}) \right ]\!, \end{equation}

can be rewritten as an integral equation problem

(3.10) \begin{equation} \boldsymbol{U}_0\boldsymbol{\cdot }\boldsymbol{\nabla }_H h(\boldsymbol{r}) = \mathcal{F} * S = \int\!\!\!\int _{\mathcal{R}} \mathcal{F}(\boldsymbol r, h(\boldsymbol{r});\, \boldsymbol{r'}, h(\boldsymbol{r'})) S(\boldsymbol{r'})\, {\textrm d}x' {\textrm d}y', \end{equation}

where

(3.11) \begin{equation} \mathcal{F} = \mathcal{G} - h_{,x} \mathcal{G}_u - h_{,y} \mathcal{G}_v. \end{equation}

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 (Reference Martin and Llewellyn Smith2012). While a complete analytical proof of this continuity – such as the one outlined in Brebbia et al. (Reference Brebbia, Telles and Wrobel2012, 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 $(\boldsymbol{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 Reference Le Dizès2024). 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 $\mathcal{C}_{\textit{ij}} = [x_i -\varDelta /2, x_{i}+\varDelta /2] \times [y_j-\varDelta /2, y_{j}+\varDelta /2]$ in the following, where $\varDelta$ is the cell side length (see figure 1). We then assume that $S$ is constant over a cell: $\forall \boldsymbol r \in \mathcal{C}_{\textit{ij}}, S(\boldsymbol r) = S_{\textit{ij}}$ . This allows us to write

(3.12) \begin{equation} \boldsymbol{U}_0\boldsymbol{\cdot }\boldsymbol{\nabla }_H h(\boldsymbol{r}) = \sum _{\textit{ij}} S_{\textit{ij}} \int\!\!\!\int _{\mathcal{C}_{\textit{ij}}} \mathcal{F}(\boldsymbol r, h(\boldsymbol{r});\, \boldsymbol{r'}, h(\boldsymbol{r'}))\, {\textrm d}x' {\textrm d}y'. \end{equation}

Integrating (3.12) over the cell $\mathcal{C}_{\textit{kl}}$ , we obtain the following matrix equation:

(3.13) \begin{equation} B = \unicode{x1D63C} S, \quad \text{with}\quad \left \{ \!\!\begin{array}{l} B_{\textit{kl}} = \int\!\!\int _{\mathcal{C}_{\textit{kl}}} \boldsymbol{U}_0\boldsymbol{\cdot }\boldsymbol{\nabla }_H h(x,y)\, {\textrm d}x\,{\textrm d}y, \\[8pt] \unicode{x1D63C}_{\textit{klij}} = \int\!\!\int _{\mathcal{C}_{\textit{kl}}} \int\!\!\int _{\mathcal{C}_{\textit{ij}}} \mathcal{F}(\boldsymbol r, h(\boldsymbol{r});\, \boldsymbol{r'}, h(\boldsymbol{r'}))\, {\textrm d}x\,{\textrm d}y\, {\textrm d}x'{\textrm d}y'. \end{array} \right . \end{equation}

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 $\unicode{x1D63C}$ 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

(3.14) \begin{equation} \boldsymbol{\tilde {J}} = \int _{-\pi }^{0} \langle \tilde {p} \boldsymbol{\tilde {u}}\rangle {\textrm d}z = \frac {1}{2} \int _{-\pi }^{0} {Re} \big(p^* \boldsymbol{u} \big)\, {\textrm d}z = \sum _{m=1}^\infty \boldsymbol{\tilde {J}}_m, \end{equation}

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

(3.15) \begin{equation} \boldsymbol{\tilde {J}}_m = {Re} \left [ \frac {\pi }{4} p_m^*(\boldsymbol r) \boldsymbol{u}_m(\boldsymbol r)\right ]\!. \end{equation}

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 \varOmega$ of height $\pi$ around the topography, namely,

(3.16) \begin{equation} \tilde {C} = \int\!\!\!\int _{\partial \varOmega } \langle \tilde {p} \boldsymbol{\tilde {u}}\rangle \boldsymbol{\cdot }{\boldsymbol{e}_r} {\textrm d}s = \int _0^{2\pi } \boldsymbol{\tilde {J}} \boldsymbol{\cdot }\boldsymbol{e}_r r {\textrm d}\theta = \sum _{m=1}^\infty \tilde {C}_m, \end{equation}

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

(3.17) \begin{equation} \tilde {C}_m = \int _0^{2\pi } \boldsymbol{\tilde {J}}_m \boldsymbol{\cdot }\boldsymbol{e}_r r {\textrm d}\theta \end{equation}

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 $\boldsymbol{\tilde {J}}$ and the conversion rate $\tilde {C}$ are scaled respectively by the quantities

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

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^\star$ , defined as the smallest mode such that the contribution $C_{M_f^\star }$ represents less than $1\,\%$ of the cumulative energy of all the lower modes, i.e. $C_{M_f^\star } \leqslant 0.01 \sum _{m=1}^{M_f^\star } C_m$ . The values of $M_f^\star$ , 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^\star$ 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.

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 \}$ .

Finally, as done by Pollmann et al. (Reference Pollmann, Nycander, Eden and Olbers2019), we are also interested in an asymptotic expression of the modal energy flux. This asymptotic expression is used in § 6 to compare the results obtained by our method and the WTA. First, we write the modal horizontal velocity in the polar coordinates $(r,\theta )$

(3.19) \begin{equation} \boldsymbol{u}_m = - i \left (\!\! \begin{array}{c} p_{m,x} + i \beta p_{m,y} \\[5pt] p_{m,y} - i \beta p_{m,x} \end{array} \!\!\right )_{(\boldsymbol{e}_x,\boldsymbol{e}_y)} = - i \left ( \!\!\begin{array}{c} p_{m,r} + i \beta p_{m,\theta }/r \\[5pt] p_{m,\theta }/r - i \beta p_{m,r} \end{array} \!\!\right )_{(\boldsymbol{e}_r,\boldsymbol{e}_{\theta })}\!. \end{equation}

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

(3.20) \begin{align} p_m &= \frac {1}{2 \pi m} \int\!\!\!\int _{\mathcal{R}}\sigma _m(\boldsymbol{r'}) \mathcal{H}_0^{(1)}(m \| \boldsymbol{r}- \boldsymbol{r'}\|)\, {\textrm d}x' {\textrm d}y', \end{align}

with $\sigma _m$ a modal forcing term defined as

(3.21) \begin{equation} \sigma _m(\boldsymbol{r})= \sin (m h(\boldsymbol{r})) S(\boldsymbol{r}). \end{equation}

Using the asymptotic expansion of the Hankel functions, similar to what is done by Pollmann et al. (Reference Pollmann, Nycander, Eden and Olbers2019) (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

(3.22) \begin{equation} p_m(\boldsymbol{r}) \sim \frac {1}{2 \pi m} \sqrt {\frac {2}{\pi m r}} {e^{i (m r - \pi /4)}} \left [ \left ( 1-\frac {i}{8 m r}\right ) f_m(\theta ) + i \frac {g_m(\theta )}{2r} \right ]\!, \end{equation}

where

(3.23) \begin{align} f_m(\theta ) &\equiv \int\!\!\!\int \sigma _m(\boldsymbol{r'}) e^{- i m \boldsymbol{e}_r\boldsymbol{\cdot }\boldsymbol{r'} } {\textrm d}x' {\textrm d}y' = \mathfrak{F}[\sigma _m](m \boldsymbol{e}_r) { = \mathfrak{F}[\sigma _m](m\cos \theta ,m\sin \theta )}, \end{align}

(3.24) \begin{align} g_m(\theta ) &\equiv \int\!\!\!\int (-i \boldsymbol{e}_r)\boldsymbol{\cdot }\boldsymbol{r'} \sigma _m(\boldsymbol{r'}) e^{- i m \boldsymbol{e}_r\boldsymbol{\cdot }\boldsymbol{r'} } {\textrm d}x' {\textrm d}y' = (\partial _{k} \mathfrak{F}[\sigma _m])(m \boldsymbol{e}_r), \end{align}

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 $\boldsymbol{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 $\mathcal{O}(r^{-2})$ of the modal energy flux

(3.25) \begin{align} \boldsymbol{\tilde {J}}_m &= \frac {1}{8 \pi ^2 m^2 r} \left \{ \left (|f_m|^2(\theta ) + \beta \frac {{\rm{Re}} \left [ f_m^*(\theta ) f_m^{\prime}(\theta )\right ]}{m r} + \frac {\textrm {Im} \left [ f_m(\theta ) g_m^*(\theta )\right ]}{r} \right ) \boldsymbol{e}_r \right . \nonumber \\ &\quad + \left . \left ( \beta \frac {|f_m|^2(\theta )}{2m r} + \frac {\textrm {Im}\left [f_m^*(\theta ) f_m^{\prime}(\theta )\right ]}{m r} \right ) \boldsymbol{e}_{\theta } \right \}\!. \end{align}

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.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 $\boldsymbol{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:

(4.1) \begin{equation} {h(x,y) = \left \{\!\! \begin{array}{cl} - \pi + \pi \delta e^{ - \frac {\| \boldsymbol r \|^2}{2 L^2}} & \text{ if $|x|,|y| \leqslant 3L$},\\[6pt] -\pi & \text{ otherwise.} \end{array} \right .} \end{equation}

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.

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

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.

The topography spans a mesh of $250 \times 250$ points, resulting in a horizontal mesh size of $\varDelta \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 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\|\lt 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 \lt 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 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.

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 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 \,{\textrm m}^{-6}$ in black and $500 \,{\textrm m}^{-3}$ in grey ; (b) cumulative sum up to mode $M_f$ scaled by model predictions: $C_{\textit{model}} = C^{\textit{WTA}}$ (subcritical) and $C_{\textit{model}} = C^{\textit{Voisin}}$ (supercritical).

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. (Reference Papoutsellis, Mercier and Grisouard2023) (see their figure 3 a). 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_{\textit{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 (Reference Llewellyn Smith and Young2002), computed over the first 20 modes

(4.2) \begin{equation} C^{\textit{WTA}} = \frac {\pi ^3}{2} \delta ^2 L^4 \sum _{m=1}^{20} C_m^{\textit{WTA}}, \,\,\,\,\,\, \text{with}\,\,\, C_m^{\textit{WTA}} = m^2 \exp \big ({-}m^2 \hat {L}^2 \big ), \end{equation}

where $C^{\textit{WTA}}$ and $C_m^{\textit{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^{\textit{WTA}} \approx 0.505$ . For the supercritical seamount, we compare our estimate with the value from Voisin (Reference Voisin2024) for an equivalent spheroid in an infinite ocean. As discussed by this author, the equivalent spheroid has the same height $\varGamma$ 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^{\prime} = \mu \varLambda / d \approx 2.143$ , which yields $C^{\textit{Voisin}} \approx 0.104$ . In both cases, the conversion increases with $M_f$ and reaches a plateau at a certain threshold: $M_f^\star = 5$ (subcritical) and $M_f^\star = 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 width $L$ of the Gaussian is adapted according to $L = \mu H_0 \delta /(\epsilon \sqrt {e})$ . All tested cases are summarised in table 2. The grid cell size $\varDelta$ 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 $\varDelta = 0.0697$ , which yields a maximum relative error on the conversion rate of approximately $5\,\%$ .

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.

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. (Reference Papoutsellis, Mercier and Grisouard2023) (see their figure 6), energy conversion rates are in good agreement with the WTA when $0.5\leqslant \epsilon \leqslant 0.9$ for a Gaussian topography with $\delta =0.3$ . When the topography becomes supercritical $\epsilon \geqslant 1$ , both our model and theirs deviate smoothly from the WTA. However, a stronger and more abrupt deviation is observed for $\epsilon \leqslant 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. (Reference Papoutsellis, Mercier and Grisouard2023), we observe a strong deviation from the WTA estimate for $\delta \geqslant 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. Reference Pétrélis, Llewellyn Smith and Young2006; Maas Reference Maas2011; Papoutsellis et al. Reference Papoutsellis, Mercier and Grisouard2023; Geoffroy et al. Reference Geoffroy, Pollmann and Nycander2024), 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 with $\beta$ for larger $\epsilon$ and smaller $\delta$ . This non-trivial effect of $\beta$ on the conversion rate was observed by Baines (Reference Baines2007, see their figure 7) for a cylindrical pillbox of height $h_m$ and radius $a$ , in the limit of high $a \omega / N H_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.

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 $\boldsymbol{\tilde {J}}$ reconstructed with $M_f^\star =2$ modes.

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 (Reference Munroe and Lamb2005, 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 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.

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 sightly deviated.

Figure 9. Energy density flux generated by an 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^\star =5$ . The dashed lines represent the locations of the maximum radial flux for a given radius.

To quantify this behaviour, we introduce the deviation angle $\theta _{\textit{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 _{\textit{max}}(r))$ , tracked for each lobe, are represented by the dashed lines in figure 9. In the following, we only focus on the right side of the topography ( $x\gt 0$ ), for which $\theta _{\textit{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 (Reference Llewellyn Smith and Young2002) and later Pollmann et al. (Reference Pollmann, Nycander, Eden and Olbers2019) derived analytical solutions for internal tide generation in the hydrostatic limit under the WTA (see equation 2.19 of Pollmann et al. Reference Pollmann, Nycander, Eden and Olbers2019). 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

(5.1) \begin{equation} \sigma _m^{\textit{WTA}} = (-1)^m m\, \boldsymbol{U}_0\boldsymbol{\cdot }\boldsymbol{\nabla }_H h(\boldsymbol{r}). \end{equation}

Figure 10 compares their mode-1 forcing terms $\sigma _1^{\textit{WTA}}$ with $\sigma _1$ from our model, obtained for a unidirectional tide along $\boldsymbol{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 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$ .

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.

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

(5.2) \begin{align} f_m^{\textit{WTA}}(\theta ) = \varsigma _m \cos (\theta )\quad \text{and}\quad g_m^{\textit{WTA}}(\theta ) = \xi _m \cos (\theta ), \end{align}

with $\varsigma _m = -i (-1)^m m^2\,\mathfrak{F}[h](m \boldsymbol{e}_r)$ and $\xi _m = -i (-1)^m m\, (1-m^2 L^2 )\mathfrak{F}[h](m \boldsymbol{e}_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 $6 L$ (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

(5.3) \begin{equation} \boldsymbol{\tilde {J}_m^{\textit{WTA}}} = \frac {|\varsigma _m|^2}{8 \pi ^2 m^2 r} \cos ^2(\theta ) \left \{ \left (1 - \frac {\beta \tan (\theta )}{m r} \right ) \boldsymbol{e}_r + \frac {\beta }{2m r} \boldsymbol{e}_{\theta } \right \}\!. \end{equation}

An interesting observation is that, at $\mathcal 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

(5.4) \begin{equation} \theta _{\textit{max}}^{\textit{WTA}}(r) \sim - \frac {\beta }{2 \chi r},\quad \text{with}\quad \chi = \frac {\sum m^{-2} |\varsigma _m|^2}{\sum m^{-3} |\varsigma _m|^2}. \end{equation}

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. (Reference Pollmann, Nycander, Eden and Olbers2019), who considered only the leading-order $\mathcal{O}(1/r)$ terms for which the flux remains radial and symmetric about the $x$ -axis.

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 _{\textit{max}}$ at which the radial flux is maximum, evolve with the distance $r$ and the Coriolis parameter $\beta$ .

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_f^\star =5$ . (b) Direction of the maximum radial flux $\theta _{\textit{max}}$ as a function of the inverse radii $1/r$ (symbols ‘+’), with the least-squares linear fit materialised by the dashed lines. The symbols ‘ $\bullet$ ’ are the deviations $\theta _\infty$ given by the asymptotic expression.

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 _{\textit{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 _{\textit{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

(6.1) \begin{equation} (r\tilde {J}_r)_\infty = \frac {1}{8 \pi ^2} \left ( \sum _m \frac {|f_m|^2(\theta )}{m^2} \right )\!. \end{equation}

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_f^\star$ 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 $60 L$ , before computing the Fourier transform $\mathfrak{F}[\sigma _m](m \boldsymbol{e}_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 ‘ $\bullet$ ’ 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 (Reference Baines2007) (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

(6.2) \begin{equation} {\theta _{\infty }} \sim \frac {\beta \,\epsilon }{2}. \end{equation}

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.

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

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.