# Internal tide generation from non-uniform barotropic body forcing

Christos E. Papoutsellis1  Matthieu J. Mercier1  Nicolas Grisouard2
<EMAIL_ADDRESS>1Institut de Mécanique des Fluides de Toulouse
(IMFT), Université de Toulouse, CNRS, Toulouse, France 2University of Toronto,
Department of Physics, Toronto, Canada

###### Abstract

We model linear, inviscid, internal tides generated by the interaction of a
barotropic tide with one-dimensional topography. Starting from the body-
forcing formulation of the hydrodynamic problem, we derive a Coupled-Mode
System (CMS) using a local eigenfunction expansion of the stream function. For
infinitesimal topography, we solve this CMS analytically, recovering the
classical Weak Topography Approximation (WTA) formula for the barotropic-to-
baroclinic energy conversion rate. For arbitrary topographies, we solve this
CMS numerically. The CMS enjoys faster convergence with respect to existing
modal solutions and can be applied in the subcritical and supercritical
regimes for both ridges and shelf profiles. We show that the non-uniform
barotropic tide affects the baroclinic field locally over topographies with
large slopes and we study the dependence of the radiated energy conversion
rate on the criticality. We show that non-radiating or weakly radiating
topographies are common in the subcritical regime. We also assess the region
of validity of the WTA approximation for the commonly used Gaussian ridge and
a compactly supported bump ridge studied here for the first time. Finally, we
provide numerical evidence showing that in the strongly supercritical regime
the energy conversion rate for a ridge (resp. shelf) approaches the value
obtained by the knife-edge (resp. step) topography.

##  1 Introduction

Internal tides (ITs) are internal waves generated in the interior of a
stratified and rotating ocean through the interaction of the astronomically
induced barotropic tidal flow with the seafloor. They are oscillatory
perturbations of the baroclinic flow at the tidal frequency, propagating away
from bottom irregularities such as ridges or continental shelves (Garrett &
Kunze, 2007). ITs are considered to be one of the main sinks of energy for the
barotropic tide, and an important contributor to ocean mixing on a global
scale (Garrett, 2003; Wunsch & Ferrari, 2004; Whalen et al., 2020). An
accurate description of the IT flow and the associated barotropic-to-
baroclinic energy conversion is necessary for reliable ocean and climate
modelling. Climate models often implement internal-wave-driven mixing through
parameterizations involving idealised estimations of the local energy
conversion rate at the IT generation sites. Moreover, the dissipation of ITs
is also parameterized by using a modal description of the flow (Klymak et al.,
2013; MacKinnon et al., 2017).

The problem of IT generation, even in its linearised time-harmonic version
considered herein, is quite complex. Simplifications are commonly based on
smallness assumptions on the relative topographic height or the so-called
criticality defined as the ratio of maximum topographic slope and the
characteristic slope of internal waves at the tidal frequency (Garrett &
Kunze, 2007).

Bell (1975) linearised the bottom boundary condition around a flat bottom and
introduced a uniform barotropic background flow, that is, a barotropic flow
that does not depend on the variable topography. This approach, also known as
the “Weak Topography Approximation" (WTA) is formally valid for topographic
features of small relative height and criticality (Llewellyn Smith & Young,
2002; Khatiwala, 2003; Vlasenko et al., 2005). Another idealised configuration
of opposing nature is to consider discontinuous (infinitely steep)
topographies such as one with zero elevation everywhere except at a single
point (“knife edge") (Larsen, 1969; Llewellyn Smith & Young, 2003; St. Laurent
et al., 2003; Nycander, 2006), or top-hat, linear ramp and step functions
(Prinsenberg et al., 1974; New, 1988; St. Laurent et al., 2003). An attractive
feature of the above simplifications is that they can lead to computationally
inexpensive methods of calculating the radiating energy, which can be applied
on a global scale, see e.g.,  Nycander (2005); Falahat et al. (2014).

The first solution to the full 2D linear problem is due to Baines (1973), who
proposed what is now called the body-forcing formulation. In this approach,
the baroclinic flow appears as a response to a non-uniform barotropic flow.
That is, a flow with horizontal and vertical velocities that depend on the
variable topography and its slope. The spatial part of the response stream
function is governed by a hyperbolic partial differential equation (PDE) with
variable forcing and homogeneous boundary conditions at the natural
boundaries. Baines (1973) proposed a numerical technique based on an integral
equation derived from the normal form of this PDE and calculated the radiated
energy for simple topographies. He noted, however, that this technique becomes
rather involved for complex topographies and proposed in Baines (1982) a
perturbative method for small-scale topographies. Gerkema et al. (2004)
obtained numerical solutions of Baines (1973)’s formulation in the time domain
by treating the radiation conditions with sponge layers. Garrett & Gerkema
(2007) showed that the body-forcing term in Baines (1973)’s formulation is
inconsistent with non-hydrostatic conditions and derived a consistent
formulation. We also adopt this formulation.

An equivalent formulation is also possible where the governing hyperbolic PDE
on the total flow is homogeneous and the forcing appears as a non-homogeneous
Dirichlet condition on the bottom boundary (boundary forcing approach)
(Sandstrom, 1975; Stashchuk & Cherkesov, 1991; Vlasenko et al., 2005); see
also Garrett & Gerkema (2007) for the relation with the body-forcing
formulation. The numerical techniques developed for this formulation are also
based on the normal form of the PDE and its characteristics. They are rather
technical, and the treatment of supercritical topographies requires additional
attention. Moreover, in order to obtain the flow fields in the physical
domain, an additional transformation is required.

Semi-analytical methods have also been developed by using the Green’s function
of free internal waves corresponding to a radiating source in a flat-strip
domain (Robinson, 1969; Pétrélis et al., 2006; Echeverri & Peacock, 2010;
Mathur et al., 2016). By construction, this approach excludes shelf and trench
geometries and has been applied only to ridges. Moreover, it is characterised
by numerical singularities associated with the solution of an integral
equation, as we develop in Section 7. Balmforth & Peacock (2009) developed a
variant of this approach in the infinite-depth case with lateral periodic
conditions.

In this work, we develop a new semi-analytical IT model based on Garrett &
Gerkema (2007)’s body-forcing formulation. The analytical step is the exact
reformulation of the hydrodynamic problem as an infinite Coupled-Mode System
(CMS) of equations accomplished by means of an exact, local eigenfunction
expansion of the stream function. This approach, also called coupled-mode
theory, has been applied to various non-uniform waveguide problems in
acoustics (Brekhovskikh & Godin, 1992; Desaubies & Dysthe, 1995; Maurel et
al., 2014; Ivansson, 2021), elasticity (Maupin, 1988; Pagneux & Maurel, 2006;
He et al., 2019) and water waves (Porter & Staziker, 1995; Athanassoulis &
Belibassakis, 1999; Papoutsellis et al., 2018) among other disciplines. In the
context of ITs, Griffiths & Grimshaw (2007) derived a CMS from Euler’s
equations using a local vertical mode decomposition and calculated 2D ITs over
a shelf topography. Similar systems are derived in Kelly (2016) and Lahaye &
Llewellyn Smith (2020). The principal difference with the present approach is
that we work with the stream function. This has the advantage that every term
in our local modal expansion satisfies exactly the bottom boundary condition
for arbitrary topography, and the solution enjoys faster convergence. We use
this CMS to calculate the flow and the barotropic-to-baroclinic energy
conversion rates for two types of ridges for a wide and finely resolved range
of maximum slopes and heights. To our knowledge, the present work is the first
attempt to perform such an extensive set of calculations in the context of the
body-forcing formulation.

The paper is organised as follows. In Section 2, we present the body-forcing
formulation and derive its energy balance equation. In Section 3, we introduce
the modal decomposition of the stream function and derive our CMS. In Section
4, we present numerical convergence and accuracy results, and in Section 5 we
visualise our solutions. In Section 6, we consider the conversion rates for
various topographies. In Section 7 we discuss the differences of the proposed
CMS with the Green’s function method and, finally, in Section 8, we present
our conclusions.

##  2 Governing Equations

###  2.1 Posing the problem

In a two-dimensional Cartesian coordinate system $Oxz$, with the vertical axis
$z$ pointing upward, we consider a horizontally infinite layer of a density-
stratified fluid bounded from above by the ocean surface, modelled as a “rigid
lid” $\\{z=0\\}$, and from below by the impermeable bottom $\\{z=-h(x)\\}$
with $h>0$. We assume that the topography is asymptotically flat, that is, its
slope vanishes at infinity, $\lim_{x\to\pm\infty}[h_{x}]=0$, and we define
$\lim_{x\to\pm\infty}h=h_{\pm}$. The latter requirement allows us to take into
account fluid domains with different depths at infinity. The parameters
associated with the topography are the characteristic depth $h_{0}$, the
characteristic height $\Lambda$ and the characteristic horizontal scaling
length $L$. We do not limit ourselves to small heights; i.e., $\Lambda$ does
not have to be much smaller than $h_{0}$ or any other characteristic vertical
length scale.

In static equilibrium, the fluid velocity is zero and the background density
profile $\rho_{\text{eq}}(z)$ weakly departs from a constant reference density
$\rho_{0}$ so that the Boussinesq approximation applies. We further assume
that $\rho_{\text{eq}}(z)$ decreases linearly with $z$ such that the Brunt-
Väisälä frequency
$N=\sqrt{(-g/{\color[rgb]{0,0,0}\rho_{0}})\text{d}{\color[rgb]{0,0,0}\rho_{\text{eq}}}/\text{d}z}$
is constant. The hydrostatic pressure $p_{0}(z)$ is then defined by its
vertical gradient
$p_{0,z}\equiv\text{d}p_{0}/\text{d}z=-{\color[rgb]{0,0,0}\rho_{0}}g$. The
hydrodynamic problem is posed on the $f$-plane, with $f$ the Coriolis
parameter. Our aim is to find perturbations of this state driven by the
interaction of a background barotropic tidal flow with the bottom topography.
The barotropic flow oscillates with an angular frequency $\omega\in(f,N)$ for
$N>f$ and is associated with a constant volume flux $Q$ corresponding to a
uniform current $Q/h_{\pm}\cos(\omega t)$ at $x\rightarrow\pm\infty$.

Figure 1: Sketch of the set-up and summary of several parameters, used in this
article. We use a supercritical topography ($\varepsilon>1$) for purposes of
illustration, although our model can also be applied to subcritical ones.

Under the Boussinesq approximation, our setup is characterised by eight
dimensional parameters ($h_{\pm}$, $L$, $\Lambda$, $f$, $N$, $\omega$, and
$Q$), which we summarise in figure 1, measured in combinations of meters and
seconds. In the case of an isolated ridge, $h_{0}$ is the depth as
$x\rightarrow\pm\infty$, and $\Lambda$ is the maximum height of the ridge,
$\Lambda=\max\\{h_{0}-h\\}$. For a shelf, we assume that $h_{-}>h_{+}$
($h_{-}$ is the oceanward depth) so that the maximum height is
$\Lambda={\color[rgb]{0,0,0}\max\\{h_{0}-h\\}}$, where we have also chosen
$h_{0}=h_{-}$ as the characteristic depth. A complete dynamical description of
our system therefore requires five non-dimensional numbers. First, we
introduce a “funnelling ratio” that measures the reduction in cross-section of
the flow,

$\delta=\frac{\max\\{{\color[rgb]{0,0,0}h_{0}-h}\\}}{h_{0}}\sim\frac{\Lambda}{h_{0}}.$
(1)

The second and third parameters are the non-dimensional frequency $\omega/f$
and the characteristic slope of the internal wave,

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

$\theta$ being the angle of the free internal wave group velocity with respect
to the horizontal plane (see Eq. 8). Note that internal waves are hydrostatic
to a good approximation when $\omega\ll N$, or equivalently, $\mu\gg 1$. Our
fourth and fifth parameters are the relative steepness $\varepsilon$,

$\varepsilon=\mu\max\\{|\partial_{x}h|\\}\sim\frac{\mu\Lambda}{L}\,,$ (3)

and the tidal excursion $\tau$ defined by

$\tau=\frac{Q}{(h_{0}-\Lambda)\omega L}.$ (4)

The parameter $\varepsilon$ represents the ratio $\tan\alpha/\tan\theta$ where
$\alpha$ is the maximum inclination of the topography with respect to the
horizontal plane. Its purpose is to measure the criticality of the topography,
with $\varepsilon<1$ ($>1$) corresponding to the subcritical (supercritical)
regime. The parameter $\tau$ compares the typical displacement amplitude of a
water parcel above topography, $Q/[(h_{0}-\Lambda)\omega]$, with the
horizontal scale $L$. If $\tau$ is finite, the curvature of the particle
trajectories at the bottom generates internal waves with frequencies other
than $\omega$ (Bell, 1975). To ensure monochromatic disturbances, we therefore
assume $\tau\ll 1$. In the ocean, for the lunar semi-diurnal tide $M_{2}$, the
tidal excursion over flat bottom $Q/(\omega h_{0})$ is $O(100\text{ m})$
(Bell, 1975), and therefore even a moderate topographic width $L=O(10\text{
km})$ would satisfy this so-called “acoustic limit”.

Of the five parameters defined above, $\omega/f=O(1)$ corresponds to a tidal
component at mid-latitudes, and we will keep $\mu^{-1}$ fixed to a small value
(for illustration purposes, $f=10^{-4}$ s-1 is the value around latitude 45∘N,
$\omega/f=1.4$ for the M2 component and we will use $N=1.5\times 10^{-3}$ s-1,
which implies $\mu\approx 15$). Starting in Section 4, $\varepsilon$ and
$\delta$ are the parameters that we vary primarily. Finally, we adopt a
standard fixed value for $Q=120$ m2 s-1 corresponding e.g. to a barotropic
velocity amplitude at $x\to-\infty$ of $U_{0}=Q/h_{0}=4$ cm s-1 and a depth
$h_{0}=3$ km. As $\varepsilon$ and $\delta$ vary, the value of $\tau$
therefore varies between calculations while remaining small.

