Internal tide generation from non-uniform barotropic body forcing

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 (respectively shelf) approaches the value obtained by the knife-edge (respectively 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 parametrisations involving idealised estimations of the local energy conversion rate at the IT generation sites. Moreover, the dissipation of ITs is also parametrised 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 the 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, Stashchuk & Hutter 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; Laurent et al. 2003; Llewellyn Smith & Young 2003; Nycander 2006), or top-hat, linear ramp and step functions (Prinsenberg, Wilmot & Rattray 1974; New 1988; 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) and Falahat et al. (2014).

The first solution to the full two-dimensional (2-D) 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, Lam & Maas (2004) obtained numerical solutions of the formulation of Baines (1973) in the time domain by treating the radiation conditions with sponge layers. Garrett & Gerkema (2007) showed that the body-forcing term in the formulation of Baines (1973) 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, 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, Llewellyn Smith & Young 2006; Echeverri & Peacock 2010; Mathur, Carter & Peacock 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 § 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 the body-forcing formulation of Garrett & Gerkema (2007). 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, Mercier & Félix 2014; Ivansson 2021), elasticity (Maupin 1988; Pagneux & Maurel 2006; He et al. 2019) and water waves (Porter & Staziker 1995; Athanassoulis & Belibassakis 1999; Papoutsellis, Charalampopoulos & Athanassoulis 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 2-D ITs over a shelf topography. Similar systems are derived by 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 § 2, we present the body-forcing formulation and derive its energy balance equation. In § 3, we introduce the modal decomposition of the stream function and derive our CMS. In § 4, we present numerical convergence and accuracy results, and in § 5, we visualise our solutions. In § 6, we consider the conversion rates for various topographies. In § 7, we discuss the differences of the proposed CMS with the Green's function method and, finally, in § 8, we present our conclusions.

2. Governing equations

2.1. Posing the problem

In a 2-D 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 $\varLambda$ and the characteristic horizontal scaling length $L$. We do not limit ourselves to small heights; i.e. $\varLambda$ 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 _{eq}(z)$ weakly departs from a constant reference density $\rho _0$ so that the Boussinesq approximation applies. We further assume that $\rho _{eq}(z)$ decreases linearly with $z$ such that the Brunt–Väisälä frequency $N = \sqrt {(-g/\rho _0)\,\text {d}\rho _{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= -\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 volume flux of constant amplitude $Q$ corresponding to a uniform current $Q/h_{\pm }\cos (\omega t)$ at $x\rightarrow \pm \infty$.

Under the Boussinesq approximation, our set-up is characterised by eight dimensional parameters ($h_{\pm }$, $L$, $\varLambda$, $f$, $N$, $\omega$ and $Q$), which we summarise in figure 1, measured in combinations of metres and seconds. In the case of an isolated ridge, $h_0$ is the depth as $x\rightarrow \pm \infty$, and $\varLambda$ is the maximum height of the ridge, $\varLambda =\max \{h_0-h\}$. For a shelf, we assume that $h_- > h_+$ ($h_-$ is the oceanward depth) so that the maximum height is $\varLambda =\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,

(2.1)\begin{equation} \delta = \frac{\max\{h_0-h\}}{h_0} \sim \frac\varLambda{h_0}. \end{equation}

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

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

$\theta$ being the angle of the free internal wave group velocity with respect to the horizontal plane (see (2.8a–c)). 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$,

(2.3)\begin{equation} \varepsilon = \mu\max\{|\partial_x h|\} \sim \frac{\mu\varLambda}L, \end{equation}

and the tidal excursion $\tau$ defined by

(2.4)\begin{equation} \tau = \frac{Q}{(h_0 - \varLambda)\omega L}. \end{equation}

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-\varLambda )\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’.

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.

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 $^\circ$N, $\omega /f = 1.4$ for the M$_2$ component and we will use $N = 1.5\times 10^{-3}$ s$^{-1}$, which implies $\mu \approx 15$). Starting in § 4, $\varepsilon$ and $\delta$ are the parameters that we vary primarily. Finally, we adopt a standard fixed value for $Q = 120$ m$^2$ 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, Ierley & Young 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 (Bühler & Muller 2007; Garrett & Kunze 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 nonlinear 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 parametrised 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 /\rho _0 - 1)$, where $\rho$ is the total density field, the governing equations are

(2.5a)\begin{gather} u_t - f v ={-}p_x,\quad v_t + f u = 0, \end{gather}

(2.5b)\begin{gather}w_t ={-}p_z + b, \end{gather}

(2.5c)\begin{gather}b_t + N^2 w = 0, \end{gather}

(2.5d)\begin{gather}u_x + w_z = 0, \end{gather}

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

(2.6a,b)\begin{equation} h_x u(x,-h) + w(x,-h) = 0\quad\text{and}\quad w(x,0) = 0,\end{equation}

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

(2.7)\begin{equation} \int_{{-}h}^0u \,\,\mathrm{d} z = Q \cos\omega t. \end{equation}

Introducing the time harmonic stream function $\psi = \text {Re}\{\phi (x,z)\exp ({-\text {i}\omega t})\}$, with $u = -\psi _z$, $w = \psi _x$, where $\phi$ is a complex amplitude and $\text {Re}$ stands for the real part, (2.5)–(2.7) lead to a hyperbolic boundary value problem (BVP) on $\phi$,

(2.8a–c)\begin{equation} \mathcal{L}_{\mu}\phi := \left(\partial_x^2 - \mu^{{-}2}\partial_z^2\right) \phi= 0,\quad \phi(x,0) =0\quad\text{and}\quad \phi(x,-h) = Q, \end{equation}

where $\mu$ is defined in (2.2) (Vlasenko et al. 2005; Garrett & Gerkema 2007). Equation (2.8a–c), 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 (2.5)–(2.6a,b), where the buoyancy in the vertical momentum equation (2.5b) is absent and (2.5c) acts as a diagnostic equation for the induced buoyancy $B$, see (A1) in Appendix A. Introducing a barotropic stream function $\varPsi = \text {Re}\{\varPhi \exp ({-\text {i} \omega t})\}$ for some time independent function $\varPhi$ with $U = -\varPsi _z$ and $W = \varPsi _x$, we obtain from (A1),

(2.9a–c)\begin{equation} \varPhi_{xx}+ \mu_0^{{-}2}\varPhi_{zz}=0,\quad \varPhi(x,0)=0\quad\text{and}\quad \varPhi(x,-h)=Q, \end{equation}

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 (2.9a) is elliptic and is also obtained by (2.8a) with $N\equiv 0$. The solution of (2.9a–c) can be written as

(2.10)\begin{equation} \varPhi = \varPhi^{(0)}+\varPhi^{r},\quad\text{with}\ \varPhi^{(0)} ={-}Qz/h, \end{equation}

where $\varPhi ^{(0)}$ represents the hydrostatic part of the barotropic flow. It is obtained by neglecting $W_t$ in (A1b) and, consequently, $\varPhi _{xx}$ in (2.9a); see Appendix A for a detailed derivation. Here, $\varPhi ^{r}$ represents the residual non-hydrostatic part which solves

(2.11a–c)\begin{equation} \varPhi^{r}_{xx}+ \mu_0^{{-}2}\varPhi^{r}_{zz}={-}\varPhi^{(0)}_{xx},\quad \varPhi^{r}(x,0)=0\quad\text{and}\quad \varPhi^{r}(x,-h)=0. \end{equation}

Note that $\varPhi$ is spatially non-uniform, that is, it depends on $z$ and $h(x)$. Also note that $\varPhi ^{r}$ vanishes at $x\rightarrow \pm \infty$. We obtain a general semi-analytical solution of (2.11a–c) that is valid for arbitrary smooth $h$ by means of a modal decomposition (§ 3.2). In Appendix A, we derive a perturbative solution that is valid if $h_0/L\ll 1$. In fact, $\varPhi ^{(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. In the IT model of Griffiths & Grimshaw (2007), $U^{(0)}$ is also used as a forcing in the case of a shelf without coastline.

2.3. The internal tide generation problem

Introducing $\phi ^{\#}=\phi -\varPhi$, we obtain from (2.8a–c) the BVP

(2.12a–c)\begin{equation} \mathcal{L}_{\mu} \phi^{\#} ={-}\mathcal{L}_{\mu}\varPhi,\quad \phi^{\#}(x,0)=0\quad\text{and}\quad \phi^{\#}(x,-h)= 0, \end{equation}

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

(2.13a–c)\begin{equation} \mathcal{L}_{\mu} \phi^{{{\dagger}}} ={-}\mathcal{L}_{\mu}\varPhi^{(0)} ={-}\varPhi^{(0)}_{xx},\quad \phi^{{{\dagger}}}(x,0)=0\quad\text{and}\quad \phi^{{{\dagger}}}(x,-h)= 0, \end{equation}

which shows that the hydrostatic part of barotropic flow, $\varPhi ^{(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 (§ 3). Here, we shall proceed with (2.13a–c) mainly because $\varPhi ^{(0)}$ is given by the simple explicit expression (2.10); $\varPhi ^{\text r}$ can be computed independently to extract the purely baroclinic field $\phi ^{\#}=\phi ^{{\dagger} }-\varPhi ^{r}$ if needed. Also, $\varPhi ^{\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. Equation (2.13a–c) is supplemented with radiation conditions ensuring that waves generated in the interior of the domain propagate outward as plane waves. Thus, we have

(2.14)\begin{equation} \phi^{{{\dagger}}} = \sum_{n=1}^{\infty}c^\pm_n\exp({\pm \text{i} k^{{\pm}}_n x})\sin \left(\frac{n{\rm \pi} z}{h_{{\pm}}}\right),\quad \text{with}\ k^{{\pm}}_n = \frac{n{\rm \pi}}{\mu h_{{\pm}}},\enspace c^{{\pm}}_n\in\mathbb{C},\quad \text{as}\ x\rightarrow\pm\infty. \end{equation}

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 = \varPhi ^{(0)}+\varPhi ^{r}+\phi ^{\#} =\varPhi ^{(0)} + \phi ^{{\dagger} }$ and we may introduce the corresponding definitions

(2.15)\begin{equation} \xi = \varXi^{(0)} + \varXi^{\text r} + \xi^\# = \varXi^{(0)} + \xi^{\dagger}, \end{equation}

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

(2.16)\begin{equation} \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} -\varPhi^{{\diamond}}_z & \dfrac{\text{i} f}\omega\varPhi^{{\diamond}}_z & \varPhi^{{\diamond}}_x & -\dfrac{\text{i} N^2}{\omega}\varPhi^{{\diamond}}_x\\ -\phi^\star_z & \dfrac{\text{i} f}\omega\phi_z^\star & \phi_x^\star & -\dfrac{\text{i} N^2}{\omega}\phi^\star_x \end{pmatrix}\exp({-\text{i} \omega t})\right\}, \end{equation}

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 }$ (respectively $b^{\star },B^{\diamond }$) follow from the second equation in (2.5a) (respectively (2.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 notation for the different flow fields in table 1. Plugging the second equality of (2.15) into (2.5)–(2.6a,b) and taking into account (A15)–(A18a,b) with $i=0$, we obtain

(2.17a)\begin{gather} u^{{{\dagger}}}_t - f v^{{{\dagger}}} ={-}p^{{{\dagger}}}_x, \quad v^{{{\dagger}}}_t + f u^{{{\dagger}}} = 0, \end{gather}

(2.17b)\begin{gather}w^{{{\dagger}}}_t ={-}p^{{{\dagger}}}_z + b^{{{\dagger}}} + \left(1-\frac{\omega^2}{N^2}\right)B^{(0)}, \end{gather}

(2.18)\begin{gather} b^{{{\dagger}}}_t + N^2 w^{{{\dagger}}} = 0, \end{gather}

(2.19)\begin{gather}u^{{{\dagger}}}_x + w^{{{\dagger}}}_z = 0, \end{gather}

(2.20a,b)\begin{gather}h_x u^{{{\dagger}}}(x,-h)+w^{{{\dagger}}}(x,-h) = 0,\quad w^{{{\dagger}}}(x,0) = 0. \end{gather}

In deriving (2.17b), we used the relation $B^{(0)} = N^2 W^{(0)}_t/ \omega ^2$, which itself is derived from (2.16). In terms of $\xi ^{\#}$, the above equations stay the same except for (2.17b), which becomes

(2.21)\begin{equation} w^{\#}_t ={-}p^{\#}_z + b^{\#} + B, \end{equation}

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

Table 1. The notation for the fields we use. The symbol $\xi$ stands for $u$, $v$, $w$, $b$, $p$ and $\varXi$ for their capitalized versions. The relations in the third column hold for $\xi$ and $\varXi$ replaced by $\phi$ and $\varPhi$.

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 (2.17)–(2.18) with $(\boldsymbol {u}^{{\dagger} },b^{{\dagger} }/N^2) \equiv (u^{{\dagger} }, v^{{\dagger} }, w^{{\dagger} },b^{{\dagger} }/N^2)$ gives

(2.22)\begin{equation} \mathcal{E}^{\dagger}_t + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(p^{\dagger}\boldsymbol{u}^{\dagger}\right) = \left(1-\frac{\omega^2}{N^2}\right)B^{(0)} w^{\dagger},\quad \text{with}\ \mathcal{E}^{\dagger}= \frac{1}{2}(\boldsymbol{u}^{{{\dagger}}})^2 +\frac{1}{2}\frac{(b^{{{\dagger}}})^2}{N^2},\end{equation}

where we have used (2.19) and where $\boldsymbol {\nabla } \equiv (\partial _x, 0,\partial _z )$. Integrating (2.22) over the domain $\varOmega = [-\infty,+\infty ]\times [-h,0]$, using the divergence theorem and (2.20a,b), we obtain

(2.23)\begin{equation} \left(\int_\varOmega \mathcal{E}^{\dagger} \,\mathrm{d} \varOmega\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_{\varOmega}B^{(0)} w^{{{\dagger}}} \,\,\mathrm{d} \varOmega, \end{equation}

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^{\rm T} {\cdot } \,\mathrm {d} t$, with $T = 2{\rm \pi} /\omega$, and taking into account the periodicity of $\mathcal {E}^{{\dagger} }(t)$, we find

(2.24)\begin{equation} C_+-C_{-}\stackrel{\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_{\varOmega} \langle B^{(0)} w^{{{\dagger}}}\rangle \,\mathrm{d} \varOmega \stackrel{\text{def}}=C_{int}, \end{equation}

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_{int}$ given as a volume integral. Note that the non-hydrostatic barotropic flow does not contribute in (2.24) because $u^{{\dagger} }=u^{\#}+U^r \rightarrow u^{\#}$ as $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 ${\rm \pi} /2$ (2.16). Thus, (2.24) remains valid with ${\dagger}$ replaced by $\#$.

We proceed by expressing $C_{\pm }$ and $C_{int}$ in terms of $\phi ^{{\dagger} }$ and $\varPhi ^{(0)}$. Using integration by parts and the fact that $\psi ^{{\dagger} }(x,0) =\psi ^{{\dagger} }(x,-h)=0$, we obtain $\int _{-h}^0 \langle p^{{\dagger} } u^{{\dagger} }\rangle {\rm d} z = \int _{-h}^0 \langle p_z^{{\dagger} } \psi ^{{\dagger} }\rangle {\rm d} z$. Then, expressing $p_z^{{\dagger} }$ from (2.17b) using (2.16), we find

(2.25)\begin{equation} \langle p_z^{{{\dagger}}}\psi^{{{\dagger}}}\rangle = \frac{\omega^2-N^2}{\omega} \langle\text{Re}\{\text{i}\phi^{{{\dagger}}}_x \exp({-\text{i}\omega t})\}\text{Re}\{\phi^{{{\dagger}}} \exp({-\text{i}\omega t})\}\rangle\quad \text{as}\ x\rightarrow\pm\infty. \end{equation}

Writing $\exp ({-\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

(2.26)\begin{equation} C_{{\pm}} = \frac{\omega^2-N^2}{2\omega}\int_{{-}h}^0 \text{Im}\{\phi^{{{\dagger}}} \overline{\phi_x^{{{\dagger}}}}\}\,\,\mathrm{d} z,\quad \text{as}\ x\rightarrow\pm\infty, \end{equation}

where the overline denotes the complex conjugate and $\text {Im}$ the imaginary part. Similarly,

(2.27)\begin{equation} C_{int} = \left(1-\frac{\omega^2}{N^2}\right)\frac{N^2}{2\omega} \int_{\varOmega} \varPhi^{(0)}_x\text{Im}\{\overline{\phi^{{{\dagger}}}_x}\}\,\,\mathrm{d} \varOmega. \end{equation}

Using the radiation conditions (2.14) in (2.26), we see that $C_+\geq 0$ (respectively $C_-\leq 0$) as $x\rightarrow +\infty$ (respectively $x\rightarrow -\infty$). Equations (2.26) and (2.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 § 4 for validation purposes.

3. Modal decomposition

3.1. Stream function modal representation

We reformulate the IT generation problem (2.13a–c)–(2.14) by representing $\phi ^{{\dagger} }$ as

(3.1)\begin{equation} \phi^{{{\dagger}}}(x,z) = \sum_{n=1}^{\infty}\phi_n(x)Z_n(z;x),\end{equation}

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 (3.1) 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 parametrised by $x$, also called the ‘reference waveguide’ (Brekhovskikh & Godin 1992),

(3.2a–c)\begin{equation} Z_{n,zz} + \frac{n^2{\rm \pi}^2}{h^2} Z_n = 0,\quad Z_n(0 ;x) = 0,\quad Z_n({-}h ;x)= 0. \end{equation}

The eigenfunctions $Z_n$ are given by

(3.3)\begin{equation} Z_n(z;x) = \sin\left(\frac{n{\rm \pi} z}{h}\right),\end{equation}

and satisfy the orthogonality relation $\int _{-h}^{0}Z_n Z_m \,{\rm d} z = h\delta _{nm}/2$, where $\delta _{nm}$ is the Kronecker delta. Note that (3.1) satisfies exactly and term by term the boundary conditions in (2.13a–c). It follows from (3.1) and (3.3) that the $\phi _n$ are defined by

(3.4)\begin{equation} \phi_n = \frac{2}{h}\int_{{-}h}^0\phi^{{{\dagger}}} Z_n \,{\rm d} z = \frac{2h^2}{{\rm \pi}^3 n^3}\left( \left[\phi^{{{\dagger}}}_{zz}Y_n\right]_{{-}h}^0 - \int_{{-}h}^0 \phi^{{{\dagger}}}_{z z z}Y_n \,\,\mathrm{d} z\right), \end{equation}

where $Y_n=\cos (n{\rm \pi} z / h)$ and the second equality is obtained after integrating by parts three times, which is allowed provided $\phi ^{{\dagger} }$ is sufficiently smooth, and using the boundary conditions in (2.13a–c). This shows that $\|\phi _n\|_{\infty }:=\max |\phi _n|=O(n^{-3})$ and that (3.1) 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, § 4) and suffice to establish the term-wise differentiability of (3.1) required for the exact modal reformulation of (2.13a–c).

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 (3.1), the truncated version of this expansion does not satisfy term by term the bottom boundary condition. 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 by Kelly & Lermusiaux (2016) for more details. Despite (3.1) being limited to 2-D flows, its major advantage is that it satisfies the bottom boundary condition in (2.13c) exactly and term by term, and exhibits faster convergence in comparison with existing approaches.

3.2. The coupled-mode system

We proceed by projecting (2.13a–c)–(2.14) onto (3.3). Substituting (3.1) into (2.13a), multiplying with $Z_n$ and integrating over the interval $[-h(x),0]$, we find that $\{\phi _n(x)\}_{n=1}^{\infty }$ solves the following CMS, for $m\geq 1$:

(3.5)\begin{equation} \phi_{m,xx}+\frac{m^2{\rm \pi}^2}{\mu^2h^2}\phi_m+\sum_{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]= 2 g_m h\left(\frac{1}{h}\right)_{xx}, \end{equation}

where $b_{mn}$, $c_{mn}$, $d_{mn}$ are $x$-independent coefficients given in Appendix B, and $g_m = Q (-1)^{m+1}/(m{\rm \pi} )$. Substituting (3.1) into (2.14), we obtain

(3.6)\begin{equation} \sum_{n=1}^{\infty}\phi_n(x)Z_n(x;z) = \sum_{n=1}^{\infty}c_n^\pm\exp({\pm} \text{i} k^{{\pm}}_n x)\sin \left(\frac{n{\rm \pi} z}{h_{{\pm}}}\right)\quad\text{as}\ x\rightarrow \pm\infty. \end{equation}

Multiplying both sides by $Z_q$, $q\geq 1$, integrating over $[-h,0]$ and taking into account that $Z_n\rightarrow \sin (n{\rm \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

(3.7)\begin{equation} \phi_{n,x} \mp\text{i} k^{{\pm}}_n\phi_n= 0,\quad\text{as}\ x\rightarrow {\pm}\infty. \end{equation}

Thus, (2.13a–c)–(2.14) are exactly reformulated as the CMS (3.5)–(3.7). Solving the latter, we may reconstruct the solution of the former by using (3.1). The energy conversion rates $C_{\pm }$ are evaluated by (2.26):

(3.8)\begin{equation} C_{{\pm}} = \frac{\omega^2-N^2}{2\omega}\frac{h}{2} \sum_{n=1}^{\infty}\text{Im} \{\phi_n^{{{\dagger}}}\overline{\phi^{{{\dagger}}}}_{\!\!\!n,x}\} ={\pm} \frac{N^2-\omega^2}{4\omega\mu}{\rm \pi} \sum_{n=1}^{\infty}n \phi_n^{{{\dagger}}}\overline{\phi^{{{\dagger}}}}_{\!\!\!n}\quad\text{as}\ x\rightarrow {\pm}\infty, \end{equation}

where (3.7) and the definition of $k_n^{\pm }$ in (2.14) have been used to obtain the second equality. The purely baroclinic field $\phi ^{\#}$ is computed using $\phi ^{\#} = \phi ^{{\dagger} }-\varPhi ^{r}$. For the computation of $\varPhi ^{r}$, we apply the same modal decomposition to (2.11a–c). The resulting CMS is given by (3.5) with $\mu ^2$ replaced by $-\mu _0^2$ and vanishing conditions at infinity, instead of (3.7); if $h_0/L\ll 1$, one could use instead the approximate asymptotic expression in (A8a,b). In the following subsection, we derive a perturbative solution of the CMS (3.5)–(3.7) for infinitesimal topography. For arbitrary topographies, we solve the CMS numerically (see § 3.4).

Remark 1 The body-forcing formulation (2.12a–c) and its coupled-mode reformulation (3.5) 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$), (2.12a–c) becomes (2.13a–c) with $\mu ^2$ replaced $\mu '^2=(\omega ^2-f^2)/N^2$ and $\phi ^{{\dagger} }$ is interpreted as a purely baroclinic response, that is, $\varPhi ^r=O(h^2/L^2)$ may be neglected (Appendix A) (Garrett & Gerkema 2007). The CMS (3.5) changes accordingly. Under this assumption and using (2.16), (2.24) reduces to the energy equation derived by 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 $\varPhi ^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$, (2.12a–c) and (3.5) 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 by 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 (3.5), we deduce that $\phi _m^{(0)}=0$ and that $\phi _m^{(1)}$ solves

(3.9)\begin{equation} \phi^{(1)}_{m,xx}+\ell_m^2\phi^{(1)}_m = 2 g_m \frac{r_{xx}}{h_0}\quad\text{for}\ m\geq 1, \end{equation}

with $\ell _m=m{\rm \pi} /(h_0\mu )$ and the radiation conditions in (3.7) with $k_m^+ = k_m^- =\ell _m$. Substituting the solution (C3) in (2.26), we obtain the conversion rate for weak topography,

(3.10)\begin{equation} C^{WTA} = \frac{F_0}{h_0^{2}}\sum_{n=1}^{\infty} |\hat{r}(\ell_n)|^2 \frac{n{\rm \pi}^2}{(\mu h_0)^2}\quad \text{with}\ F_0 = \frac{\left[(N^2-\omega^2)(\omega^2-f^2)\right]^{1/2}}{2{\rm \pi}\omega}U_0^2h_0^2, \end{equation}

where $\hat {r}(\xi ) = \int _{-\infty }^{+\infty }\exp (-\text {i} x\xi )r(s)\,{\rm d} s$ is the Fourier transform of $r$ (Appendix C). This is the formula given by Laurent et al. (2003) (up to multiplication with $\rho _0$), which was first derived by 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 (2.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, (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$th mode of the fundamental frequency, respectively, may be replaced by its asymptotic value $U_0k_{j1}/(2\omega _0)$ leading to our (3.10) with $\omega =\omega _0$.

3.4. Numerical solution for arbitrary topography

We truncate the infinite CMS (3.5)–(3.7) 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 fourth-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,\ldots,N_X$, $n=1,\ldots,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 (3.1). We then compute the baroclinic fields and conversion rates using (2.16) and (2.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_{G}$ and the ‘bump’ $h=h_0-h_{B}$, where

(4.1a,b)\begin{equation} h_{G}=\varLambda\exp\left(-\frac{x^2}{2L^2}\right)\quad\text{and}\quad h_{\text{B}}=\varLambda\exp\left(1-\frac{1}{1-x^2/L^2}\right)\mathbb{1}_{({-}L,L)}, \end{equation}

with $\mathbb {1}_{(-L,L)}=1\ \text {in}\ (-L,L)$ and $\mathbb {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 that in Griffiths & Grimshaw (2007, § 5), namely, $h=h_{S}$ with

(4.2)\begin{equation} h_{S} = \left\{\begin{array}{@{}l} h_-\quad\text{for}\quad x\leq 0,\\ h_-+ (h_+-h_-) \sin^2({\rm \pi} x/2L) \quad\text{for}\ 0\leq x\leq L,\quad \text{and}\\ h_+\quad\text{for}\ x\geq L. \end{array}\right. \end{equation}

We recall from § 2.1 that two important parameters will be considered, namely, the relative topography height $\delta$ and the criticality $\varepsilon$ 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}\,\text {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_+ = \varLambda$).

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 as $\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 § 3.4) to ensure that the numerical solution does not depend on them. We show the results on $\|\phi _n\|_{\infty }$ in figure 3(a). We see that the decay rate drops from $n^{-3}$, the theoretically expected rate in the case of a smooth solution (§ 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).

Figure 3. (a) Decay of $\|\phi _n\|_{\infty }$ in log-scale for the case of a bump topography with $\delta = 0.5$ and $\varepsilon = 0.7$, $0.85$, $0.90$, $0.95$ and $1.00$. (b) Normalised absolute error of the energy balance equation, $E$ = $|C_+-C_--C_{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 =0.7$ (blue surface) and $\varepsilon =1.0$ (orange surface).

Next, we consider the normalised error of the energy equation (2.24), $E = \lvert C_+- C_- - C_{int}\rvert /F_0$, where $C_{\pm }$ are calculated via (2.26) and $C_{int}$ via (2.27). Here, $C_{\pm }$ depends only on boundary values of $\phi _n$, whereas $C_{int}$ depends on their values in the entire computational domain. For an exact solution, $E$ would be zero. Therefore, monitoring $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$th internal wave mode over the ridge and $s= 4,6,8,10$ and $12$. We show results for $\varepsilon =0.7$ and $\varepsilon =1.0$ in figure 3(b). We obtain the expected $s^{-4}$ decay in both cases, verifying the fourth-order accuracy of the spatial discretisation scheme. For $\varepsilon =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 =1.0$, the error decays as $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$ (respectively, $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$ (respectively, $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. However, 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,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.

5. Visualisation of flow fields

5.1. Gaussian ridge

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

Figure 4. Gaussian ridges with $\delta =0.5$; (a–c) $\varepsilon =0.8$ and (d–f) $\varepsilon =1.2$. From top to bottom, non-dimensional baroclinic stream function $\psi ^{\#}/Q$, energy density $\mathcal {E}^{\#}/U_0^2$ and body-forcing term $\varepsilon ^2b^2\varPhi _{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$. Magenta 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 colour maps used in this and the subsequent figures are from Thyng et al. (2016).

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 nonlinearity 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 (2.13a), which couples the baroclinic motion with the non-uniform barotropic current and scales with $(1/h)_{xx}$ (figure 4c,f). This effect becomes more apparent if we write the PDE (2.13a) as a first-order system,

(5.1)\begin{equation} \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 ={-}\varPhi^{(0)}_{xx}. \end{equation}

Note that the above equation also holds with $\phi ^{{\dagger} }$ replaced by $\phi ^{\#}$ and $\varPhi ^{(0)}$ by $\varPhi$. Equation (5.1) shows that the solution is not constant along a given characteristic unless $\varPhi ^{(0)}_{xx}=0$. This is captured by our solution as an adjustment of the iso-energy curves near the ridge crest (figure 4b,e). In the insets of figure 4(b,e), we also plot the field $\phi ^{\#}_x+ \phi ^{\#}_z/\mu$ to highlight the variation of the beams. 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.

5.2. A non-radiating ridge

For further validation, we reproduce a case investigated by Maas (2011, § 4.2, ‘Ridge’). Maas considered the dimensionless version of (2.13a–c), obtained by the scaling $(x,z,\phi ^{{\dagger} },h) = (L\tilde {x},L/\mu \tilde {z},Q\widetilde {\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 Maas (2011, (4.15)), which is very similar to a Gaussian ridge with $\varLambda =0.19$ and $L=e/\sqrt {2}$ in our (4.1a,b), 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 of the order of $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 $\varPhi ^{(0)}_{xx}$. All variables are dimensionless.

6. Energy conversion rates

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

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 (3.10), using a log-log plot 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, (2.12)), using a linear plot to focus on the largest $\varepsilon$ values. For $\delta =0.1$, $C$ increases rapidly up to approximately $\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 nonlinear calculations. For $\delta =0.5$, $C$ monotonically increases with $\varepsilon$, approaching $C_{KE}$. For $\delta =0.9$, a local minimum is attained at approximately $\varepsilon = 1.4$, after which $C$ increases, also approaching $C_{KE}$ (inset of figure 6a). 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 that 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 precise. The horizontal domain considered describes between $99\,\%$ and $99.5\,\%$ of the variation of height of the topography. We refer to § 7 for a further discussion.

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^{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^{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 panels (c,d). The hatched region in panel (b) corresponds to topographies for which conversion is negligible for both the full solution and the WTA.

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

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^{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$ of 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 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 (4.2) with $h_-=2000$ m and $h_+=1000$ m ($\delta = 0.5$) for which calculations based on another modal decomposition method are reported by 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 figure 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 (Laurent et al. 2003). We obtained this result from Laurent et al. (2003, (32)–(36)) for 2000 modes. Note however that for $\delta =0.5$, the matrix in (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 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 the prediction of Laurent et al. (2003) 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 Green's function of Robinson (1969) 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 ^{G}(x,z) =\varPhi _0(z) + \vartheta (x,z)$, where $\varPhi _0=-Qz/h_0$ represents the barotropic flow corresponding to the uniform strip $[-\infty,+\infty ]\times [-h_0,0]$. From (2.8a–c), the field $\vartheta$ solves

(7.1a–c)\begin{equation} \mathcal{L}_{\mu}\vartheta = 0,\quad \vartheta(x,0)=0\quad\text{and}\quad \vartheta(x,-h)= Q(1-h/h_0), \end{equation}

together with conditions at $x\rightarrow \pm \infty$. Next, to solve (7.1a–c), $\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)$,

(7.2)\begin{equation} \vartheta(x,z)=\int_{{-}a}^{a}\gamma(x')G(x,x',z,-h(x'))\,\,\mathrm{d}\kern0.7pt x', \end{equation}

where it is additionally assumed that $-h_0+h$ is supported on $[-a,a]$. Note that (7.2) 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 (7.2) satisfies the first two equations in (7.1a–c) as well as the radiation conditions as $a\rightarrow +\infty$. The third equation in (7.1a–c) yields

(7.3)\begin{equation} \int_a^{{-}a}\gamma(x')G(x,x',-h(x),-h(x'))\,\,\mathrm{d}\kern0.7pt x' = Q(1-h/h_0), \end{equation}

which must be solved for $\gamma (x)$. Equation (7.3) is a Fredholm integral equation of the first kind, which is known to be ill-posed (see e.g. Groetsch (1984, § 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 (7.3) 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 $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 § 5, the baroclinic wave field varies significantly in regions where the body-forcing term is important (recall (5.1)). However, the ansatz (7.2) 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 (7.1a–c). We can reconcile the two approaches by verifying that the purely baroclinic field, satisfying (2.12a–c), is recovered by $\phi ^{\#} = \varPhi _0+\vartheta -\varPhi$, with $\varPhi = \varPhi ^{(0)}+\varPhi ^r$.

The second issue is associated with the the so-called ‘knife-edge limit’ established by Pétrélis et al. (2006), see also Echeverri & Peacock (2010). The conversion rate obtained via the ansatz (7.2) tends to this limit because (7.2) 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 ($\ldots$) 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 (2.13a–c) (or, equivalently, in (2.12a–c)) is not well defined because (2.9a–c) 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_0h_0$ everywhere except at the position of the knife-edge where it equals $U_0(h_0-\varLambda )$, with $\varLambda$ 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 2-D using the body-forcing formulation of Garrett & Gerkema (2007) 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 in (2.13a), 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 (respectively bump ridge) is typical for infinite-support (respectively 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 solution fields of the two formulations differ over sloping regions, and the solution obtained 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 have provided new numerical evidence that the conversion rate approaches the one obtained in the case of a step (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.