With the parameters defined above, we can discuss the linearity of our
equations. A first way to define it is to estimate the susceptibility of
radiated internal waves to undergo instabilities, which, in the $\tau\ll 1$,
$\varepsilon\ll 1$ regime, is small when $\varepsilon\tau\ll 1$ (e.g.,
Balmforth et al., 2002; Garrett & Kunze, 2007). In other words, the $\tau\ll
1$, $\varepsilon\ll 1$ regime is linear by construction. However, regardless
of the value of $\tau$, linearity breaks down as $\varepsilon$ increases,
implying that this parameter is a better measure of linearity (Garrett &
Kunze, 2007; Bühler & Muller, 2007; Grisouard & Bühler, 2012). In this
article, we adopt the common approach of letting $\varepsilon$ be
significantly supercritical in some cases while always solving a linear set of
equations (Pétrélis et al., 2006; Griffiths & Grimshaw, 2007; Balmforth &
Peacock, 2009; Echeverri & Peacock, 2010; Mathur et al., 2016). Indeed, we are
primarily interested in predicting conversion from a given large-scale
barotropic forcing to a topography-scale response. The subsequent non-linear
evolution of this response, which could take the form of different
instabilities (see Dauxois et al., 2018, for a review), could be addressed by
separate parameterised procedures such as that of Muller & Bühler (2009), but
would be beyond the scope of this article.

Based on the discussion above, we neglect nonlinear effects and diffusion of
momentum and buoyancy, and we model the flow by the linearised, inviscid
Boussinesq equations. Introducing the buoyancy
$b=-g(\rho/{\color[rgb]{0,0,0}\rho_{0}}-1)$, where $\rho$ is the total density
field, the governing equations are

$\displaystyle u_{t}-fv$ $\displaystyle=-p_{x},\quad v_{t}+fu=0,$ (5a)
$\displaystyle w_{t}$ $\displaystyle=-p_{z}+b,$ (5b) $\displaystyle
b_{t}+N^{2}w$ $\displaystyle=0,$ (5c) $\displaystyle u_{x}+w_{z}$
$\displaystyle=0,$ (5d)

where $(u,v,w)$ are the velocity components, $p$ is the associated pressure
perturbation divided by $\rho_{0}$, and subscripts denote partial derivatives.
On the boundaries $z=-h$ and $z=0$, we have the impermeability conditions

$h_{x}u(x,-h)+w(x,-h)=0\quad\text{and}\quad w(x,0)=0,$ $None$

and we also assume that the total flux through a vertical cross-section is
oscillating with frequency $\omega$ and constant amplitude $Q$,

$\displaystyle\int_{-h}^{0}u\>\mathrm{d}z=Q\cos\omega t.$ (7)

Introducing the time harmonic stream function
$\psi=\text{Re}\\{\phi(x,z)\text{e}^{-i\omega t}\\}$, with $u=-\psi_{z}$,
$w=\psi_{x}$, where $\phi$ is a complex amplitude and Re stands for the real
part, Eqs. (5)–(7) lead to a hyperbolic boundary value problem (BVP) on
$\phi$,

$\displaystyle\mathcal{L}_{\mu}\phi:=\left(\partial_{x}^{2}-\mu^{-2}\partial_{z}^{2}\right)\phi$
$\displaystyle=0,\quad\phi(x,0)=0\quad\text{and}\quad\phi(x,-h)=Q,$ (8)

where $\mu$ is defined in Eq. (2) (Vlasenko et al., 2005; Garrett & Gerkema,
2007). Eq. (8), also known as the boundary forcing formulation, is completed
by the requirement that $\phi$ has the form of a barotropic flow plus a flow
radiating away from the topography in the form of internal waves. We introduce
the barotropic flow in the following subsection.

###  2.2 Barotropic flow

The barotropic velocity components $U$, $V$, $W$ and scaled pressure $P$
satisfy Eqs. (5)-($None$) where the buoyancy in the vertical momentum equation
(5b) is absent, and Eq. (5c) acts as a diagnostic equation for the induced
buoyancy $B$, see Eqs. (44) in appendix A. Introducing a barotropic stream
function $\Psi=\text{Re}\\{\Phi\text{e}^{-\text{i}\omega t}\\}$ for some time
independent function $\Phi$ with $U=-\Psi_{z}$ and $W=\Psi_{x}$, we obtain
from (44),

$\displaystyle\Phi_{xx}+\mu_{0}^{-2}\Phi_{zz}=0,\quad\Phi(x,0)=0\quad\text{and}\quad\Phi(x,-h)=Q,$
(9)

with $\mu_{0}^{-2}=1-f^{2}/\omega^{2}$ (Garrett & Gerkema, 2007) where the
boundary conditions ensure that the total mass transport through any vertical
cross-section is solely due to the barotropic flow. Note that the PDE in (9)
is elliptic and is also obtained by (8) with $N\equiv 0$. The solution of (9)
can be written as

$\displaystyle\Phi=\Phi^{(0)}+\Phi^{\text{r}},\quad\text{with}\quad\Phi^{(0)}=-Qz/h,$
(10)

where $\Phi^{(0)}$ represents the hydrostatic part of the barotropic flow. It
is obtained by neglecting $W_{t}$ in (44b) and, consequently, $\Phi_{xx}$ in
(9); see Appendix A for a detailed derivation. $\Phi^{\text{r}}$ represents
the residual non-hydrostatic part which solves

$\displaystyle\Phi^{\text{r}}_{xx}+\mu_{0}^{-2}\Phi^{\text{r}}_{zz}=-\Phi^{(0)}_{xx},\quad\Phi^{\text{r}}(x,0)=0\quad\text{and}\quad\Phi^{\text{r}}(x,-h)=0.$
(11)

Note that $\Phi$ is spatially non-uniform, that is, it depends on $z$ and
$h(x)$. Also note that $\Phi^{\text{r}}$ vanishes at $x\rightarrow\pm\infty$.
We obtain a general semi-analytical solution of (11) that is valid for
arbitrary smooth $h$ by means of a modal decomposition (subsection 3.2). In
Appendix A, we derive a perturbative solution that is valid if $h_{0}/L\ll 1$.

In fact, $\Phi^{(0)}$ coincides with the leading order part of this solution.
The corresponding velocities are given by
$[U^{(0)},W^{(0)}]=[Q/h,Qz(1/h)_{x}]\cos(\omega t)$. As
$x\rightarrow\pm\infty$, $W^{(0)}=0$ and $U^{(0)}=Q/h_{\pm}\cos(\omega t)$
coincide with the spatially uniform (depth-averaged) barotropic currents far
from the topography. $U^{(0)}$ is also used as a forcing in the IT model of
Griffiths & Grimshaw (2007) in the case of a shelf without coastline.

###  2.3 The internal tide generation problem

Introducing $\phi^{\\#}=\phi-\Phi$, we obtain from (8) the BVP

$\displaystyle\mathcal{L}_{\mu}\phi^{\\#}=-\mathcal{L}_{\mu}\Phi,\quad\phi^{\\#}(x,0)=0\quad\text{and}\quad\phi^{\\#}(x,-h)=0,$
(12)

which shows that the barotropic flow $\Phi$ forces a purely baroclinic
response $\phi^{\\#}$. Alternatively, introducing
$\phi^{\dagger}=\phi^{\\#}+\Phi^{\text{r}}$ and exploiting the linearity of
$\mathcal{L}_{\mu}$, the BVP (12) becomes

$\displaystyle\mathcal{L}_{\mu}\phi^{\dagger}=-\mathcal{L}_{\mu}\Phi^{(0)}=-\Phi^{(0)}_{xx},\quad\phi^{\dagger}(x,0)=0\quad\text{and}\quad\phi^{\dagger}(x,-h)=0,$
(13)

which shows that the hydrostatic part of barotropic flow, $\Phi^{(0)}$, forces
a baroclinic response plus a non-hydrostatic barotropic one (Garrett &
Gerkema, 2007). This formulation is referred to as the body-forcing
formulation since the forcing appears only in the wave equation (Garrett &
Gerkema, 2007; Garrett & Kunze, 2007). An advantage of working with such a
formulation is that the unknown field satisfies homogeneous Dirichlet
conditions that make the application of a coupled-mode approach
straightforward (Section 3). Here, we shall proceed with (13) mainly because
$\Phi^{(0)}$ is given by the simple explicit expression (10);
$\Phi^{\text{r}}$ can be computed independently in order to extract the purely
baroclinic field $\phi^{\\#}=\phi^{\dagger}-\Phi^{\text{r}}$ if needed. Also,
$\Phi^{\text{r}}$ results in a spatially trapped correction to the flow;
namely, it vanishes as $x\rightarrow\pm\infty$ and thus does not influence the
far-field energy. Eqs. (13) are supplemented with radiation conditions
ensuring that waves generated in the interior of the domain propagate outward
as plane waves. Thus, we have

$\displaystyle{\color[rgb]{0,0,0}\phi^{\dagger}=\sum\limits_{n=1}^{\infty}c^{\pm}_{n}\text{e}^{\pm\text{i}k^{\pm}_{n}x}\sin\left(\frac{n\pi
z}{h_{\pm}}\right),\ \ \text{with}\ \ k^{\pm}_{n}=\frac{n\pi}{\mu h_{\pm}},\
c^{\pm}_{n}\in\mathbb{C},\ \ \text{as}\ \ x\rightarrow\pm\infty}.$ (14)

It is useful to write down expressions for the flow fields in terms of the
amplitudes of the stream functions. From the above analysis,
$\phi=\Phi^{(0)}+\Phi^{\text{r}}+\phi^{\\#}=\Phi^{(0)}+\phi^{\dagger}$ and we
may introduce the corresponding definitions

$\displaystyle\begin{aligned}
\xi&=\Xi^{(0)}+\Xi^{\text{r}}+\xi^{\\#}=\Xi^{(0)}+\xi^{\dagger},\end{aligned}$
(15)

where $\xi$ is a placeholder for any of $u$, $v$, $w$, $b$ or $p$, and $\Xi$
is a placeholder for any of $U$, $V$, $W$, $B$ or $P$. The flow components
appearing in Eq. (2.15) satisfy the relations

$\begin{pmatrix}U^{\diamond}&V^{\diamond}&W^{\diamond}&B^{\diamond}\\\
u^{\star}&v^{\star}&w^{\star}&b^{\star}\end{pmatrix}=\text{Re}\left\\{\begin{pmatrix}-\Phi^{\diamond}_{z}&\frac{\text{i}f}{\omega}\Phi^{\diamond}_{z}&\Phi^{\diamond}_{x}&-\frac{\text{i}N^{2}}{\omega}\Phi^{\diamond}_{x}\\\
-\phi^{\star}_{z}&\frac{\text{i}f}{\omega}\phi_{z}^{\star}&\phi_{x}^{\star}&-\frac{\text{i}N^{2}}{\omega}\phi^{\star}_{x}\end{pmatrix}\text{e}^{-\text{i}\omega
t}\right\\},$ (16)

where the superscript $\star$ stands for either $\dagger$ or $\\#$, and the
superscript $\diamond$ stands for either $(0)$ or r. The relations for
$v^{\star},V^{\diamond}$ (resp. $b^{\star},B^{\diamond}$) follow from the
second equation in (5a) (resp. (5c)) and the assumption that all fields have
the same time periodicity. A similar expression can be derived for the
pressure containing additionally boundary terms. For easy reference, we
summarise the notations for the different flow fields in table 1. Plugging the
second equality of (15) into (5)–($None$) and taking into account (57)–(60)
with $i=0$, we obtain

$\displaystyle u^{\dagger}_{t}-fv^{\dagger}$
$\displaystyle=-p^{\dagger}_{x},\quad v^{\dagger}_{t}+fu^{\dagger}=0,$ (17a)
$\displaystyle w^{\dagger}_{t}$
$\displaystyle=-p^{\dagger}_{z}+b^{\dagger}+\left(1-\frac{\omega^{2}}{N^{2}}\right)B^{(0)},$
(17b) $\displaystyle b^{\dagger}_{t}+N^{2}w^{\dagger}$ $\displaystyle=0,$ (18)
$\displaystyle u^{\dagger}_{x}+w^{\dagger}_{z}$ $\displaystyle=0,$ (19)
$\displaystyle h_{x}u^{\dagger}(x,-h)+w^{\dagger}(x,-h)$
$\displaystyle=0,\quad w^{\dagger}(x,0)=0.$ (20)

In deriving (17b), we used the relation $B^{(0)}=N^{2}W^{(0)}_{t}/\omega^{2}$,
which itself derives from (16). In terms of $\xi^{\\#}$, the above equations
stay the same except for (17b), which becomes

$\displaystyle w^{\\#}_{t}=-p^{\\#}_{z}+b^{\\#}+B,$ (21)

showing that the baroclinic flow is forced by the buoyancy force created by
the barotropic flow (Garrett & Gerkema, 2007).

Notation | Field type & Equations | Related definitions
---|---|---
$\phi$, $\xi$ | Total, (8), (5)–($None$) | $\xi=\xi^{\\#}+\Xi$
$\Phi$, $\Xi$ | Barotropic, (9), (44)–($None$) | $\Xi=\Xi^{(0)}+\Xi^{\text{r}}$
$\Phi^{(0)}$, $\Xi^{(0)}$ | Hydrostatic barotropic, (10), (57)–(60) |
$\Phi^{\text{r}}$, $\Xi^{\text{r}}$ | Non-hydrostatic barotropic, (11) | $\Xi^{\text{r}}=\Xi-\Xi^{(0)}$
$\phi^{\dagger}$, $\xi^{\dagger}$ | Baroclinic-barotropic(non-hydrostatic), (12), (17)–(20) | $\xi^{\dagger}=\xi^{\\#}+\Xi^{\text{r}}$
$\phi^{\\#}$, $\xi^{\\#}$ | Baroclinic, (12), (17a), (21), (18)–(20) |
Table 1: The notations for the fields we use. The symbol $\xi$ stands for $u$,
$v$, $w$, $b$, $p$, and $\Xi$ for their capitalized versions. The relations in
the third column hold for $\xi$ and $\Xi$ replaced by $\phi$ and $\Phi$.

###  2.4 Energy equation and conversion rate

We derive here the energy equation for the above IT generation problem and use
it to define the energy conversion rates.

The dot product of (17)–(18) with
$\left(\mathbf{u}^{\dagger},b^{\dagger}/N^{2}\right)\equiv(u^{\dagger},v^{\dagger},w^{\dagger},b^{\dagger}/N^{2})$
gives

$\displaystyle\mathcal{E}^{\dagger}_{t}+\bm{\nabla\cdot}\left(p^{\dagger}\mathbf{u}^{\dagger}\right)=\left(1-\frac{\omega^{2}}{N^{2}}\right)B^{(0)}w^{\dagger},\quad\text{with}\quad\mathcal{E}^{\dagger}=\frac{1}{2}(\mathbf{u}^{\dagger})^{2}+\frac{1}{2}\frac{(b^{\dagger})^{2}}{N^{2}},$
(22)

where we have used (19) and where
$\bm{\nabla}\equiv\left(\partial_{x},0,\partial_{z}\right)$. Integrating (22)
over the domain $\Omega=[-\infty,+\infty]\times[-h,0]$, using the divergence
theorem and (20), we obtain

$\displaystyle\left(\int_{\Omega}\mathcal{E}^{\dagger}\>\mathrm{d}\Omega\right)_{t}+\left[\int_{-h}^{0}p^{\dagger}u^{\dagger}\>\mathrm{d}z\right]^{+\infty}_{-\infty}=\left(1-\frac{\omega^{2}}{N^{2}}\right)\int_{\Omega}B^{(0)}w^{\dagger}\>\mathrm{d}\Omega,$
(23)

with
$[\,\cdot\,]^{+\infty}_{-\infty}=\lim_{x\to\infty}(\cdot)-\lim_{x\to-\infty}(\cdot)$.

Applying the time average operator
$\langle\,\cdot\,\rangle=1/T\int_{0}^{T}\cdot\,\>\mathrm{d}t$, with
$T=2\pi/\omega$, and taking into account the periodicity of
$\mathcal{E}^{\dagger}(t)$, we find

$\displaystyle\begin{aligned}
C_{+}-C_{-}\stackrel{{\scriptstyle\text{def}}}{{=}}\left[\int_{-h}^{0}\langle
p^{\dagger}u^{\dagger}\rangle\>\mathrm{d}z\right]^{+\infty}_{-\infty}&=\left(1-\frac{\omega^{2}}{N^{2}}\right)\int_{\Omega}\langle
B^{(0)}w^{\dagger}\rangle\>\mathrm{d}\Omega\stackrel{{\scriptstyle\text{def}}}{{=}}C_{\text{int}},\end{aligned}$
(24)

where we have defined the energy conversion rates $C_{\pm}$ that represent the
rates at which energy is radiated at $\pm\infty$, and the total energy
convergence rate $C_{\text{int}}$ given as a volume integral. Note that the
non-hydrostatic barotropic flow does not contribute in (24) because
$u^{\dagger}=u^{\\#}+U^{r}\rightarrow u^{\\#}$ as
${\color[rgb]{0,0,0}x\rightarrow\pm\infty}$ and $\langle
B^{(0)}w^{\dagger}\rangle=\langle B^{(0)}w^{\\#}\rangle+\langle
B^{(0)}W^{r}\rangle=\langle B^{(0)}w^{\\#}\rangle$ since $B^{(0)}$ and $W^{r}$
are out of phase by $\pi/2$ (Eq. (16)). Thus, (24) remains valid with
$\dagger$ replaced by $\\#$. We proceed by expressing $C_{\pm}$ and
$C_{\text{int}}$ in terms of $\phi^{\dagger}$ and $\Phi^{(0)}$. Using
integration by parts and the fact that
$\psi^{\dagger}(x,0)=\psi^{\dagger}(x,-h)=0$ we obtain
$\int_{-h}^{0}\left\langle p^{\dagger}u^{\dagger}\right\rangle
dz=\int_{-h}^{0}\left\langle p_{z}^{\dagger}\psi^{\dagger}\right\rangle dz$.
Then, expressing $p_{z}^{\dagger}$ from (17b) using (16), we find

$\displaystyle\langle
p_{z}^{\dagger}\psi^{\dagger}\rangle=\frac{\omega^{2}-N^{2}}{\omega}\langle\text{Re}\\{\text{i}\phi^{\dagger}_{x}\text{e}^{-\text{i}\omega
t}\\}\text{Re}\\{\phi^{\dagger}\text{e}^{-\text{i}\omega
t}\\}\rangle\quad\text{as}\quad x\rightarrow\pm\infty.$ (25)

Writing $\text{e}^{-\text{i}\omega t}=\cos\omega t-\text{i}\sin\omega t$ and
noting that $\langle\cos\omega t\sin\omega t\rangle=0$, $\langle\cos^{2}\omega
t\rangle=\langle\sin^{2}\omega t\rangle=1/2$ we obtain

$\displaystyle
C_{\pm}={\color[rgb]{0,0,0}\frac{\omega^{2}-N^{2}}{2\omega}}\int_{-h}^{0}\text{Im}\\{\phi^{\dagger}\overline{\phi_{x}^{\dagger}}\\}\>\mathrm{d}z,\quad\text{as}\quad
x\rightarrow\pm\infty,$ (26)

where the overline denotes the complex conjugate and Im the imaginary part.
Similarly,

$\displaystyle C_{\text{int}}$
$\displaystyle=\left(1-\frac{\omega^{2}}{N^{2}}\right)\frac{N^{2}}{2\omega}\int_{\Omega}\Phi^{(0)}_{x}\text{Im}\\{\,\overline{\phi^{\dagger}_{x}}\,\\}\>\mathrm{d}\Omega.$
(27)

Using the radiation conditions (14) in Eq. (26), we see that $C_{+}\geq 0$
(resp. $C_{-}\leq 0$) as $x\rightarrow+\infty$ (resp. $x\rightarrow-\infty$).
Eqs. (26) and (27) provide us with two ways of calculating the total
conversion rate, either by using the the far-field baroclinic flow or by using
a barotropic-baroclinic interaction term defined in the entire fluid domain.
This fact will be used in Section 4 for validation purposes.

##  3 Modal decomposition

###  3.1 Stream function modal representation

We reformulate the IT generation problem (13)–(14) by representing
$\phi^{\dagger}$ as

$\displaystyle\phi^{\dagger}(x,z)=\sum\limits_{n=1}^{\infty}\phi_{n}(x)Z_{n}(z;x),$
(28)

where $\\{Z_{n}(z;x)\\}_{n=0}^{\infty}$ are prescribed vertical basis
functions with a parametric dependence on $x$ and
$\\{\phi_{n}(x)\\}_{n=0}^{\infty}$ are unknown complex modal amplitudes to be
determined. For the expansion (28) to be exact, the set
$\\{Z_{n}(z;x)\\}_{n=0}^{\infty}$ must be complete. In the present constant
stratification case, this set is obtained as the set of eigenfunctions of a
Sturm-Liouville problem parameterised by $x$, also called the “reference
waveguide” (Brekhovskikh & Godin, 1992),

$\displaystyle
Z_{n,zz}+{\color[rgb]{0,0,0}\frac{n^{2}\pi^{2}}{h^{2}}}Z_{n}=0,\quad
Z_{n}(0;x)=0,\quad Z_{n}(-h;x)=0.$ (29)

The eigenfunctions $Z_{n}$ are given by

$\displaystyle{\color[rgb]{0,0,0}Z_{n}(z;x)=\sin\left(\frac{n\pi
z}{h}\right),}$ (30)

and satisfy the orthogonality relation
$\int_{-h}^{0}Z_{n}Z_{m}dz=h\delta_{nm}/2$, where $\delta_{nm}$ is the
Kronecker delta. Note that (28) satisfies exactly and term by term the
boundary conditions in (13). It follows from (28) and (30) that the
$\phi_{n}$’s are defined by

$\displaystyle\phi_{n}=\frac{2}{h}\int_{-h}^{0}\phi^{\dagger}Z_{n}dz={\color[rgb]{0,0,0}\frac{2h^{2}}{\pi^{3}n^{3}}}\left(\left[\phi^{\dagger}_{zz}Y_{n}\right]_{-h}^{0}-\int_{-h}^{0}\phi^{\dagger}_{zzz}Y_{n}\>\mathrm{d}z\right),$
(31)

where $Y_{n}=\cos({\color[rgb]{0,0,0}n\pi z/h})$ and the second equality is
obtained after three integrations by parts, which are allowed provided
$\phi^{\dagger}$ is sufficiently smooth, and use of the boundary conditions in
(13). This shows that $\|\phi_{n}\|_{\infty}:=\max|\phi_{n}|=O(n^{-3})$ and
that (28) converges uniformly in this case. Similar estimates are obtained for
$\|\phi_{n,x}\|_{\infty}$ and $\|\phi_{n,xx}\|_{\infty}$ by adapting the
procedure developed in Athanassoulis & Papoutsellis (2017, Section 4) and
suffice to establish the term-wise differentiability of (28) required for the
exact modal reformulation of (13).

Griffiths & Grimshaw (2007) use a similar expansion for $u$ and derive a
series representation for $w$ by using the incompressibility and the bottom-
boundary conditions. The maximum decay rate in this case is $O(n^{-2})$ due to
the non-vanishing of $u$ on the boundaries. Note also that in contrast to
(28), the truncated version of this expansion does not satisfy term by term
the bottom boundary condition. On the other hand, Kelly (2016) uses two sets
of eigenfunctions: one for $u$ and $p$, and the other for $w$ and $b$. In this
approach, $w$ vanishes identically on the bottom which is not the case for
arbitrary $h$. Consequently, this approach should be regarded as an
approximation, see the discussion in Kelly & Lermusiaux (2016) for more
details. Despite (28) being limited to 2D flows, its major advantage is that
it satisfies the bottom boundary condition in (13) exactly and term by term,
and exhibits faster convergence in comparison with existing approaches.

###  3.2 The Coupled-Mode System

We proceed by projecting Eqs. (13)–(14) onto (30). Substituting (28) into
(13), multiplying with $Z_{n}$ and integrating over the interval $[-h(x),0]$,
we find that $\\{\phi_{n}(x)\\}_{n=1}^{\infty}$ solves the following Coupled-
Mode System (CMS), for $m\geq 1$,

$\displaystyle\phi_{m,xx}+{\color[rgb]{0,0,0}\frac{m^{2}\pi^{2}}{\mu^{2}h^{2}}}\phi_{m}+\sum\limits_{n=1}^{\infty}\left[\frac{b_{mn}h_{x}}{h}\phi_{n,x}+\left(\frac{c_{mn}h_{x}^{2}}{h^{2}}+\frac{d_{mn}h_{xx}}{h}\right)\phi_{n}\right]=2g_{m}h\left(\frac{1}{h}\right)_{xx},$
(32)

where $b_{mn}$, $c_{mn}$, $d_{mn}$ are $x$-independent coefficients given in
Appendix B, and $g_{m}=Q(-1)^{{\color[rgb]{0,0,0}m+1}}/(m\pi)$. Substituting
Eq. (28) into (14), we obtain

$\displaystyle\sum\limits_{n=1}^{\infty}\phi_{n}(x)Z_{n}(x;z)=\sum\limits_{n=1}^{\infty}c_{n}^{\pm}\exp(\pm\text{i}k^{\pm}_{n}x)\sin\left({\color[rgb]{0,0,0}\frac{n\pi
z}{h_{\pm}}}\right)\quad\text{as}\quad x\rightarrow\pm\infty.$ (33)

Multiplying both sides by $Z_{q}$, $q\geq 1$, integrating over $[-h,0]$ and
taking into account that $Z_{n}\rightarrow\sin({\color[rgb]{0,0,0}n\pi
z/h_{\pm}})$ as $x\rightarrow\pm\infty$, we obtain
$\phi_{n}=c_{n}^{\pm}\exp(\pm\text{i}k^{\pm}_{n}x)$, and therefore

$\displaystyle\phi_{n,x}{\color[rgb]{0,0,0}\mp}\text{i}k^{\pm}_{n}\phi_{n}=0,\quad\text{as}\quad
x\rightarrow{\pm}\infty.$ (34)

Thus, Eqs. (13)–(14) are exactly reformulated as the CMS (32)–(34). Solving
the latter, we may reconstruct the solution of the former by using (28). The
energy conversion rates $C_{\pm}$ are evaluated by (26):

$\displaystyle
C_{\pm}=\frac{\omega^{2}-N^{2}}{2\omega}\frac{h}{2}\sum\limits_{n=1}^{\infty}\text{Im}\\{\phi_{n}^{\dagger}\overline{\phi^{\dagger}}_{\\!\\!\\!n,x}\\}=\pm\frac{N^{2}-\omega^{2}}{4\omega\mu}\pi\sum\limits_{n=1}^{\infty}n\phi_{n}^{\dagger}\overline{\phi^{\dagger}}_{\\!\\!\\!n}\quad\text{as}\quad
x\rightarrow{\pm}\infty,$ (35)

where Eqs. (34) and the definition of $k_{n}^{\pm}$ in (14) have been used to
obtain the second equality. The purely baroclinic field $\phi^{\\#}$ is
computed using $\phi^{\\#}=\phi^{\dagger}-\Phi^{\text{r}}$. For the
computation of $\Phi^{\text{r}}$, we apply the same modal decomposition to
(11). The resulting CMS is given by (32) with $\mu^{2}$ replaced by
$-\mu_{0}^{2}$ and vanishing conditions at infinity, instead of (34); if
$h_{0}/L\ll 1$, one could use instead the approximate asymptotic expression in
(51). In the following subsection, we derive a perturbative solution of the
CMS (32)–(34) for infinitesimal topography. For arbitrary topographies, we
solve the CMS numerically (see subsection 3.4).

###### Remark 1

The body-forcing formulation (12) and its coupled-mode reformulation (32) are
valid for non-hydrostatic conditions. If the hydrostatic approximation (HA) is
invoked for both the baroclinic and the barotropic flow ($\omega^{2}\ll N^{2}$
and $h_{0}^{2}\ll L^{2}$), Eq. (12) becomes (13) with $\mu^{2}$ replaced
$\mu^{\prime 2}=(\omega^{2}-f^{2})/N^{2}$ and $\phi^{\dagger}$ is interpreted
as a purely baroclinic response, that is, $\Phi^{r}=O(h^{2}/L^{2})$ may be
neglected (Appendix A) (Garrett & Gerkema, 2007). The CMS (32) changes
accordingly. Under this assumption and using (16), Eq. (24) reduces to the
energy equation derived in Gerkema et al. (2004) starting from the hydrostatic
governing equations. If instead we assume $\omega^{2}\ll N^{2}$ and
$h_{0}^{2}/L^{2}=O(1)$ then $\Phi^{r}$ is not negligible but does not affect
the energy flux at infinity. In the reverse situation, $\omega^{2}<N^{2}$ and
$h_{0}^{2}/L^{2}\ll 1$, Eqs. (12) and (32) hold as is but once again
$\phi^{\dagger}$ is interpreted as a purely baroclinic response. For more
details on the relevance of the HA we refer to the discussion in Garrett &
Gerkema (2007). Unless otherwise stated, we do not invoke the HA because the
obtained mathematical simplification is minute and the precise quantification
of the induced error for different values of $\mu$ is out of the scope of this
work.

###  3.3 Infinitesimal topography solution

Let $h(x)=h_{0}-\epsilon\,r(x)$ with $|\epsilon|\ll 1$ and $r=O(1)$, for some
characteristic depth $h_{0}$. Introducing the asymptotic expansion
$\phi_{m}=\sum_{i=0}^{K}\epsilon^{i}\phi_{m}^{(i)}$ and Taylor-expanding $1/h$
in (32), we deduce that $\phi_{m}^{(0)}=0$ and that $\phi_{m}^{(1)}$ solves

$\displaystyle\phi^{(1)}_{m,xx}+\ell_{m}^{2}\phi^{(1)}_{m}=2g_{m}\frac{r_{xx}}{h_{0}},\quad\text{for}\quad
m\geq 1,$ (36)

with $\ell_{m}={\color[rgb]{0,0,0}m}\pi/(h_{0}\mu)$, and the radiation
conditions in (34) with $k_{m}^{+}=k_{m}^{-}=\ell_{m}$. Substituting the
solution (64) in (26) we obtain the conversion rate for weak topography,

$\displaystyle
C^{\text{WTA}}=\frac{F_{0}}{h_{0}^{2}}\sum\limits_{n=1}^{\infty}|\hat{r}(\ell_{n})|^{2}\frac{n\pi^{2}}{(\mu
h_{0})^{2}}\,\,\,\,\,\text{with}\,\,\,\,\,F_{0}=\frac{\left[(N^{2}-\omega^{2})(\omega^{2}-f^{2})\right]^{1/2}}{2\pi\omega}U{{}_{0}}^{2}h_{0}^{2},$
(37)

where $\hat{r}(\xi)=\int_{-\infty}^{+\infty}\exp(-\text{i}x\xi)r(s)ds$ is the
Fourier transform of $r$ (Appendix C). This is the formula given in St.
Laurent et al. (2003) (up to multiplication with $\rho_{0}$), which was first
derived in Llewellyn Smith & Young (2002) in the case of hydrostatic internal
waves. Khatiwala (2003) derived a different formula starting from a multi-
frequency representation of the response fields in terms of Bessel functions.
In the acoustic limit (recall discussion following Eq. 4), the dominant
contribution to the far-field energy comes from the fundamental frequency,
say, $\omega_{0}$ (Bell, 1975). The quantity $J_{1}(U_{0}k_{j1}/\omega_{0})$
appearing in Khatiwala (2003)’s Eq. (28), where $J_{1}$ is the order-one
Bessel function of the first kind, and $U_{0}$ and $k_{j1}$ are his barotropic
tidal amplitude and wavenumber for the $j^{\text{th}}$ mode of the fundamental
frequency, respectively, may be replaced by its asymptotic value
$U_{0}k_{j1}/(2{\color[rgb]{0,0,0}\omega_{0}})$ leading to our Eq. (37) with
$\omega=\omega_{0}$.

###  3.4 Numerical solution for arbitrary topography

We truncate the infinite CMS (32)–(34) by keeping the first $M$ equations and
replacing the infinite domain by a finite interval $X=[x_{L},x_{R}]$ of length
$L_{X}=x_{R}-x_{L}$. We discretise $X$ with a uniform spacing $\delta x$,
$\\{x_{i}$, $i=1,N_{X}\\}$ and we approximate the derivatives using
$4^{th}$-order finite differences up to the boundary points $x=x_{L},x_{R}$.
The corresponding formulae can be found in Papoutsellis et al. (2019, Appendix
C). We thus obtain a sparse square linear system of dimension $(N_{X}M)^{2}$
for the grid values of each modal amplitude $\phi_{n}(x_{i})$,
$i=1,\dots,N_{X}$, $n=1,\dots,M$ which we solve by means of a $LU$
decomposition. Using this solution, we reconstruct the field $\phi^{\dagger}$
by means of a truncated version of (28). We then compute the baroclinic fields
and conversion rates using Eqs. (16) and (26). We also provide an estimation
of the free-surface elevation due to the IT motion given by
$\eta=[p^{\\#}]_{z=0}/g$, where $[p^{\\#}]_{z=0}$ is the pressure induced on
$z=0$ by the baroclinic flow (Appendix D).

##  4 Convergence, accuracy and singularity formation

In this section, we introduce the topographic profiles we use in our
calculations and present results showing the good performance of our semi-
analytical solution.

###  4.1 Topographic profiles

We consider two ridge profiles, namely the “Gaussian" $h=h_{0}-h_{\text{G}}$
and the “bump", $h=h_{0}-h_{\text{B}}$, where

$\displaystyle
h_{\text{G}}=\Lambda\exp\left(-\frac{x^{2}}{2L^{2}}\right)\quad\text{and}\quad
h_{\text{B}}=\Lambda\exp\left(1-\frac{1}{1-x^{2}/L^{2}}\right)\mathds{1}_{(-L,L)},$
(38)

with $\mathds{1}_{(-L,L)}=1\,\,\text{in}\,\,(-L,L)$ and
$\mathds{1}_{(-L,L)}=0$ otherwise. Note that the “bump" profile has a compact
support, in contrast with the Gaussian, and all its derivatives are continuous
at $x=\pm L$; see figure 2. We also consider the case of a shelf connecting
two different constant depths. The shelf profile is the same as in Griffiths &
Grimshaw (2007, Section 5), namely, $h=h_{\text{S}}$ with

$\displaystyle h_{\text{S}}=\left\\{\begin{aligned} &h_{-}\quad\text{for}\quad
x\leq 0,\\\ &h_{-}+(h_{+}-h_{-})\sin^{2}(\pi x/2L)\quad\text{for}\quad 0\leq
x\leq L,\quad\text{and}\\\ &h_{+}\quad\text{for}\quad x\geq
L.\end{aligned}\right.$ (39)

We recall from subsection 2.1 that two important parameters will be
considered, namely, the relative topography height $\delta$ (Eq. 1) and the
criticality $\varepsilon$ (Eq. 3) of the bottom slope. We also recall that,
for illustration purposes, $(\omega,U_{0})$ corresponds to a typical $M_{2}$
tide, i.e., $(\omega,U_{0})=(1.4\times 10^{-4}\
\text{s}^{-1},0.04\,\text{m\,s}^{-1})$ and that the Brunt-Väisälä and Coriolis
frequencies are $N=1.5\times 10^{-3}\ \text{s}^{-1}$ and $f=10^{-4}\
\text{s}^{-1}$, respectively. For all ridges, the depth far from the
topography is $h_{0}=3000$ m.

Figure 2: The three topographic profiles we consider (in the Shelf case, we
have used $h_{-}=0$, $h_{+}=\Lambda$).

###  4.2 Results

We first examine the rate of decay of $\phi_{n}$ with $n$ obtained by solving
the CMS. It is known that the baroclinic field becomes singular in the
infinite-depth horizontally periodic case as $\varepsilon\rightarrow 1$
(Balmforth et al., 2002) and for a finite-depth shelf for
$\varepsilon\gtrapprox 1$ (Griffiths & Grimshaw, 2007). In these works, the
presence of the singularity is identified by a decrease in the rate of decay
of the coefficients in the respective modal solutions. We verify that this is
also the case for the present CMS. We use the “bump" profile with
$\delta=0.5$, $\varepsilon\in[0.9(0.05)1.5]$ and solve the CMS for
sufficiently large values of the parameters $(M,N_{X},L_{X})$ (see Section
3.4) in order to ensure that the numerical solution does not depend on them.
We show the results on $\|\phi_{n}\|_{\infty}$ in Figure 3$(b)$. We see that
the decay rate drops from $n^{-3}$, the theoretically expected rate in the
case of a smooth solution (Section 3.1), to $n^{-3/2}$ as $\varepsilon$
approaches 1. The latter rate suggests the presence of a square root
singularity on $\phi^{\dagger}$ and an inverse square root singularity on its
derivatives, i.e., on the velocities (Salem, 1939; Raisbeck, 1955).

Next, we consider the normalised error of the energy equation (24),
${\color[rgb]{0,0,0}E}=\lvert
C_{+}-C_{-}-C_{\text{int}}\rvert/{\color[rgb]{0,0,0}F_{0}}$, where $C_{\pm}$
are calculated via (26) and $C_{\text{int}}$ via (27). $C_{\pm}$ depends only
on boundary values of $\phi_{n}$’s, whereas $C_{\text{int}}$ depends on their
values in the entire computational domain. For an exact solution,
${\color[rgb]{0,0,0}E}$ would be zero. Therefore, monitoring
${\color[rgb]{0,0,0}E}$ as a function of $(M,N_{X})$ is a good indication of
the accuracy of the numerical solution. We use $\delta x=L_{M}/s$, where
$L_{M}$ is the horizontal wavelength of the $M^{\text{th}}$ internal wave mode
over the ridge and $s=4,6,8,10$ and $12$. We show results for
$\varepsilon={\color[rgb]{0,0,0}0.7}$ and
$\varepsilon={\color[rgb]{0,0,0}1.0}$ in figure 3$(b)$. We obtain the expected
$s^{-4}$ decay in both cases, verifying the $4^{\text{th}}$-order accuracy of
the spatial discretisation scheme. For $\varepsilon={\color[rgb]{0,0,0}0.7}$,
the error decays rapidly as $M^{-4}$. As an example, we mention that for
$(M,s)=(30,6)$, $E=3.1\times 10^{-7}$. For
$\varepsilon={\color[rgb]{0,0,0}1.0}$, the error decays as
${\color[rgb]{0,0,0}M^{-1/2}}$ for $s\geq 6$ and $M\geq 80$, demonstrating the
slower convergence of the modal solution when the underlying field becomes
singular. Nevertheless, the numerical solution accurately satisfies the energy
balance even in this case. For example, for $(M,s)=(120,10)$, $E=1.6\times
10^{-6}.$

A final issue that must be addressed is the truncation of the infinite domain.
In other words, we must ensure that the lateral boundary conditions applied at
the ends of the computational domain $X$ are effective as radiation conditions
and the size of $X$ does not affect the solution. If $h$ (resp., $h_{x}$ in
the shelf case) is not compactly supported, then the length of the
computational domain, $L_{X}$, is chosen so that $h_{x}$ at the boundaries is
negligible and further increasing $L_{X}$ does not change the calculated
conversion rate. For compactly supported $h$ (resp., $h_{x}$), we have
examined the sensitivity of the conversion rate on the choice of $L_{X}$ and
found that choosing $X$ as the support is sufficient for a convergent
solution.

Concerning the choice of $M$ and $N_{X}$ (or $s$), some remarks are in order.
As $\delta$ decreases for constant $\varepsilon$, the horizontal resolution
must increase to adequately represent the topography of the ridge. On the
other hand, as $\delta$ increases, the topography becomes longer and the
computational domain must increase. Moreover, as $\varepsilon$ increases
beyond the subcritical regime, both $M$ and the horizontal resolution must
increase to achieve a good representation of the singular solution. Thus, the
extent of the $(\varepsilon,\delta)$ values that can be considered depends on
the given computational resources. In this work, we let $\delta$ vary in
$[0.1,0.9]$, which is sufficient for our purposes. In this range, the choices
$(M,s)=(64,6)$ and $(M,s)=(128,{\color[rgb]{0,0,0}12})$ lead to convergent
solutions for $\varepsilon\leq 1$ and $1\leq\varepsilon\leq 2$, respectively,
while for $\varepsilon>2$ we increase $(M,s)$ further until the solution
becomes independent of them.

Figure 3: $(a)$: Decay of $\|\phi_{n}\|_{\infty}$ in log-scale for the case of
a bump topography with $\delta=0.5$ and $\varepsilon={\color[rgb]{0,0,0}0.7}$,
$0.85$, $0.90$, $0.95$, and $1.00$. $(b)$: Normalised absolute error of the
energy balance equation, $E$ = $|C_{+}-C_{-}-C_{\text{int}}|/F_{0}$, in log-
scale as a function of the order of truncation $M$ and the spatial
discretisation parameter $s=L_{M}/\delta x$ for
$\varepsilon={\color[rgb]{0,0,0}0.7}$ (blue surface) and
$\varepsilon={\color[rgb]{0,0,0}1.0}$ (orange surface).

##  5 Visualisation of flow fields

###  5.1 Gaussian ridge

Here we consider solutions of the CMS for a subcritical
($\varepsilon=0.{\color[rgb]{0,0,0}8}$) and a supercritical
($\varepsilon=1.{\color[rgb]{0,0,0}2}$) Gaussian ridge. In figure 4, we show
the purely baroclinic stream function $\phi^{\\#}$, the energy density
$\mathcal{E}^{\\#}=(\mathbf{u}^{\\#})^{2}/2+(b^{\\#})^{2}/2/N^{2}$, the
reconstructed free-surface elevation, and the body-forcing term
$\Phi^{(0)}_{xx}$.

We clearly observe the beam-like structure of ITs, which is finer and more
intense in the supercritical case. In this case, the solution approximates a
field that is continuous along the beams, where the energy density attains
very large values (theoretically infinite). The free surface also transforms
from smooth to cusp-like as $\varepsilon$ exceeds $1$. Such fine-scale large-
amplitude features are regularised in the ocean via non-linearity and
viscosity, which are beyond the present theory. Nevertheless, such flow-field
calculations are readily obtained using the CMS and give us useful information
for observable quantities. For example, in the supercritical case, the global
maximum on the horizontal velocity is attained at the point of intersection of
the beams above the ridge. The following maxima are attained at the first
points of reflection at the free surface. Right at these points, the free
surface attains its maximum slope.

Interestingly, the solution varies considerably between the sloping regions
and the flatter regions. This is due to the body-forcing term in (13), which
couples the baroclinic motion with the non-uniform barotropic current and
scales with ${\color[rgb]{0,0,0}(1/h)_{xx}}$ (Figures 4$(c)$ and $(f)$). This
effect becomes more apparent if we write the PDE in (13) as a first-order
system,

$\displaystyle\left(\phi^{\dagger}_{x}\mp\frac{1}{\mu}\phi^{\dagger}_{z}\right)_{x}\pm\frac{1}{\mu}\left(\phi^{\dagger}_{x}\mp\frac{1}{\mu}\phi^{\dagger}_{z}\right)_{z}$
$\displaystyle=-\Phi^{(0)}_{xx}.$ (40)

Note that the above equation also holds with $\phi^{\dagger}$ replaced by
$\phi^{\\#}$ and $\Phi^{(0)}$ by $\Phi$. Eq. (40) shows that the solution is
not constant along a given characteristic unless $\Phi^{(0)}_{xx}=0$. This is
captured by our solution as an adjustment of the iso-energy curves near the
ridge crest, figures 4$(b)$ and $(e)$. In the insets of these figures, we also
plot the field $\phi^{\\#}_{x}+\phi^{\\#}_{z}/\mu$ in order to highlight the
variation of the beams.4). This effect seems to have gone unnoticed in the
literature, although it is clearly present in the calculations of Stashchuk &
Cherkesov (1991). It can also be inferred from the calculations of Baines
(1973), even though he does not visualise flow fields in the physical domain.

Figure 4: Gaussian ridges with $\delta=0.5$;
$\varepsilon={\color[rgb]{0,0,0}0.8}$ (left) and
$\varepsilon={\color[rgb]{0,0,0}1.2}$ (right). From top to bottom, non-
dimensional baroclinic stream function $\psi^{\\#}/Q$, energy density
$\mathcal{E}^{\\#}/U_{0}^{2}$, and body forcing term
$\varepsilon^{2}b^{2}\Phi_{xx}^{(0)}/Q$. Free-surface elevation $3\times
10^{4}\eta/h_{0}$ is also shown in light blue. Black lines correspond to iso-
energy curves of $\mathcal{E}^{\\#}/U_{0}^{2}=0.36$. Magenda dashed lines
correspond to characteristic lines of slope $\pm\mu^{-1}$. Vertical dashed
lines correspond to the critical points of the topography. Insets show the
field $\phi^{\\#}_{x}+\phi^{\\#}_{z}/\mu$ in the region
$[-30,+7]\times[-1.05,0]$. The colormaps used in this and the subsequent
figures are from Thyng et al. (2016).

###  5.2 A non-radiating ridge

For further validation, we reproduce a case investigated by Maas (2011,
Section 4.2, “Ridge”). Maas considered the dimensionless version of Eq. (13),
obtained by the scaling
$(x,z,\phi^{\dagger},h)=(L\tilde{x},L/\mu\tilde{z},Q\tilde{\phi^{\dagger}},L/\mu\tilde{h})$,
and showed that it has an exact solution, denoted here by $\phi^{ex}$, that
radiates no energy at infinity for a specially constructed ridge; see also
Wunsch & Wunsch (2022) for an alternative construction of non-radiating
topographies. For the case in Eq. (4.15) of Maas (2011), which is very
similar to a Gaussian ridge with $\Lambda=0.19$ and $L=e/\sqrt{2}$ in our
Eq.(38), the exact non-radiating solution $\phi^{ex}$ is shown in Figure
5$(a)$. Our calculation, $\phi^{\dagger}$, is in agreement with $\phi^{ex}$ as
shown in figure 5$(b)$ where the difference $\phi^{ex}-\phi^{\dagger}$ is
plotted. The field $\phi^{ex}$ (or $\phi^{\dagger}$) is a combination of a
baroclinic and a non-hydrostatic barotropic response trapped near the ridge.
The purely baroclinic response is shown in Figure 5$(c)$. The hydrostatic
body-forcing term is shown in Figure 5$(d)$. The calculated non-dimensional
energy conversion rate is on the order $10^{-17}$.

Figure 5: The non-radiating ridge of Maas (2011); Non-dimensional $(a)$ exact
baroclinic response $\phi^{ex}$, and $(b)$ difference
$\phi^{ex}-\phi^{\dagger}$ where $\phi^{\dagger}$ is calculated using the CMS,
$(c)$ purely baroclinic response $\phi^{\\#}$ calculated with the CMS, and
$(d)$ body-forcing term $\Phi^{(0)}_{xx}$. All variables are dimensionless.

##  6 Energy conversion rates

In this section, we calculate the total energy conversion rate $C=C_{+}-C_{-}$
(Eq. 24) for the topographies presented in Section 4. Our primary objective is
to examine the dependence of $C$ on $\varepsilon$. We also examine the region
of validity of the WTA prediction (Eq. 37).

We first consider Gaussian profiles with $\delta=0.1$, $0.5$, $0.9$, and
$\varepsilon$ ranging from $0.1$ to $5$. We show in Figure 6$(a)$ the non-
dimensional calculated conversion rate $C/F_{0}$, where $F_{0}$ is given in
(37), using a log-log plot in order to capture large relative variations. In
the inset, we show $C/C_{KE}$, where $C_{KE}$ is the knife-edge prediction of
(Llewellyn Smith & Young, 2003, Eq. (2.12)), using a linear plot to focus on
the largest $\varepsilon$ values. For $\delta=0.1$, $C$ increases rapidly up
to about $\varepsilon=0.2$ and more slowly afterwards. The WTA prediction is
in excellent agreement with the full solution in the entire subcritical
regime. For $\delta=0.5$ and $0.9$, we clearly observe abrupt drops in $C$ for
discrete values of $\varepsilon$. However, note that $C/F_{0}$ never strictly
drops to zero. Nevertheless, calculations for values of $\varepsilon$ closely
straddling a local conversion minimum in the case $\delta=0.5$ reached
conversion values approaching machine precision (not shown). The WTA presents
strong qualitative and quantitative differences and clearly fails to predict
these abrupt drops. As $\varepsilon$ exceeds 1, the transition looks smooth
for all $\delta$. For $\delta=0.1$, $C$ is initially increasing until
$\varepsilon\approx 2.0$, after which it remains very close to $C_{KE}$.
Khatiwala (2003) also reports this increase for a small-height Gaussian
profile up to $\varepsilon=1.6$ in his non-linear calculations. For
$\delta=0.5$, $C$ monotonically increases with $\varepsilon$, approaching
$C_{KE}$. For $\delta=0.9$, a local minimum is attained at about
$\varepsilon=1.4$, after which $C$ increases, also approaching $C_{KE}$ (inset
of figure 6$(a)$). It should also be noted that we obtained a similar trend
for the Witch of Agnesi profile (not shown). In figure 6$(a)$, we also show
that the calculated conversion rate is in overall agreement with the one
obtained by using the Green’s function method (Echeverri & Peacock, 2010;
Mercier et al., 2012). The differences appearing for small $\varepsilon$ in
the case $\delta=0.9$ are due to the domain truncation needed to avoid
singularities in the Green-function solution. Furthermore, the modal
convergence being slower than the CMS approach, these values are also less
precise111The horizontal domain considered describes between $99\%$ and
$99.5\%$ of the variation of height of the topography.. We refer to section 7
for a further discussion.

We also performed the same calculations as above for the bump profile, figure
6$(b)$. For $\delta=0.1$, $C/F_{0}$ reaches a single local minimum, which is
adequately predicted by the WTA. For larger $\delta$, similarly as before, we
observe local extrema in the full solution, although in this case there are
many more in any given range of $\varepsilon$ values. The WTA gives
qualitatively similar results, predicting certain local extrema with values
comparable to the full solution, but not at the right locations. A similar
behaviour is also reported in (Vlasenko et al., 2005, Section 2.3.1) in terms
of the amplitude of the first mode for a small-height sinusoidal bump. As in
the Gaussian profile, the transition to the supercritical regime is smooth and
for larger $\varepsilon$ and all $\delta$ the conversion rates approach
$C_{KE}$. In the extreme case $\delta=0.9$, local extrema in $C$ are obtained
up to $\varepsilon=3.2$. After this value, $C$ increases slowly towards
$C_{KE}$. In this case, we have extended our calculations up to
$\varepsilon=11$ and found that $C/C_{KE}\approx 0.98$.

In order to assess the region of validity of the WTA prediction, we performed
more than 20,000 calculations per ridge in the parameter space
$(\varepsilon,\delta)=[0.1,1.5]\times[0.1,0.9]$ and calculated the relative
error $|(C-C^{\text{WTA}})/C|$. We show the results in the right panels of
Figure 6. Note that in the case of the Gaussian ridge, there is a hatched
region (Figure 6$(b)$, upper left corner) that corresponds to a practically
negligible conversion rate $C/F_{0}$ on the order of $10^{-8}$. The
corresponding values obtained by the WTA reach machine precision. Based on
these figures, it is clear that in both cases we may identify coherent regions
in the parameter space for which the WTA is valid according to some relevant
criterion (e.g., conversion rates that differ by less than 10% with respect to
the full solution). Interestingly, these regions include unexpectedly large
values of $\varepsilon$ and $\delta$. However, it should be stressed out that
the region of validity of the WTA depends strongly on the type of the ridge
and is clearly smaller in the case of the bump profile.

We close this section by demonstrating the applicability of the CMS in the
case of a shelf profile. We consider the topography in Eq. (39) with
$h_{-}=2000$ m and $h_{+}=1000$ m ($\delta=0.5$) for which calculations based
on another modal decomposition method are available in Griffiths & Grimshaw
(2007), GG07, for $\varepsilon\in[0.1,4]$. As in that work, we use the
hydrostatic version of the CMS (Remark 1). We plot our results on the
conversion rates $C_{\pm}$ for $\varepsilon\in[0.1,5]$ in figure 7, together
with the results digitised from Fig. 11 of GG07. The results from the two
methods are in agreement for all $\varepsilon$. In the inset of figure 7, we
additionally compare in a linear scale our calculations with the conversion
rate for a step profile, $C_{st}$, obtained by using eigenfunction matching at
the topographic discontinuity (St. Laurent et al., 2003)222We obtained this
result from Eqs. (32)-(36) of St. Laurent et al. (2003) for 2000 modes. Note
however that for $\delta=0.5$, the matrix in Eq. (32) is singular and a value
$\delta=0.500001$ is used here.. As is also noted in GG07, $C$ is slowly
increasing up to $\varepsilon=4$. This increase is also reported in the early
calculations of Craig (1987) for a linear step up to $\varepsilon=2$. Our
calculations for $\varepsilon>4$ suggest that $C$ slowly approaches the value
$C_{st}$ from below. In fact, we have observed that $C/C_{st}$ goes from 0.98
for $\varepsilon=5$, to 0.99 for $\varepsilon=10$.

Figure 6: $(a)$ Calculated non-dimensional energy conversion rate $C/F_{0}$ in
log-log scale as a function of $\varepsilon$ (solid lines) for $\delta=0.1$,
$0.5$ and $0.9$, for a Gaussian ridge. Dashed lines correspond to the WTA
prediction $C^{\text{WTA}}/F_{0}$. Circles correspond to the results obtained
by using the Green’s function method. The inset corresponds to the same
calculation in linear scale, normalised by $C_{KE}$ which denotes the knife-
edge prediction of Llewellyn Smith & Young (2003). $(b)$ Relative error
$|(C-C^{\text{WTA}})/C|$ for
$(\varepsilon,\delta)\in[0.1,1.5]\times[0.1,0.9]$. The same calculations for
the bump ridge are shown in $(c)$ and $(d)$. The hatched region in $(b)$
corresponds to topographies for which conversion is negligible for both the
full solution and the WTA.

Figure 7: Non-dimensional conversion rates $C_{-}/Q^{2}$ (blue solid line) and
$C_{-}/Q^{2}$ (blue dashed line) for the shelf profile with $h_{-}/h_{+}=0.5$.
Circles and crosses correspond to the calculations of Griffiths & Grimshaw
(2007). The inset corresponds to the same calculation in linear scale and the
dotted dashed blue line shows the total calculated conversion rate,
$C=C_{+}-C_{-}$, normalised by St. Laurent et al. (2003)’s prediction for a
step profile, $C_{st}$.

##  7 Discussion on the Green’s function method

We briefly discuss here how our method compares with an existing method for IT
calculations based on Robinson (1969)’s Green’s function for the homogeneous
internal wave equation in a uniform strip (Pétrélis et al., 2006; Echeverri &
Peacock, 2010). This approach is limited to domains of the same depth at
infinity, say $h_{0}$, and the starting point is to write the total flow as
$\phi^{\text{G}}(x,z)=\Phi_{0}(z)+\vartheta(x,z)$, where $\Phi_{0}=-Qz/h_{0}$
represents the barotropic flow corresponding to the uniform strip
$[-\infty,+\infty]\times[-h_{0},0]$. From (8), the field $\vartheta$ solves

$\displaystyle\mathcal{L}_{\mu}\vartheta=0,\quad\vartheta(x,0)=0\quad\text{and}\quad\vartheta(x,-h)=Q(1-h/h_{0}),$
(41)

together with conditions at $x\rightarrow\pm\infty$. Next, in order to solve
(41), $\vartheta$ is expressed as an ansatz defined as the integral of the
product of the Green’s function evaluated at $z=-h$ with an unknown
distribution $\gamma(x)$,

$\displaystyle\vartheta(x,z)=\int_{-a}^{a}\gamma(x^{\prime})G(x,x^{\prime},z,-h(x^{\prime}))\>\mathrm{d}x^{\prime},$
(42)

where it is additionally assumed that $-h_{0}+h$ is supported on $[-a,a]$.
Note that (42) implies that $\vartheta$ is defined on $[-a,a]\times[-h_{0},0]$
therefore trench topographies (for which $h\geq h_{0}$) are also excluded.
Expression (42) satisfies the first two equations in (41) as well as the
radiation conditions as $a\rightarrow+\infty$. The third equation in (41)
yields

$\displaystyle\int_{a}^{-a}\gamma(x^{\prime})G(x,x^{\prime},-h(x),-h(x^{\prime}))\>\mathrm{d}x^{\prime}=Q(1-h/h_{0}),$
(43)

which must be solved for $\gamma(x)$. Eq. (43) is a Fredholm integral
equation of the first kind, which is known to be ill-posed (see e.g. Groetsch
(1984, Section 1.1), Groetsch (1990)). This has unpleasant consequences, as is
recognised by Pétrélis et al. (2006) and Echeverri & Peacock (2010); the
linear system obtained via the series representation of the discontinuous
Green’s function and the spatial discretisation of Eq. (43) becomes ill-
conditioned or even singular if $a$ is large or if $h_{x}\approx 0$ in some
region, among other special cases (Echeverri et al., 2011). Moreover, as noted
by Maas (2011), the physical meaning of $\vartheta$ is clear only as
${\color[rgb]{0,0,0}x\rightarrow\pm\infty}$, in which case
$\vartheta(x,-h)\rightarrow 0$ and $\vartheta$ may be interpreted as the far-
field baroclinic response. Therefore, the present CMS is more general and
numerically advantageous. There exist, however, two issues that should be
discussed further.

The first issue concerns the  local dynamics above the topography. As shown in
Section 5, the baroclinic wave field varies significantly in regions where the
body-forcing term is important (recall Eq. (40)). On the other hand, the
ansatz (42) postulates that $\vartheta$ can be represented as a superposition
of sources, all of which ignore the interaction between the baroclinic and the
non-uniform barotropic background flow. The intensities of the sources are
instead determined by an ill-posed equation stemming from the bottom boundary
condition in (41). We can reconcile the two approaches by verifying that the
purely baroclinic field, satisfying Eq. (12), is recovered by
$\phi^{\\#}=\Phi_{0}+\vartheta-\Phi$, with $\Phi=\Phi^{(0)}+\Phi^{r}$.

The second issue is associated with the the so-called “knife-edge limit”
established in Pétrélis et al. (2006), see also Echeverri & Peacock (2010).
The conversion rate obtained via the ansatz (42) tends to this limit because
(42) is actually a direct extension of the knife-edge solution of Llewellyn
Smith & Young (2003) in the smooth topography case. In fact, Pétrélis et al.
(2006) state that only in the knife-edge case is “the kernel of the integral
equation (…) so singular that the resulting linear system is well
conditioned”.

The knife-edge solution corresponds to a discontinuous topography
($\varepsilon=\infty$) for which the body-forcing term in Eq. (13) (or,
equivalently, in Eq. (12)) is not well-defined because Eq. (9) requires the
domain to be smooth. Moreover, in the knife-edge formulation, the amplitude of
the horizontal barotropic velocity is constant, say $U_{0}$. This implies that
the total flux through a vertical cross-section is $U_{0}h_{0}$ everywhere
except at the position of the knife-edge where it equals
$U_{0}(h_{0}-\Lambda)$, with $\Lambda$ the height of the knife-edge. This is
not consistent with the key assumption in the body-forcing formulation that
the barotropic flow accounts for a constant total flux through any vertical
cross-section. Therefore, the present numerical evidence showing that the
conversion rate, obtained by the body forcing formulation, approaches the
knife-edge (or step) limit for large $\varepsilon$ seems unexpected. Although
this limiting behaviour of the singular solution is interesting, its
theoretical analysis is beyond the scope of this paper and remains a topic for
future investigation.

##  8 Conclusions

We developed a new infinite CMS describing the generation of linear ITs in 2D
using Garrett & Gerkema (2007)’s body-forcing formulation and an exact, local
eigenfunction expansion of the stream function. In the weak topography limit,
we recover from this CMS the classical formula of Llewellyn Smith & Young
(2002) for the conversion rate. For general topographies, we solve the
truncated CMS numerically using fourth-order finite differences and
demonstrate the convergence and accuracy of the semi-analytical scheme. The
formation of singularities is detected as a decrease in the rate of decay of
the modal amplitudes. We additionally derive the energy equation for the body-
forcing formulation and show that our solution verifies it with very good
accuracy, even when the field becomes singular. Further, we show how one can
obtain the purely baroclinic response and also propose a method to estimate
the free-surface elevation induced by the baroclinic motion within the rigid-
lid approximation.

By reconstructing the flow field, we showed that the interaction between ITs
and the background non-uniform barotropic flow affects the dynamics locally
around the sloping topography. We calculated the conversion rate for two
different ridges. For the commonly used Gaussian profile, we find distinct
points of practically zero energy radiation in the subcritical regime for
$\delta>0.1$. We have checked that similar null-points exist also for the
Witch of Agnesi profile (not shown). Our calculations verify the hypothesis of
Maas (2011) that non-radiating topographies are common and possible for
subcritical topographies with sufficiently large values of $\delta$. For the
compactly supported bump profile, local extrema of the conversion rate exist
for all $\delta$ and persist even in the supercritical regime for large
$\delta$. These local minima do not correspond to null points, but rather to
points of weak radiation. These differences are due to the qualitative
differences between a Gaussian and a bump profile. The baroclinic response
depends not only on the profile, but also on its first and second derivatives
through the body forcing term Eq. (13), which can differ significantly for an
infinite-support decaying topography and a compactly supported one. In the
context of the WTA approximation, the different behaviour of the conversion
rates can also be explained by the different decays of the Fourier transforms
of the topographies. We have also tested the Witch of Agnesi profile and a
compactly supported polynomial ridge and found that these two different
behaviours persist. Thus, it appears that the trend of the conversion rate for
a Gaussian ridge (resp. bump ridge) is typical for infinite-support (resp.
compactly supported) ridges. We have also quantified the error of the WTA with
respect to the full linear solution. Our general conclusion is that the WTA is
adequate even for large $\varepsilon$ for sufficiently small $\delta$ but the
region of its validity strongly depends on the type of topography. Similarly,
our results show that the knife-edge limit agrees with the linear singular
solution for sufficiently large $\varepsilon$, although the value of
$\varepsilon$ for which this limit is attained depends on $\delta$ and the
topography.

Our comparisons show that the CMS and the Green’s function methods yield the
same conversion rate for ridges. However, the CMS can be applied to more
general domains and is numerically more convenient. The solutions fields of
two methods differ over sloping regions, the solution from the Green’s
function method has to be corrected to properly describe the purely baroclinic
field. This may be of importance for the horizontal phase propagation of
internal tides above topographies, and could have an influence for sea-surface
observations and analysis.Our calculations in the case of a shelf agree  with
the modal solution of Griffiths & Grimshaw (2007). The main advantage of the
present CMS is that it exhibits faster convergence, thus it is numerically
more efficient. By extending our calculations to the strongly supercritical
regime we provide new numerical evidence that the conversion rate approaches
the one obtained in the case of a step (St. Laurent et al., 2003).

The present CMS can be extended in several directions. It can be applied, for
instance, to the investigation of topographically trapped flows for arbitrary
topographies extending the study of Maas & Zimmerman (1989). The assumption of
constant stratification can be easily lifted, at the cost of using a set of
local basis functions that is calculated numerically in the case of smooth
stratification. The linearised free-surface condition can be taken into
account by using basis functions defined in terms of a local transcendental
equation. The solution of the CMS for other topographies such as trenches,
multiple ridges or realistic topographic profiles is straightforward. A Matlab
script implementing the solution of the CMS is freely available through the
following link https://github.com/chpapoutsellis/InternalTidesCMS .

## Declaration of Interests

The authors report no conflict of interest.

## Acknowledgments

The authors acknowledge the support of the IRT/STAE Foundation (funding for
the post-doctoral contract of Ch.P.), and also L.R.M. Maas for fruitful
discussions. We would like to thank the reviewers and J. Nycander for their
comments and suggestions, which greatly improved the quality of this paper.
M.M. thanks Lucile Planes and Martin Benebig for testing computations with
iTides. N.G. acknowledges the support of the Natural Sciences and Engineering
Research Council of Canada (NSERC) [funding reference numbers RGPIN-2015-03684
and RGPIN-2022-04560] and of the Canadian Space Agency [14SUSWOTTO].

##  Appendix A Asymptotic analysis of the barotropic equations

The barotropic flow solves the system (Baines, 1973; Garrett & Gerkema, 2007)

$\displaystyle U_{t}-fV$ $\displaystyle=-P_{x},\quad V_{t}+fU=0,$ (44a)
$\displaystyle W_{t}$ $\displaystyle=-P_{z},$ (44b) $\displaystyle
B_{t}+N^{2}W$ $\displaystyle=0,$ (44c) $\displaystyle U_{x}+W_{z}$
$\displaystyle=0,$ (44d)

where $(U,V,W)$, $P$ and $B$ denote respectively the velocities, the scaled
pressure, and the buoyancy induced by the barotropic flow. On the boundaries,
we have

$h_{x}U(x,-h)+W(x,-h)=0,\quad W(x,0)=0.$ $None$

As described in section 2.2, the above system reduces to the BVP (9) in terms
of $\Phi$. In order to derive an asymptotic solution for (9) we introduce the
scaling,

$\displaystyle\tilde{z}=z/h_{0},\quad\tilde{x}=x/L,\quad\tilde{h}=h/h_{0},\quad\tilde{\Phi}=\Phi/Q,$
(46)

and the parameter $\sigma=h_{0}^{2}/L^{2}$, where $L$ is the horizontal scale
of the topography. The dimensionless version of (9) is written as, after
dropping the tildes,

$\displaystyle\sigma\Phi_{xx}+\mu_{0}^{-2}\Phi_{zz}=0,\quad\Phi(x,0)=0,\quad\Phi\left(x,-h\right)=1.$
(47)

If $\sigma\ll 1$, we can solve (47) by using the asymptotic expansion
$\Phi^{\text{app}}=\sum_{i=0}^{K}\sigma^{i}\Phi^{(i)}$, for some $K\geq 0$.
Substitution of $\Phi^{\text{app}}$ in the first equation of (47) yields

$\displaystyle\sum_{j=0}^{K}\sigma^{j}\left(\Phi^{(j-1)}_{xx}+\mu_{0}^{-2}\Phi^{(j)}_{zz}\right)$
$\displaystyle=O(\sigma^{K+1}).$ (48)

where the convention $\Phi^{(-1)}=0$ is used. By requiring the residual
$O(\sigma^{K+1})$ to be cancelled we obtain the following recurrence relation:

$\displaystyle\begin{aligned}
(j=0):&\quad\Phi^{(0)}_{zz}=0,\quad\Phi^{(0)}(x,0)=0,\quad\Phi^{(0)}(x,-h)=1\\\
(1\leq j\leq K):&\quad\begin{aligned}
\Phi^{(j)}_{zz}&=-\mu_{0}^{2}\Phi^{(j-1)}_{xx},\quad\Phi^{(j)}(x,0)&=0,\quad\Phi^{(j)}(x,-h)&=0.\end{aligned}\end{aligned}$
(49)

Solving the above BVPs we obtain

$\displaystyle\begin{aligned}
\Phi^{(0)}=\frac{z}{-h},\,\,\,\,\Phi^{(j)}=\frac{z\mu_{0}^{2}}{h}\int_{-h}^{0}\int_{-h}^{z^{\prime}}\Phi^{(j-1)}_{xx}\>\mathrm{d}z^{\prime\prime}\>\mathrm{d}z^{\prime}+\mu_{0}^{2}\int^{0}_{z}\int_{-h}^{z^{\prime}}\Phi^{(j-1)}_{xx}\>\mathrm{d}z^{\prime\prime}\>\mathrm{d}z^{\prime}.\end{aligned}$
(50)

Performing the computation for $j=1$ and using (46) we find

$\displaystyle\Phi^{(0)}=Q\frac{z}{-h},\quad\Phi^{(1)}=-Q\mu_{0}^{2}\left(\frac{1}{-h}\right)_{xx}\frac{1}{6}\left(z^{3}-h^{2}z\right).$
(51)

We turn now to the system (44)–($None$). Introducing the scaling
$\tilde{t}=\omega t$ and

$\displaystyle\tilde{U}=U/U_{0},\quad\tilde{V}=V/U_{0},\quad\tilde{W}=W/W_{0},\quad\tilde{P}=P/(\omega
U_{0}L),\quad\tilde{B}=B/g,$ (52)

with $U_{0}=\sqrt{gh_{0}}$ and $W_{0}=U_{0}h_{0}/L$, (44)–($None$) become,
after dropping the tildes,

$\displaystyle U_{t}-\frac{f}{\omega}V$ $\displaystyle=-P_{x},\quad
V_{t}+\frac{f}{\omega}U=0,$ (53a) $\displaystyle\sigma W_{t}$
$\displaystyle=-P_{z},$ (53b) $\displaystyle B_{t}+\frac{N^{2}}{\omega g}W$
$\displaystyle=0,\quad\quad\quad\,\,\,$ (54) $\displaystyle U_{x}+W_{z}$
$\displaystyle=0,\quad\quad\,\,$ (55) $\displaystyle
h_{x}U(x,-h)+W(x,-h)=0,\quad W(x,0)=0.$ (56)

Plugging $\Xi^{\text{app}}=\sum_{i=0}^{K}\sigma^{i}\Xi^{(i)}$, where $\Xi$ is
any of the fields $U$, $V$, $W$, $B$, $P$, into (53)–(56), we see that
$\Xi^{(i)}$ satisfies (53a), (54)–(56) while (53b) gives

$\displaystyle\sum\limits_{i=0}^{K}\sigma^{i}\left(W^{(i-1)}_{t}+P_{z}^{(i)}\right)=O(\sigma^{K+1}).$

Then $O(\sigma^{K+1})$ is cancelled if $W^{(i-1)}_{t}+P_{z}^{(i)}=0$, and in
dimensional form we have

$\displaystyle U^{(i)}_{t}-fV^{(i)}$ $\displaystyle=-P^{(i)}_{x},\quad
V^{(i)}_{t}+fU^{(i)}=0,$ (57a) $\displaystyle W^{(i-1)}_{t}$
$\displaystyle=-P^{(i)}_{z},$ (57b) $\displaystyle B^{(i)}_{t}+N^{2}W^{(i)}$
$\displaystyle=0,\quad\quad\quad\,\,\,$ (58) $\displaystyle
U^{(i)}_{x}+W^{(i)}_{z}$ $\displaystyle=0,\quad\quad\,\,$ (59) $\displaystyle
h_{x}U^{(i)}(x,-h)+W^{(i)}(x,-h)=0,\quad W^{(i)}(x,0)=0.$ (60)

For $i=0$, Eq. (57b) is $P^{(0)}_{z}=0$, thus, at leading order the flow is
hydrostatic.

##  Appendix B Matrix coefficients of the CMS

The matrix coefficients appearing in (32) are given by

$\displaystyle b_{mn}$ $\displaystyle=\left\\{\begin{aligned}
4(-1)^{m+n}\frac{mn}{m^{2}-n^{2}},\,m\neq n\\\ 1,\,m=n\end{aligned}\right.,$
(61) $\displaystyle c_{mn}$ $\displaystyle=\left\\{\begin{aligned}
-4(-1)^{m+n}\frac{mn(m^{2}+n^{2})}{(m^{2}-n^{2})^{2}},\,m\neq n\\\
-\frac{1}{2}-\frac{1}{3}n^{2}\pi^{2},\,m=n\end{aligned}\right.$ (62)
$\displaystyle d_{mn}$ $\displaystyle=\left\\{\begin{aligned}
2(-1)^{m+n}\frac{mn}{m^{2}-n^{2}},\,m\neq n\\\
\frac{1}{2},\,m=n\end{aligned}\right.$ (63)

##  Appendix C Infinitesimal topography solution of the CMS

For a ridge topography, we have that $r\rightarrow 0$ as
$x\rightarrow\pm\infty$ and $h_{0}$ is chosen as the far-field depth. For a
shelf, $h$ is assumed as a depth transition from $h_{0}-\epsilon\,r_{-}$ to
$h_{0}-\epsilon\,r_{+}$, with $r\rightarrow r_{\pm}$ as
$x\rightarrow\pm\infty$ and $r_{+}+r_{-}=0$ and $h_{0}$ represents the mean
depth. Using the Taylor expansion

$\displaystyle\frac{1}{h}=\frac{1}{h_{0}-\epsilon
r}=\frac{1}{h_{0}}+\frac{1}{h_{0}^{2}}\epsilon r+\dots,$

we note that the various $h$-dependent terms in (32) scale as follows

$\displaystyle\frac{h_{x}}{h}$ $\displaystyle=-\epsilon
r_{x}\left(\frac{1}{h_{0}}+\frac{1}{h_{0}^{2}}\epsilon
r+\dots\right)=O(\epsilon),$ $\displaystyle\frac{h_{x}^{2}}{h^{2}}$
$\displaystyle=O(\epsilon^{2})$ $\displaystyle\frac{h_{xx}}{h}$
$\displaystyle=-\epsilon
r_{xx}\left(\frac{1}{h_{0}}+\frac{1}{h_{0}^{2}}\epsilon
r+\dots\right)=O(\epsilon),$ $\displaystyle h\left(\frac{1}{h}\right)_{xx}$
$\displaystyle=\epsilon\frac{r_{xx}}{h_{0}}+O(\epsilon^{2})$

Substituting the above expressions together with
$\phi_{m}=\sum_{i=0}^{K}\epsilon^{i}\phi_{m}^{(i)}$ in (32), we obtain
$\phi_{m}^{(0)}=0$ and (36) for $\phi_{m}^{(1)}$. The solution of (36)
satisfying the radiation conditions (34) with $k_{n}^{\pm}=\ell_{n}$ is

$\displaystyle\phi^{(1)}_{n}=\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\left[A_{n}(x;r)\text{e}^{\text{i}\ell_{n}x}-B_{n}(x;r)\text{e}^{-\text{i}\ell_{n}x}\right],\quad
n\geq 1,$ (64)

with

$\displaystyle
A_{n}(x;r)=\int_{-\infty}^{x}\text{e}^{-\text{i}\ell_{n}s}r_{s}(s)\>\mathrm{d}s,\,\,B_{n}(x;r)=\int_{x}^{+\infty}\text{e}^{\text{i}\ell_{n}s}r_{s}(s)\>\mathrm{d}s;$
(65)

see, e.g., Gerkema & Zimmerman (2008, Chapter 7)). The response field at first
order is reconstructed by (28) with $\phi_{n}$ given by (64) and $Z_{n}$ by
(30) with $h=h_{0}$. The conversion rate is obtained using (64) and (26),

$\displaystyle
C^{\text{WTA}}_{\pm}=\frac{\omega^{2}-N^{2}}{2\omega}\frac{h_{0}}{2}\sum\limits_{n=1}^{\infty}\text{Im}\left\\{\phi^{(1)}_{n}{\color[rgb]{0,0,0}\overline{\phi^{(1)}_{n,x}}}\right\\},\quad
x\rightarrow\pm\infty.$ (66)

In order to compute $C^{\text{WTA}}_{-}$, we first note that as
$x\rightarrow-\infty$, we have $A_{n}\rightarrow 0$ and

$\displaystyle\phi^{(1)}_{n}$
$\displaystyle=-\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}B_{n}(-\infty;r)\,\text{e}^{-\text{i}\ell_{n}x}$
$\displaystyle=-\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\left(\int_{-\infty}^{\infty}\text{e}^{\text{i}\ell_{n}s}r_{s}(s)\>\mathrm{d}s\right)\,\text{e}^{-\text{i}\ell_{n}x}$
$\displaystyle=-\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\left(\left[\text{e}^{\text{i}\ell_{n}s}r(s)\right]_{-\infty}^{+\infty}-\text{i}\ell_{n}\int_{-\infty}^{\infty}\text{e}^{\text{i}\ell_{n}s}r(s)\>\mathrm{d}s\right)\text{e}^{-\text{i}\ell_{n}x}$
$\displaystyle=-\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\left(\left[\text{e}^{\text{i}\ell_{n}s}r(s)\right]_{-\infty}^{+\infty}-\text{i}\ell_{n}\overline{\hat{r}}(\ell_{n})\right)\text{e}^{-\text{i}\ell_{n}x},$

where $\hat{r}(\xi)=\int_{-\infty}^{+\infty}\exp(-\text{i}x\xi)r(s)ds$ is the
Fourier transform of $r$ and $\hat{r}(-\xi)=\overline{\hat{r}}(\xi)$, since
$r\in\mathbb{R}$. Similarly, we find

$\displaystyle\phi^{(1)}_{n,x}$
$\displaystyle=\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\text{i}\ell_{n}B_{n}(-\infty;r)\,\text{e}^{-\text{i}\ell_{n}x}$
$\displaystyle=\frac{g_{n}}{{\color[rgb]{0,0,0}h_{0}}}\left(\text{i}\ell_{n}\left[\text{e}^{\text{i}\ell_{n}s}r(s)\right]_{-\infty}^{+\infty}+\ell_{n}^{2}\overline{\hat{r}}(\ell_{n})\right)\text{e}^{-\text{i}\ell_{n}x},\quad
x\rightarrow-\infty.$

Combining the above results, we compute

$\displaystyle\phi^{(1)}_{n}\overline{\phi^{(1)}_{n,x}}=-\frac{g_{n}^{2}}{{\color[rgb]{0,0,0}h_{0}^{2}}}\Pi_{n},\quad
x\rightarrow-\infty,$

where

$\displaystyle\Pi_{n}$
$\displaystyle=\left(\left[\text{e}^{\text{i}\ell_{n}s}r(s)\right]^{+\infty}_{-\infty}-\text{i}{\color[rgb]{0,0,0}\ell_{n}}\overline{\hat{r}}(\ell_{n})\right)\left(-\text{i}\ell_{n}\left[\text{e}^{-\text{i}\ell_{n}s}r(s)\right]^{+\infty}_{-\infty}+\ell_{n}^{2}\hat{r}(\ell_{n})\right)$
$\displaystyle=-\text{i}\ell_{n}\left[r^{2}\right]^{+\infty}_{-\infty}+\left[\text{e}^{\text{i}\ell_{n}s}r\right]^{+\infty}_{-\infty}\ell_{n}^{2}\hat{r}(\ell_{n})+\ell_{n}^{2}\overline{\hat{r}}(\ell_{n})\left[\text{e}^{-\text{i}\ell_{n}s}r\right]^{+\infty}_{-\infty}-\text{i}\ell_{n}^{3}\overline{\hat{r}}(\ell_{n})\hat{r}(\ell_{n})$
$\displaystyle=-\text{i}\ell_{n}\left[r^{2}\right]^{+\infty}_{-\infty}+\left[\text{e}^{\text{i}\ell_{n}s}r\right]^{+\infty}_{-\infty}\ell_{n}^{2}\hat{r}(\ell_{n})+\overline{\ell_{n}^{2}\hat{r}(\ell_{n})\left[\text{e}^{\text{i}\ell_{n}s}r\right]^{+\infty}_{-\infty}}-\text{i}\ell_{n}^{3}\overline{\hat{r}}(\ell_{n})\hat{r}(\ell_{n})$
$\displaystyle=-\text{i}\ell_{n}\left[r^{2}\right]^{+\infty}_{-\infty}+2\text{Re}\left\\{\left[\text{e}^{\text{i}\ell_{n}s}r\right]^{+\infty}_{-\infty}\ell_{n}^{2}\hat{r}(\ell_{n})\right\\}-\text{i}\ell_{n}^{3}\overline{\hat{r}}(\ell_{n})\hat{r}(\ell_{n})$

The first term of the above right hand side vanishes for both the ridge and
the shelf cases. It then follows from (66) that

$\displaystyle
C^{\text{WTA}}_{-}=-\frac{N^{2}-\omega^{2}}{2\omega}\frac{h_{0}}{2}\sum\limits_{n=1}^{\infty}\frac{g_{n}^{2}}{{\color[rgb]{0,0,0}h_{0}^{2}}}\ell_{n}^{3}\overline{\hat{r}}(\ell_{n})\hat{r}(\ell_{n}).$

Repeating the same procedure for $x\rightarrow+\infty$, we obtain
$C_{+}^{\text{WTA}}=-C_{-}^{\text{WTA}}$. Recalling that
$g_{n}=Q(-1)^{n{\color[rgb]{0,0,0}+1}}/(n\pi)$ and $Q=Uh_{0}$, we easily find
that $\mathcal{C}^{\text{WTA}}=C_{+}^{\text{WTA}}-C_{-}^{\text{WTA}}$ is given
by (37).

##  Appendix D Reconstruction of the free-surface elevation

It is possible to reconstruct the free-surface elevation $\eta$ induced by ITs
within the rigid-lid approximation by using the baroclinic surface pressure at
$z=0$, $p_{s}(x,t)=p^{\\#}(x,0,t)$; $\eta=p_{s}/g$. To find $p_{s}$, we
evaluate $p_{x}^{\\#}$ and $p_{xx}^{\\#}$ on $z=0$, using (16),

$\displaystyle\left[p^{\\#}_{x}\right]_{z=0}=\left[fv^{\\#}-u^{\\#}_{t}\right]_{z=0}:=g(x,t),\,\,\left[p^{\\#}_{xx}\right]_{z=0}=\left[fv^{\\#}_{x}-u^{\\#}_{xt}\right]_{z=0}:=F(x,t),$
(67)

and we formulate and solve the following BVP on $[x_{L},x_{R}]$,

$\displaystyle p^{\\#}_{s,xx}=F(x,t),\quad
p^{\\#}_{s,x}(x_{L})=g(x_{L},t),\quad p^{\\#}_{s,x}(x_{R})=g(x_{R},t).$ (68)

## References

  * (1)
  * Athanassoulis & Belibassakis (1999) Athanassoulis, G. A. & Belibassakis, K. A. (1999), ‘A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions’, J. Fluid Mech. 389, 275–301.
  * Athanassoulis & Papoutsellis (2017) Athanassoulis, G. A. & Papoutsellis, Ch. E. (2017), ‘Exact semi-separation of variables in waveguides with non-planar boundaries’, Proc. R. Soc. A. 473(2201), 20170017.
  * Baines (1973) Baines, P. (1973), ‘The generation of internal tides by flat-bump topography’, Deep-Sea Res. 20(2), 179 – 205.
  * Baines (1982) Baines, P. (1982), ‘On internal tide generation models’, Deep Sea Res. Part I Oceanogr. Res. Pap. 29(3), 307 – 338.
  * Balmforth et al. (2002) Balmforth, N. J., Ierley, G. R. & Young, W. R. (2002), ‘Tidal conversion by subcritical topography’, J. Phys. Oceanogr. 32(10), 2900–2914.
  * Balmforth & Peacock (2009) Balmforth, N. J. & Peacock, T. (2009), ‘Tidal conversion by supercritical topography’, J. Phys. Oceanogr. 39(8), 1965–1974.
  * Bell (1975) Bell, T. H. (1975), ‘Lee waves in stratified flows with simple harmonic time dependence’, J. Fluid Mech. 67, 705–722.
  * Brekhovskikh & Godin (1992) Brekhovskikh, L. M. & Godin, O. (1992), Acoustics of Layered Media II: point sources and bounded beams, Springer.
  * Bühler & Muller (2007) Bühler, O. & Muller, C. J. (2007), ‘Instability and focusing of internal tides in the deep ocean’, J. Fluid Mech. 588, 1–28.
  * Craig (1987) Craig, P. (1987), ‘Solutions for internal tidal generation over coastal topography’, J. Mar. Res. 45, 83–105.
  * Dauxois et al. (2018) Dauxois, T., Joubaud, S., Odier, P. & Venaille, A. (2018), ‘Instabilities of Internal Gravity Wave Beams’, Annu. Rev. Fluid Mech. 50(1), 131–156.
  * Desaubies & Dysthe (1995) Desaubies, Y. & Dysthe, K. (1995), ‘Normal-mode propagation in slowly varying ocean waveguides’, J. Acoust. Soc 97(2), 933–946.
  * Echeverri & Peacock (2010) Echeverri, P. & Peacock, T. (2010), ‘Internal tide generation by arbitrary two-dimensional topography’, J. Fluid Mech. 659, 247–266.
  * Echeverri et al. (2011) Echeverri, P., Yokossi, T., Balmforth, N. J. & Peacock, T. (2011), ‘Tidally generated internal-wave attractors between double ridges’, J. Fluid Mech. 669, 354–374.
  * Falahat et al. (2014) Falahat, S., Nycander, J., Roquet, F. & Zarroug, M. (2014), ‘Global calculation of tidal energy conversion into vertical normal modes’, J. Phys. Oceanogr. 44(12), 3225 – 3244.
  * Garrett (2003) Garrett, C. (2003), ‘Internal tides and ocean mixing’, Science 301(5641), 1858–1859.
  * Garrett & Gerkema (2007) Garrett, C. & Gerkema, T. (2007), ‘On the Body-Force Term in Internal-Tide Generation’, J. Phys. Oceanogr. 37(8), 2172–2175.
  * Garrett & Kunze (2007) Garrett, C. & Kunze, E. (2007), ‘Internal tide generation in the deep ocean’, Annu. Rev. Fluid Mech. 39(1), 57–87.
  * Gerkema et al. (2004) Gerkema, T., Lam, F.-P. A. & Maas, L. R. (2004), ‘Internal tides in the Bay of Biscay: conversion rates and seasonal effects’, Deep Sea Res. Part II Top. Stud. Oceanogr. 51(25), 2995 – 3008.
  * Gerkema & Zimmerman (2008) Gerkema, T. & Zimmerman, J. (2008), An introduction to internal waves: Lecture Notes, R. Neth. Inst. for Sea Res., Den Burg.  https://www.vliz.be/imisdocs/publications/ocrd/60/307760.pdf.
  * Griffiths & Grimshaw (2007) 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(3), 428–451.
  * Grisouard & Bühler (2012) Grisouard, N. & Bühler, O. (2012), ‘Forcing of oceanic mean flows by dissipating internal tides’, J. Fluid Mech. 708, 250–278.
  * Groetsch (1990) Groetsch, C. (1990), ‘Convergence analysis of a regularized degenerate kernel method for fredholm integral equations of the first kind.’, Integr. Equ. Oper. theory 13, 67 – 75.
  * Groetsch (1984) Groetsch, C. W. (1984), The theory of Tikhonov regularization for Fredholm equations of the first kind, Pitman Advanced Publishing Program.
  * He et al. (2019) He, J., Homa, D., Pickrell, G. & Wang, A. (2019), ‘Coupled mode analysis for 3d stress-free elastic acoustic waveguide’, IEEE Access 7, 117796–117803.
  * Ivansson (2021) Ivansson, S. M. (2021), ‘Coupled-mode field computations for media with locally reacting irregular boundaries’, J. Acoust. Soc 150(4), 2985–2998.
  * Kelly (2016) 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(12), 3777 – 3788\.
  * Kelly & Lermusiaux (2016) Kelly, S. M. & Lermusiaux, P. F. J. (2016), ‘Internal-tide interactions with the Gulf Stream and Middle Atlantic Bight shelf break front’, J. Geophys. Res. Oceans 121(8), 6271–6294.
  * Khatiwala (2003) Khatiwala, S. (2003), ‘Generation of internal tides in an ocean of finite depth: analytical and numerical calculations’, Deep Sea Res. Part I Oceanogr. Res. Pap. 50(1), 3 – 21.
  * Klymak et al. (2013) Klymak, J. M., Buijsman, M., Legg, S. & Pinkel, R. (2013), ‘Parameterizing surface and internal tide scattering and breaking on supercritical topography: The one- and two-ridge cases’, J. Phys. Oceanogr. 43(7), 1380 – 1397.
  * Lahaye & Llewellyn Smith (2020) Lahaye, N. & Llewellyn Smith, S. G. (2020), ‘Modal Analysis of Internal Wave Propagation and Scattering over Large-Amplitude Topography’, J. Phys. Oceanogr. 50(2), 305–321.
  * Larsen (1969) Larsen, L. (1969), ‘Internal waves incident upon a knife edge barrier’, Deep Sea Res. Part I Oceanogr. Res. Pap. 16(5), 411–419.
  * Llewellyn Smith & Young (2002) Llewellyn Smith, S. G. & Young, W. R. (2002), ‘Conversion of the barotropic tide’, J. Phys. Oceanogr. 32 (5), 1554–1566.
  * Llewellyn Smith & Young (2003) Llewellyn Smith, S. G. & Young, W. R. (2003), ‘Tidal conversion at a very steep ridge’, J. Fluid Mech. 495, 175–191.
  * Maas (2011) Maas, L. R. M. (2011), ‘Topographies lacking tidal conversion’, J. Fluid Mech. 684, 5–24.
  * Maas & Zimmerman (1989) Maas, L. R. M. & Zimmerman, J. T. F. (1989), ‘Tide-topography interactions in a stratified shelf sea II. Bottom trapped internal tides and baroclinic residual currents’, Geophys. Astrophys. Fluid Dyn. 45, 37–69.
  * MacKinnon et al. (2017) MacKinnon, J. A., Zhao, Z., Whalen, C. B., Waterhouse, A. F., Trossman, D. S., Sun, O. M., St. Laurent, L. C., Simmons, H. L., Polzin, K., Pinkel, R., Pickering, A., Norton, N. J., Nash, J. D., Musgrave, R., Merchant, L. M., Melet, A. V., Mater, B., Legg, S., Large, W. G., Kunze, E., Klymak, J. M., Jochum, M., Jayne, S. R., Hallberg, R. W., Griffies, S. M., Diggs, S., Danabasoglu, G., Chassignet, E. P., Buijsman, M. C., Bryan, F. O., Briegleb, B. P., Barna, A., Arbic, B. K., Ansong, J. K. & Alford, M. H. (2017), ‘Climate process team on internal wave–driven ocean mixing’, Bull. Am. Meteorol. Soc. 98(11), 2429–2454.
  * Mathur et al. (2016) 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(7), 2157–2168.
  * Maupin (1988) Maupin, V. (1988), ‘Surface waves across 2-D structures: a method based on coupled local modes’, Geophys. J. Int. 93(1), 173–185.
  * Maurel et al. (2014) Maurel, A., Mercier, J.-F. & Félix, S. (2014), ‘Propagation in waveguides with varying cross section and curvature: a new light on the role of supplementary modes in multi-modal methods’, Proc. Math. Phys. Eng. 470(2166), 20140008.
  * Mercier et al. (2012) Mercier, M. J., Ghaemsaidi, S. J., Echeverri, P., Mathur, M. & Peacock, T. (2012), ‘iTides’, Zenodo doi:10.5281/zenodo.4421548.  https://zenodo.org/record/4421548#.Yx2kxHZBxPY.
  * Muller & Bühler (2009) Muller, C. J. & Bühler, O. (2009), ‘Saturation of the Internal Tides and Induced Mixing in the Abyssal Ocean’, J. Phys. Oceanogr. 39(9), 2077–2096.
  * New (1988) New, A. (1988), ‘Internal tidal mixing in the bay of biscay’, Deep Sea Res. Part I Oceanogr. Res. Pap. 35(5), 691–709.
  * Nycander (2005) Nycander, J. (2005), ‘Generation of internal waves in the deep ocean by tides’, J. Geophys. Res. Oceans 110, C10028.
  * Nycander (2006) Nycander, J. (2006), ‘Tidal generation of internal waves from a periodic array of steep ridges’, J. Fluid Mech. 567, 415 – 432.
  * Pagneux & Maurel (2006) Pagneux, V. & Maurel, A. (2006), ‘Lamb wave propagation in elastic waveguides with variable thickness’, Proc. Math. Phys. Eng. 462(2068), 1315–1339.
  * Papoutsellis et al. (2018) Papoutsellis, Ch. E., Charalampopoulos, A. & Athanassoulis, G. A. (2018), ‘Implementation of a fully nonlinear Hamiltonian Coupled-Mode Theory, and application to solitary wave problems over bathymetry’, Eur. J. Mech. B Fluids 72, 199–224.
  * Papoutsellis et al. (2019) Papoutsellis, Ch. E., Yates, M. L., Simon, B. & Benoit, M. (2019), ‘Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model’, Coast. Eng. 154, 103579.
  * Porter & Staziker (1995) Porter, D. & Staziker, D. J. (1995), ‘Extensions of the mild-slope equation’, J. Fluid Mech. 300, 367–382.
  * Prinsenberg et al. (1974) Prinsenberg, S., Wilmot, W. & Rattray, M. (1974), ‘Generation and dissipation of coastal internal tides’, Deep Sea Res. Part I Oceanogr. Res. Pap. 21(4), 263–281.
  * Pétrélis et al. (2006) Pétrélis, F., Llewellyn Smith, S. & Young, W. R. (2006), ‘Tidal conversion at a submarine ridge’, J. Phys. Oceanogr. 36(6), 1053–1071.
  * Raisbeck (1955) Raisbeck, G. (1955), ‘The order of magnitude of the fourier coefficients in functions having isolated singularities.’, Am. Math. Mon. 62(3), 149–154.
  * Robinson (1969) Robinson, R. (1969), ‘The effects of a vertical barrier on internal waves’, Deep Sea Res. Part I Oceanogr. Res. Pap. 16(5), 421 – 429.
  * Salem (1939) Salem, R. (1939), Essais sur les séries trigonométriques, number 230 in ‘Thèses de l’entre-deux-guerres’.  http://www.numdam.org/item/THESE_1939__230__1_0/.
  * Sandstrom (1975) Sandstrom, H. (1975), ‘On topographic generation and coupling of internal waves’, Geophys. fluid. dyn. 7(1), 231–270.
  * St. Laurent et al. (2003) St. Laurent, L., Stringer, S., Garrett, C. & Perrault-Joncas, D. (2003), ‘The generation of internal tides at abrupt topography’, Deep Sea Res. Part I Oceanogr. Res. Pap. 50(8), 987 – 1003.
  * Stashchuk & Cherkesov (1991) Stashchuk, N. & Cherkesov, L. (1991), ‘Generation of internal waves resulting from the interaction of a barotropic tide with a horizontally- inhomogeneous density held and bottom topography’, Sov. J. Phys. Oceanogr., 2, 79–87.
  * Thyng et al. (2016) Thyng, K. M., Greene, C. A., Hetland, R. D., Zimmerle, H. M. & DiMarco, S. F. (2016), ‘True colors of oceanography: Guidelines for effective and accurate colormap selection’, Oceanography .
  * Vlasenko et al. (2005) Vlasenko, V., Stashchuk, N. & Hutter, K. (2005), Baroclinic Tides, Cambridge University Press.
  * Whalen et al. (2020) Whalen, C. B., de Lavergne, C., Garabato, A. C. N., 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.
  * Wunsch & Ferrari (2004) Wunsch, C. & Ferrari, R. (2004), ‘Vertical mixing, energy, and the general circulation of the oceans’, Annu. Rev. Fluid Mech. 36(1), 281–314.
  * Wunsch & Wunsch (2022) Wunsch, C. & Wunsch, J. (2022), ‘Baroclinic tidal conversion: note on a paper of L.R.M. Maas’, J. Fluid Mech. 946, A47.