2.1. Governing equations

We start from the incompressible Euler equations under the Boussinesq approximation. The governing equations are

(2.1a) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} & = - \frac {\boldsymbol{\nabla }p + g \rho \boldsymbol{e}_z}{\rho ^\star } ,\end{align}

(2.1b) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho & = 0, \end{align}

(2.1c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} & = 0 ,\\[10pt]\nonumber \end{align}

where $\boldsymbol{u}(x,z,t) = (u, w)$ is the velocity vector, $p(x,z,t)$ is the pressure, $\rho (x,z,t)$ is the density, $\rho ^\star$ is the reference density, $g$ is the gravitational acceleration and $\boldsymbol{e}_z$ is the upward unit vector. For the sake of conciseness, differentiation is also denoted as $\partial u / \partial t = \partial _t u = u_t$ in the following. Fluid is bounded in the vertical direction by a flat top ( $z = H$ ) and a corrugated bottom ( $z = h(x)$ ) boundary, such that

(2.2) \begin{align} w (x, H, t) = 0, \quad w (x, h, t) = h_x u(x, h, t) \end{align}

are understood. Although the equation for the density (2.1b ) does not require boundary conditions in general, to simplify the problem, we shall assume that density is homogeneous along each boundary such that $\rho = \rho _{min}$ at $z = H$ and $\rho = \rho _{max}$ at $z = h(x)$ hold, where $\rho _{min}$ and $\rho _{max}$ are the minimum and maximum values of the density in the system, respectively. If this condition is initially fulfilled, it remains valid all the time.

We define the reference state as a combination of density stratification, vertically sheared horizontal flow and homogeneous oscillating tidal flow; i.e. $\rho = \rho _0 (z)$ , $u = u_0(z, t) \equiv U(z) - U_T \cos \omega _T t$ and $w = 0$ . The stratification is statically stable and the buoyancy frequency is defined as $N = \sqrt {(- g / \rho ^\star )\, {\rm d} \rho _0 / {\rm d} z}$ . Figure 1 illustrates the situation under consideration.

Figure 1.

Situation considered in this study. Internal gravity waves are generated over a spatially localised small bottom topography forced by an oscillatory barotropic flow in the presence of a sheared background current.

For a flat bottom boundary case, $h = 0$ , the reference state $(\boldsymbol{u}, \rho ) = (\boldsymbol{u}_0, \rho _0)$ is an exact solution of (2.1), but generally, the bottom topography acts as an obstacle for the current and a flow rising up the sloping boundary generates internal gravity waves. Therefore, we may write the solution as a superposition of the reference components and wave disturbances,

(2.3) \begin{align} u = u_0 + u', \quad w = w', \quad \rho = \rho _0 + \rho ', \quad p = p_0 + p' . \end{align}

Inserting these expressions into (2.1) and (2.2), and regarding the reference pressure as

(2.4) \begin{align} p_0 = g \int ^H_z \rho _0 (\zeta ) \,\text{d}\zeta - \rho ^\star \omega _T U_T x \sin \omega _T t , \end{align}

we derive a set of governing equations for the wave motion,

(2.5a) \begin{align} \frac {\partial \boldsymbol{u}'}{\partial t} + u_0 \frac {\partial \boldsymbol{u}'}{\partial x} + u_{0z} w' \boldsymbol{e}_x + \boldsymbol{u}' \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}' & = - \frac {\boldsymbol{\nabla }\!p' + g \rho ' \boldsymbol{e}_z}{\rho ^\star }, \end{align}

(2.5b) \begin{align} \frac {\partial \rho '}{\partial t} + u_0 \frac {\partial \rho '}{\partial x} - \frac {\rho ^\star N^2}{g} w' + \boldsymbol{u}' \boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho' & = 0, \end{align}

(2.5c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}' & = 0 ,\\[8pt]\nonumber \end{align}

as well as their boundary conditions,

(2.6) \begin{align} w' (x, H, t) = 0, \quad w' (x, h, t) = h_x \left [ u_0(h, t) + u'(x, h, t) \right ] . \end{align}

In the following, we shall assume that the bottom topography and the disturbance at the initial time are horizontally localised such that $\lim _{x \to \pm \infty } \{ h, u', w' \} = 0$ at $t = 0$ is understood. This condition, as well as the incompressible constraint (2.1c ), makes the net horizontal volume transport governed by the reference state and independent of the bottom-generated disturbance. This point is made explicit by a formula, $\int _{h(x)}^H u(x, z, t)\,\text{d}z = \int _0^H u_0(z, t) \,\text{d}z$ , or equivalently,

(2.7) \begin{align} \int ^H_h u' \,\text{d}z = \int _0^h u_0 \,\text{d}z , \end{align}

whose horizontal derivative conforms to the bottom boundary condition in (2.6).

2.2. Scaling and dimensionless constants

Up to now, we have not made any approximations; the set of equations governing the disturbance motions, (2.5) and (2.6), is equivalent to the original ones, (2.1) and (2.2). A key step to making the problem more tractable is linearising the equations based on proper scaling and asymptotic expansion. We write the typical velocity as $U^\star$ , the typical buoyancy frequency as $N^\star$ and the typical topographic height as $h^\star$ . We then scale and redefine the variables as

(2.8) \begin{align} x + \int ^t_0 U_T \cos \omega _T t' \text{d}t' = H \tilde {x}, \quad z = H \tilde {z}, \quad t = \frac {\tilde {t}}{N^\star } , \nonumber \\ (u', w') = U^\star \left (- \frac {\partial \tilde {\psi }}{\partial \tilde {z}}, \frac {\partial \tilde {\psi }}{\partial \tilde {x}} \right ), \quad \rho ' = \frac {\rho ^\star U^\star N^\star \tilde {\rho }}{g}, \quad p' = \rho ^\star U^\star N^\star H \tilde {p} , \nonumber \\ \left ( U, U_T \right ) = U^\star \big( \tilde {U}, \tilde {U}_T \big), \quad \omega _T = N^\star \tilde {\omega }_T, \quad N = N^\star \tilde {N}, \quad h = h^\star \tilde {h} , \end{align}

where the horizontal coordinate is changed to follow the background oscillating current, and $\tilde {\psi }$ and $\tilde {\rho }$ represent the streamfunction and the density perturbation for the wave component, respectively.

In the present setting, two dimensionless constants naturally arise as

(2.9) \begin{align} \epsilon \equiv \frac {N^\star h^\star }{U^\star }, \quad \mu \equiv \frac {U^\star }{N^\star H} , \end{align}

in which $\epsilon$ measures the nonlinearity of the topographically generated internal waves and $\mu$ represents the Froude number of the reference flow.

2.3. Asymptotic regime and expansion

To apply the asymptotic expansion, we assume a small topography regime $\epsilon \ll 1$ so that the nonlinear effects are negligible in the leading order. If we regard the full depth $H$ as the typical length scale of the reference velocity, $\mu ^{-2}$ roughly represents the gradient Richardson number. Accordingly, the assumption that $\mu$ is small would allow the use of the WKB approximation in the vertical direction. However, we do not rely on such an approximation here; from now on, we set $\mu =1$ , which means that the velocity variables are scaled by $N^\star H$ . For simplicity, we will omit the tilde on all symbols in the following.

We expand the unknown variables in terms of a small parameter $\epsilon$ as

(2.10) \begin{align} (\psi , \rho , p) = \epsilon (\psi _1, \rho _1, p_1) + \epsilon ^2 (\psi _2, \rho _2, p_2) + \cdots . \end{align}

Accordingly, from (2.5) to (2.7), the equations governing the first-order terms are derived as

(2.11a) \begin{align} \frac {\partial {\nabla} ^2 \psi _1}{\partial t} & = - \frac {\partial \rho _1}{\partial x} - U \frac {\partial {\nabla} ^2 \psi _1}{\partial x} + U_{zz} \frac {\partial \psi _1}{\partial x}, \end{align}

(2.11b) \begin{align} \frac {\partial \rho _1}{\partial t} & = N^2 \frac {\partial \psi _1}{\partial x} - U \frac {\partial \rho _1}{\partial x}, \\[10pt]\nonumber\end{align}

with the boundary conditions,

(2.12) \begin{align} \psi _1 (x, 1, t) = 0, \quad \psi _1 (x, 0, t) = \left [ U(0) - U_T \cos \omega _T t \right ] h \left ( x - \int ^t_0 U_T \cos \omega _T t' \,\text{d}t' \right ) . \end{align}

Note that the bottom boundary condition is now set at $z = 0$ instead of $z = \epsilon h$ . This treatment is common for topographic wave generation theory and valid as far as we consider linear processes. However, it requires some modifications in the energy budget equation. We will discuss this point in § 5.

Compared with the original nonlinear equations, the reduced linear equations, (2.11) and (2.12), are simple enough to be solved analytically. In the following, we assume $N(z)$ and $U(z)$ are sufficiently smooth functions and satisfy the condition for the Richardson number $N^2 / U_z^2 \gt 1/4$ everywhere in $0 \leq z \leq 1$ , which ensures the flow stability (Howard Reference Howard1961; Miles Reference Miles1961). In addition, we also assume $U_z\gt 0$ to let $N^2 / U_z^2$ be finite.

3.1. Horizontal Fourier transform

Since the coefficients in (2.11) do not depend on $x$ , it is useful to take the Fourier transform in the horizontal direction,

(3.1) \begin{align} \left [ \psi _1(x,z,t), \rho _1(x,z,t), h(x) \right ] = \frac {1}{2\pi } \int ^\infty _{-\infty } \big[ \hat {\psi }_1(k,z,t), \hat {\rho }_1(k,z,t), \hat {h}(k) \big ] e^{\textrm {i} kx} \,\text{d}k . \end{align}

Following Bell (Reference Bell1975), we apply the Jacobi–Anger expansion to the Fourier transform of the bottom boundary condition in (2.12) in terms of the harmonics of the tidal frequency to obtain

(3.2) \begin{align} \hat {\psi }_1 (k, 0, t) = \sum _{n = -\infty }^\infty \frac {(k U_0 - \omega _n) \hat {h}}{k} J_n \left ( \frac {k U_T}{\omega _T} \right ) e^{- \textrm {i} \omega _n t} , \end{align}

where $U_0 \equiv U(0)$ , $\omega _n \equiv n \omega _T$ and $J_n$ is the $n$ th-order Bessel function of the first kind (Equation (3.2) involves an infinite number of terms corresponding to the harmonics of the tidal frequency. Their relative contributions depend on the parameters $k U_T / \omega_T$ and $k U_0 / \omega_T$ . If the steady current is absent at the bottom ( $U_0 = 0$ ), it is classically known that the tidal frequency $\vert n \vert = 1$ is dominant when $k U_T / \omega_T \ll 1$ , while high harmonics, $\vert n \vert \gt 1$ , are important when $k U_T / \omega_T \gtrsim 1$ . On the other hand, if we take into account the steady current, $U_0 \neq 0$ , a stationary component, $n = 0$ , also appears (Shakespeare Reference Shakespeare2020) and is dominant when $k U_T / \omega_T \ll 1$ .). The governing equations (2.11) are transformed into

(3.3a) \begin{align} \frac {\partial {\nabla} _k^2 \hat {\psi }_1}{\partial t} & = - \textrm {i} k \hat {\rho }_1 - \textrm {i} k U {\nabla} _k^2 \hat {\psi }_1 + \textrm {i} k U_{zz} \hat {\psi }_1 ,\end{align}

(3.3b) \begin{align} \frac {\partial \hat {\rho }_1}{\partial t} & = \textrm {i} k N^2 \hat {\psi }_1 - \textrm {i} k U \hat {\rho }_1,\\[7pt]\nonumber \end{align}

with ${\nabla} ^2_k \equiv -k^2 + \partial _z^2$ . Equations (3.2) and (3.3), as well as the upper boundary condition, $\hat {\psi }_1(k, 1, t) = 0$ , redefine the problem to be solved. Classically, there are two strategies to address this type of problem. The first is to directly find the solution that satisfies the inhomogeneous boundary condition. The second is to decompose the solution into two parts: one that satisfies the inhomogeneous boundary condition and is easily detectable, and the other that satisfies homogeneous boundary conditions and can be solved separately. This section employs the second approach because it provides a better insight into the basic spectral properties of waves inherent in shear flow problems. The first strategy will be applied in § 4, where we discuss the stationary response solution.

3.2. Potential and vortical flow decomposition

To handle the inhomogeneous boundary condition at the bottom, we employ the Helmholtz decomposition for the velocity field to write the streamfunction as $\hat {\psi }_1 = \hat {\psi }_1^p + \hat {\psi }_1^v$ . Here, $\hat {\psi }_1^p$ and $\hat {\psi }_1^v$ represent the potential and vortical flow parts, respectively. The potential flow solves ${\nabla} ^2_k \hat {\psi }_1^p = 0$ and satisfies the inhomogeneous bottom boundary condition (3.2) and the homogeneous upper boundary condition. An analytical solution is immediately derived as

(3.4) \begin{align} \hat {\psi }^p_1 = \frac {\sinh \left [ k(1 - z) \right ]}{\sinh k} \sum _{n = -\infty }^\infty \frac {(k U_0 - \omega _n) \hat {h}}{k} J_n \left ( \frac {k U_T}{\omega _T} \right ) e^{- \textrm {i} \omega _n t} . \end{align}

The vertical structure of the potential flow depends on the magnitude of the horizontal wavenumber. When the wavenumber is small, $\vert k \vert \ll 1$ , the horizontal component of the potential flow, $-\partial _z \hat {\psi }^p_1$ , is homogeneous in the vertical direction, representing the barotropic response in the hydrostatic regime. If the wavenumber is large, $\vert k \vert \gg 1$ , the potential flow is localised close to the bottom and exponentially weak in the upper part of the fluid.

Now that a solution of the potential flow part is obtained, the remaining vortical flow part $\hat {\psi }_1^v$ as well as the density perturbation $\hat {\rho }_1$ are determined from (3.3) as the solutions of a set of inhomogeneous equations,

(3.5a) \begin{align} \frac {\partial {\nabla} _k^2 \hat {\psi }^v_1}{\partial t} & = - \textrm {i} k \hat {\rho }_1 - \textrm {i} k U {\nabla} _k^2 \hat {\psi }^v_1 + \textrm {i} k U_{zz} \hat {\psi }^v_1 + \textrm {i} k U_{zz} \hat {\psi }^p_1, \end{align}

(3.5b) \begin{align} \frac {\partial \hat {\rho }_1}{\partial t} & = \textrm {i} k N^2 \hat {\psi }^v_1 - \textrm {i} k U \hat {\rho }_1 + \textrm {i} k N^2 \hat {\psi }^p_1 ,\\[10pt]\nonumber \end{align}

with the homogeneous boundary conditions, $\hat {\psi }^v_1 (k, 1, t) = \hat {\psi }^v_1 (k, 0, t) = 0$ . Our present problem is thus to investigate the linear response of prognostic variables $(\hat {\psi }^v_1, \hat {\rho }_1)$ to time-dependent external forcing originating from $\hat {\psi }_1^p$ as governed by (3.5). Equations (3.5) are valid for each $k$ independently. Since all the terms on the right-hand sides vanish when $k=0$ , we may set $\hat {\psi }_1^v(0, z, t) = \hat {\rho }_1(0, z, t) = 0$ .

3.3. Formal solutions derived from Laplace transform

To derive the general solution of (3.5), it is convenient to introduce a linear operator $\unicode{x1D648}_k$ and rewrite the system as a forced Schrödinger-type evolution equation

(3.6) \begin{align} \textrm {i} \frac {\partial \hat {\boldsymbol{v}}}{\partial t} = \unicode{x1D648}_k \hat {\boldsymbol{v}} + \sum _{n = -\infty }^\infty \hat {\boldsymbol{f}}_{\!n} e^{- \textrm {i} \omega _n t} , \end{align}

where

(3.7) \begin{align} \hat {\boldsymbol{v}}(k, z, t) = \left (\! \begin{array}{c} {\nabla} _k^2 \hat {\psi }^v_1 \\ \hat {\rho }_1 \end{array} \!\right ), \quad \hat {\boldsymbol{f}}_{\!n} (k, z) = - \frac {(k U_0 - \omega _n) \hat {h} \sinh \left [ k(1 - z)\right ]}{\sinh k} J_n \left ( \frac {k U_T}{\omega _T} \right ) \left (\!\begin{array}{c} U_{zz} \\ N^2 \end{array}\!\right ) \end{align}

and

(3.8) \begin{align} \unicode{x1D648}_k = \left (\!\begin{array}{cc} k U - k U_{zz} {\nabla} _k^{-2} & k \\[3pt] - k N^2 {\nabla} ^{-2}_k & k U \end{array} \!\right )\!. \end{align}

Here, the inverse of the Laplacian is defined by ${\nabla} _k^{-2} A = B \Leftrightarrow A = {\nabla} _k^2 B$ with $\left . B \right |_{z = 0,1} = 0$ for a prescribed $A$ , and its explicit representation is given in (A2). For an arbitrary initial condition $\hat {\boldsymbol{v}}(k,z,0) = \hat {\boldsymbol{v}}_0(k, z)$ , the solution of (3.6) can be formally obtained using a temporal Laplace transform. We define

(3.9) \begin{align} \tilde {\boldsymbol{v}} (k, z; \sigma ) = \int _0^{\infty } \hat {\boldsymbol{v}} (k, z, t) e^{\textrm {i} \sigma t}\, \text{d} t , \end{align}

where $\sigma$ is a complex number with $\Im \sigma$ sufficiently large for the integral to converge. Applying this transform to (3.6) yields

(3.10) \begin{align} \sigma \tilde {\boldsymbol{v}} = \unicode{x1D648}_k \tilde {\boldsymbol{v}} + \textrm {i} \hat {\boldsymbol{v}}_0 + \sum _{n = -\infty }^\infty \frac {\textrm {i} \hat {\boldsymbol{f}}_{\!n}}{(\sigma - \omega _n)} . \end{align}

The inverse transform then gives the general solution

(3.11) \begin{align} \hat {\boldsymbol{v}} = \frac {\textrm {i}}{2 \pi } \int _{-\infty + \textrm {i} \gamma }^{\infty + \textrm {i} \gamma } \frac {e^{- \textrm {i} \sigma t} \hat {\boldsymbol{v}}_0}{\sigma \unicode{x1D644} - \unicode{x1D648}_k} \,\text{d}\sigma + \frac {\textrm {i}}{2 \pi } \sum _{n = - \infty }^\infty \int _{-\infty + \textrm {i} \gamma }^{\infty + \textrm {i} \gamma } \frac {e^{- \textrm {i} \sigma t} \hat {\boldsymbol{f}}_{\!n}}{(\sigma - \omega _n) (\sigma \unicode{x1D644} - \unicode{x1D648}_k)} \,\text{d}\sigma , \end{align}

where $\unicode{x1D644}$ is the identity matrix and $\gamma \gt 0$ is chosen so that (3.9) converges for all $\sigma$ with $\Im \,\sigma \gt \gamma$ .

Equation (3.11) expresses the solution as an inverse Laplace transform whose integrand contains $(\sigma \unicode{x1D644}-\unicode{x1D648}_k)^{-1}$ , the resolvent of the operator $\unicode{x1D648}_k$ . In general, such inverse transforms are evaluated by deforming the contour and accounting for the singular behaviour of the integrand in the complex $\sigma$ -plane. It is therefore essential to identify the set of values of $\sigma$ for which the resolvent fails to be well defined. For systems with finitely many degrees of freedom, the relevant singularities typically reduce to isolated poles, so the residue theorem yields the familiar decomposition into free modal oscillations and a particular forced response. In the present shear-flow problem, however, the unknowns are functions of the continuous vertical coordinate $z$ and the resolvent can exhibit non-isolated singular behaviour that cannot be captured by residues alone. We refer to this set as the spectrum of $\unicode{x1D648}_k$ , which we determine in the next subsection.

3.4. Spectrum of the operator $\unicode{x1D648}_k$

To make the previous statement concrete, we regard the resolvent $(\sigma \unicode{x1D644}-\unicode{x1D648}_k)^{-1}$ as a function of the complex parameter $\sigma$ . For each fixed $(k,\sigma )$ , (3.10) is a boundary-value problem in $z$ . When this problem admits a regular solution in the entire domain $0 \leq z \leq 1$ , the resolvent is well defined. Conversely, values of $\sigma$ for which the boundary-value problem fails to have a regular solution correspond to singularities of the resolvent. We denote the set of such values of $\sigma$ for a given $k$ by $\varSigma _k$ and call it the spectrum of $\unicode{x1D648}_k$ .

The spectrum $\varSigma _k$ is directly assessed by making explicit the resolvent of $\unicode{x1D648}_k$ , or equivalently, solving the boundary-value problem (3.10) for each $(k,\sigma )$ . In practice, this programme reduces to solving a second-order differential equation. This point is made clear by reorganising (3.10) in terms of the streamfunction as

(3.12) \begin{align} \mathcal{T}_{k,\sigma } \tilde {\psi } = F \quad \text{with} \quad \mathcal{T}_{k,\sigma } \equiv \frac{\text{d}^2}{\text{d} z^2} + \frac {k^2 N^2}{(k U - \sigma )^2} - \frac {k U_{zz}}{k U - \sigma } - k^2 , \end{align}

where $F(z)$ collects the inhomogeneous terms involving $\hat {\boldsymbol{v}}_0$ and $\hat {\boldsymbol{f}}_{\!n}$ . The homogeneous counterpart of this equation, $\mathcal{T}_{k,\sigma }\tilde {\psi }=0$ , is the Taylor–Goldstein equation governing unforced internal waves in shear flows. For the inhomogeneous problem (3.12), the solution can be written as an integral representation in terms of the corresponding Green function (e.g. Bender & Orszag Reference Bender and Orszag2013). We summarise this Green-function solution in Appendix A; for our purposes, the key point is that it becomes singular for particular values of $\sigma$ and these values constitute $\varSigma _k$ .

As elucidated in Appendix A, the spectrum $\varSigma _k$ consists of two parts with distinct physical origins. First, isolated values of $\sigma \in \varSigma _k$ correspond to regular eigenvalue solutions of the Taylor–Goldstein equation. We collect such values of $\sigma$ to define a subset of the spectrum $\varSigma ^d_k = \{\sigma \,\vert \, D(k,\sigma )=0\}$ , where $D(k,\sigma )$ specifies the dispersion relation of vertically standing internal waves and is obtained by enforcing the boundary conditions on the Taylor–Goldstein equation (see Appendix A.2.1). Since these solutions are regular throughout the domain, the Miles–Howard stability theorem applies and ensures that $\varSigma ^d_k \subset \mathbb{R}$ under the assumption $N^2/U_z^2 \gt 1/4$ for all $0 \leq z \leq 1$ .

Second, when $\sigma /k$ lies within the range of the background flow $U(z)$ , the Taylor–Goldstein equation admits singular solutions associated with a critical level at which $kU=\sigma$ (Appendix A.2.5). The corresponding values of $\sigma$ then form a continuous interval, $\varSigma ^c_k = \left \{\sigma \,\vert \, U(0) \leq \sigma /k \leq U(1)\right \}$ . In the following, we refer to $\varSigma ^d_k$ and $\varSigma ^c_k$ as the discrete and continuous spectra, respectively.

Figure 2.

Integration contour of the inverse Laplace transform for a positive $k$ . We shift the contour originally defined in the upper half of the complex plane in the negative direction along the imaginary axis. There exist an infinite number of poles, which involve the eigenvalues of $\unicode{x1D648}_k$ , i.e. the discrete spectrum $\varSigma ^d_k$ , and the forcing frequency $\omega _n$ . An element of the discrete spectrum, $\sigma \in \varSigma ^d_k$ , is either in a range $\sigma \lt k U(0)$ or $\sigma \gt k U(1)$ , while the forcing frequency $\omega _n$ may be located anywhere on the real axis. There also exist three branch points on the real axis. Two branch points $kU(0)$ and $kU(1)$ originate from the factors $\gamma ^\pm (0)$ and $\gamma ^\pm (1)$ , and the intermediate point $kU(z)$ from $\gamma ^\pm (z)$ , which are involved in the solutions of the Taylor–Goldstein equation, (A4). Since the location of the intermediate point depends on $z$ and thus ranges from $k U(0)$ to $k U(1)$ , the continuous spectrum, $\varSigma ^c_k$ , is composed of $\sigma \in [ k U(0), kU(1)]$ .

To examine how each element in the spectrum $\varSigma _k$ affects the solution (3.11), we consider the integration contour illustrated in figure 2; we shift the contour in the negative direction of the imaginary axis, but circumventing the singular points (elements of $\varSigma _k$ ). The net contribution to the integration is constituted of the residues from poles and integration along branch cuts. In the first term in (3.11), which represents the free wave solution, the residues of $\sigma \in \varSigma ^d_k$ represent standing waves that are the discrete eigenmodes of the system, in agreement with the fact that $D(k, \omega ) = 0$ represents the dispersion relations of these modes. In contrast, integration along a branch cut that extends from a point in $\varSigma ^c_k$ originates from the steady shear flow. There exist three branch points as indicated in figure 2, among which the middle one dependent on $z$ reflects the critical level. This point appears in the integrand in the form

(3.13) \begin{align} \left [ \sigma - k U(z) \right ]^{m + 1/2 \pm \textrm {i} \nu } , \end{align}

where $m$ is an integer, and $\nu$ is a function of $\sigma$ and $k$ (see Appendix A). Because the exponent is not an integer, the residue theorem is not applicable there. An asymptotic evaluation shows that the contributions from the branch cuts behave algebraically in time and the wave energy will not grow from a smooth initial profile in the long-term limit (Brown & Stewartson Reference Brown and Stewartson1980; Jose et al. Reference Jose, Roy, Bale and Govindarajan2015). The unforced system is, therefore, stable.

The integrand of the second term in (3.11) possesses another singularity at the forcing frequency, $\sigma = \omega _n$ . If $\omega _n$ does not belong to the spectrum $\varSigma _k$ , the residue of this pole yields harmonic oscillation and, hence, the system is still stable. However, if the forcing frequency matches an element of the spectrum, integration around this point can exhibit algebraic growth in time. To focus on this resonant process, it is convenient to add a small relaxation term that mimics viscosity to the time derivative, setting $\partial _t \to \partial _t + r$ with $r \gt 0$ in (3.6). In this system, $\varSigma _k$ is slightly lowered along the imaginary axis so that the contributions to the integral in (3.11) from the spectrum are exponentially damped with time. The residue from $\sigma = \omega _n$ is solely responsible for the long-term behaviour. Accordingly, the stationary response solution is derived as

(3.14) \begin{align} \hat {\boldsymbol{v}} = \sum _{n = - \infty }^\infty \frac {e^{- \textrm {i} \omega _n t} \hat {\boldsymbol{f}}_{\!n}}{ (\omega _n + \textrm {i} r) \unicode{x1D644} - \unicode{x1D648}_k} . \end{align}

This solution is regular everywhere for a finite $r$ , but when taking the inviscid limit $r \to 0$ , singularity can arise. Specifically, if the forcing frequency $\omega _n$ belongs to the discrete spectrum $\varSigma ^d_k$ , the solution diverges almost everywhere in $z$ . On the other hand, if $\omega _n$ is an interior point of the continuous spectrum $\varSigma ^c_k$ , the solution is singular at a particular level (Appendix A.2.5). In either case, wave energy $E$ diverges – energy is unboundedly supplied from the forcing term. The spectrum inherent in the unforced system determines the condition of the permanent energy supply from the external tidal force.

Although we are mainly interested in the inviscid problem, keeping a small but finite $r \gt 0$ throughout the paper is useful in that it moves the location of a singular point slightly from a real axis in the $z$ - or $k$ -plane. This enables us to: (i) display the vertical structure of a solution that possesses a critical level (figure 3 d); (ii) determine in which half-plane a pole or branch point exists (figure 4); (iii) learn the behaviour of a solution in the genuine resonant situation (figure 5); and (iv) derive the explicit formula for the energy production rate (§ 5.3).

4.1. Single-wavenumber solution

Now, we inherit the idea of the previous section that the model involves a small relaxation term specified by $r \gt 0$ and the system reaches a stationary state. Consequently, the solution in the inviscid limit consists of the linear superposition of harmonic oscillations as

(4.1) \begin{align} \big [ \hat {\psi }_1(k, z, t), \hat {\rho }_1(k, z, t) \big ] = \lim _{r \to 0} \sum _{n = -\infty }^\infty \big [\hat {\psi }_{1n} (z; k, r), \hat {\rho }_{1n} (z; k, r) \big ] e^{- \textrm {i} \omega _n t}. \end{align}

Each element of the streamfunction solves the homogeneous Taylor–Goldstein equation, $\mathcal{T}_{k,\omega _n + \textrm {i} r} \hat {\psi }_{1n} = 0$ , with the inhomogeneous bottom boundary condition (3.2) and the homogeneous upper boundary condition. The density perturbation is then related via (3.3b ). We readily write down the solutions as

(4.2a) \begin{align} \hat {\psi }_{1n} (z; k, r) & = \frac {(k U_0 - \omega _n) \hat {h}(k)}{k} J_n \left ( \frac {k U_T}{\omega _T} \right ) \frac {A_1(z; k, \omega _n + \textrm {i} r)}{D(k, \omega _n + \textrm {i} r)}, \end{align}

(4.2b) \begin{align} \hat \rho _{1 n}(z; k, r) & =\frac {k N^2}{kU(z)-(\omega _n+ \textrm {i} r)} \hat \psi _{1 n} (z; k, r) ,\\[10pt]\nonumber \end{align}

where $A_1$ represents a vertical structure function that is a solution for $\mathcal{T}_{k,\omega _n + \textrm {i} r} A_1 = 0$ with a one-side boundary condition $\left . A_1 \right \vert _{z=1} = 0$ and $D$ plays the role of normalisation such that $\left . A_1 \right \vert _{z=0} = D$ is understood. Since the solutions in physical space, $\psi$ and $\rho$ , are real functions, their Fourier transform should satisfy the symmetric property as to changes in signs of wavenumber and frequency, i.e. $\hat {\psi }_{1n} (z; k, r) = \hat {\psi }^\dagger _{1-n} (z; - k, r)$ and $\hat {\rho }_{1n} (z; k, r) = \hat {\rho }^\dagger _{1-n} (z; - k, r)$ , where ${}^\dagger$ denotes the complex conjugate. Therefore, it is enough to consider the cases of $n \geq 0$ in the following.

Figure 4.

Integration contour of the Fourier integral (3.1) for positive $x$ . The contour originally defined on the real axis is shifted over a finite distance on the complex plane along the imaginary axis. For a wavenumber $k = k_r + \textrm {i} k_i$ , the exponential factor in the integrand becomes $e^{\textrm {i} k x} = e^{\textrm {i} k_r x - k_i x}$ . Therefore, for a large $x$ , contributions from the line segments in the upper half of the complex plane become exponentially small, and the integration is dominated by those around the poles and branch points, which approach the real axis in the inviscid limit $r \to 0$ . For negative $x$ , the integration contour is shifted in the opposite direction. Note that the branch point around which the solution (4.11) is evaluated is $k' = (\omega _n + \textrm {i} r) / U(z)$ . Currently, we are considering the case when $U\gt 0$ for any $z \in [0, 1]$ . For $U \lt 0$ , the corresponding branch point exists in the lower half-plane.

The equation $\mathcal{T}_{k,\sigma } A_1 = 0$ can be solved for a general set of $N(z)$ and $U(z)$ numerically, but for the simplest case of uniform shear and stratification, $N=1$ and $U = U_0 +{} (U_1 - U_0) z$ , an analytical expression is available (Appendix B). Figure 3 shows the spectrum of $\unicode{x1D648}_k$ as well as some examples of the vertical structure function $A_1 (z)$ where we set $U_0 = 0$ and $U_1 = 0.15$ . In these plots, we augment the structure function by multiplying $e^{\textrm {i} k x}$ to illustrate both the vertical and horizontal structures. Depending on the horizontal wavenumber and the forcing frequency, the wave structure exhibits distinct characteristics. For $(k, \sigma ) = (-5, 1.5)$ as in figure 3(b), because the generated wave’s frequency is greater than the background buoyancy frequency, the solution exhibits an evanescent structure localised close to the bottom. If we choose a wavenumber and frequency close to a dispersion curve of a discrete mode (i.e. the discrete spectrum $\varSigma ^d_k$ ), specifically $(k, \sigma ) = (2.5, 0.6)$ in figure 3(c), the solution becomes a vertically standing wave structure (we cannot define a stationary solution if the forcing frequency lies in the dispersion curves since a denominator $D$ in (4.2) vanishes in that limit). There exist an infinite number of such standing modes in the vicinity of the edges of the continuous spectrum (Appendix A.2.3).

By contrast, when the forcing frequency lies within the continuous spectrum, $\omega _n \in \varSigma ^c_k$ , the inviscid stationary solution exhibits a singularity at a critical level where $k U$ coincides with $\omega _n$ . More specifically, according to (A4) and (A8), the vertical structure function $A_1$ admits a local representation

(4.3) \begin{align} A_1 (z; k, \omega _n + \textrm {i} r) & = (z - z^c)^{1/2 - \textrm {i} \nu } \alpha _1(z) + (z - z^c)^{1/2 + \textrm {i} \nu } \alpha _2(z) , \end{align}

where $z^c$ and $\nu$ are functions of $k$ and $\sigma = \omega _n + \textrm {i} r$ defined through

(4.4) \begin{align} k U(z^c) = \sigma \quad \text{and} \quad \nu = \sqrt {\left . \frac {N^2}{U_z^2} \right \vert _{z = z^c} - \frac {1}{4}} , \end{align}

and $\alpha _1(z)$ and $\alpha _2(z)$ are smooth functions. In the inviscid limit $r \to 0$ , the branch point $z = z^c$ approaches the real axis and the two local solutions in (4.3) are connected across $z^c$ according to the rule (A10). The upper-boundary condition then imposes $\alpha _2 = \alpha _1^\dagger$ . Above the critical level, since the two terms in (4.3) have the same amplitude, $A_1$ exhibits a standing wave structure. Below the critical level, the amplitudes of the two terms differ by a factor $e^{2 \pi \nu }$ , so that the upward-propagating component dominates the downward-propagating component. The wave amplitude in the upper layer is smaller than that of the lower layer by a factor $e^{\pi \nu }$ . These features are evident in figure 3(d), $(k, \sigma ) = (5, 0.5)$ .

This amplitude jump across the critical level is a classical signature of wave focusing and absorption into the background flow at a critical layer, most transparently discussed in initial-value problems (Booker & Bretherton Reference Booker and Bretherton1967). The previous discussion, however, assumes a single horizontal wavenumber, which corresponds to a sinusoidal topography. For localised topography, the forcing excites a continuum of horizontal wavenumbers, and the manifestation of the critical-layer behaviour is qualitatively different, as we discuss in the following.

4.2. Wave generation from localised topography

We now turn to wave generation by bottom topography that is localised in physical space. In this case, the full solution is obtained by performing the inverse Fourier transform, i.e. inserting (4.2) and (4.1) into (3.1). This step requires some care compared with the inverse Laplace transform discussed earlier. Even when the integrand is analytic in $k$ along the real axis, shifting the contour arbitrarily far into the complex plane is generally inappropriate because the analytic structure of $\hat {h}$ must also be taken into account. For instance, for a Gaussian ridge $h(x) \propto e^{-x^2/2}$ , the Fourier transform behaves as $\hat {h}(k) \propto e^{-k^2/2}$ and diverges as $k \to \pm \textrm {i} \infty$ . Closed-form expressions are therefore available only for special choices of $U$ , $N$ and $h$ . Here, we instead adopt an asymptotic approach to describe the solution far from the localised topography. As illustrated in figure 4, the integration contour is shifted a finite distance into the complex plane while circumventing the relevant singularities. In the limit $|x| \to \infty$ , the dominant contributions arise from singularities located near the real axis, since the integrand decays exponentially away from it. This contour-deformation argument applies to a broad class of spatially localised topographies.

Figure 5.

Spectrum of an operator $\unicode{x1D648}_k$ for a constant shear and stratification case, $N = 1$ and $U = 0.05 +{} 0.15 z$ . On the dispersion curves of discrete modes, there exist points where group velocity $c^g_j$ vanishes.

In the present problem, the relevant singularities are poles and branch points, which originate from the discrete and continuous spectra of the system, respectively. We first investigate the contributions from the poles. For a prescribed real $k$ , the location of the zeros of $D(k, \sigma )$ is specified by the discrete spectrum as $\varSigma ^d_k = \{ \sigma ^\pm _m \}_{m = 1, 2, \ldots }$ , where $\sigma ^+$ and $\sigma ^-$ are distinguished in terms of the sign of $\sigma / k - U$ (Appendix A.2). Similarly, for a prescribed real $\sigma$ , we shall write the location of the zeros of $D(k, \sigma )$ on the real axis as $K^d_\sigma = \{ k_j \}_{j = 1, 2, \ldots }$ . When $r \gt 0$ is sufficiently small, in the vicinity of $k = k_j \in K^d_{\omega _n}$ , $D$ is locally expanded as

(4.5a) \begin{align} D(k, \omega _n + \textrm {i} r) & = - \left [ k - k_j - \frac {\textrm {i} r}{c^g_j} + \mathcal{R}(k, r) \right ] c^g_j D_\sigma\\[-10pt]\nonumber \end{align}

(4.5b) \begin{align} \mbox{with} \quad c^g_j & \equiv - \frac {D_k}{D_\sigma } = \left . \frac {\partial \sigma }{\partial k} \right \vert _{D=0} ,\\[10pt]\nonumber \end{align}

where $\mathcal{R}$ is a residual factor that vanishes in the limits of $r \to 0$ and $k \to k_j$ faster than the other terms, i.e. $\mathcal{R} = o(k - k_j, r)$ in Landau notation. Expression (4.5b ) represents the slope of the dispersion curve, i.e. the horizontal group velocity of an eigenmode. According to (4.5a ), the integrand of the Fourier integral (3.1) involves a pole at $k = k_j + \textrm {i} r / c_j^g + o(r) \equiv k'_j$ . If $c_j^g \gt 0$ , this pole is located slightly above the real axis and the circular integration around it contributes to the solution for $x \to + \infty$ . If $c_j^g \lt 0$ , by contrast, the pole is located below the real axis and it contributes when $x \to - \infty$ . We shall write these factors involved in the Fourier integral in (3.1) as

(4.6) \begin{align} \lim _{r \to 0} \frac {1}{2 \pi } \oint _{k \sim k'_j} \hat {\psi }_{1n} e^{\textrm {i} k x}\,\text{d}k = \frac {(k_j U_0 - \omega _n) \hat {h}(k_j)}{\textrm {i} k_j \big\vert c^g_j \big\vert D_\sigma (k_j, \omega _n)} J_n \left ( \frac {k_j U_T}{\omega _T} \right ) A_1(z; k_j, \omega _n) e^{\textrm {i} k_j x} , \end{align}

where integration is performed in the counterclockwise (clockwise) direction for $c^g_j \gt 0$ ( $c^g_j \lt 0$ ). Equation (4.6) represents a vertically standing-mode wave train that travels far away from the topography while not changing its amplitude. This is the solution discussed by Lamb & Dunphy (Reference Lamb and Dunphy2018).

Equations (4.5a ) and (4.6) are invalid when the group velocity vanishes, $c^g_j = 0$ . An example of this situation is illustrated in figure 5. An alternative expansion of $D$ in this case is

(4.7) \begin{align} D(k, \omega _n + \textrm {i} r) = - \left [ \frac {(k - k_j)^2}{2} - \frac {\textrm {i} r}{c^g_{j k}} + \mathcal{R}(k, r) \right ] c^g_{j k} D_\sigma , \end{align}

where $\mathcal{R} = o((k - k_j)^2, r)$ and $c^g_{j k} = - D_{kk} / D_\sigma$ is understood. Supposing $c^g_{jk} \lt 0$ , the integrand in (4.7) has a pair of poles at

(4.8) \begin{align} k = k_j \pm (1 - \textrm {i}) \sqrt {\frac {r}{\big\vert c^g_{jk} \big\vert }} + o\big( r^{1/2} \big) \equiv k^\pm _j . \end{align}

Circular integration around these poles for small $r$ is evaluated as

(4.9) \begin{align} \frac {1}{2 \pi } \oint _{k \sim k^\pm _j} \hat {\psi }_{1n} e^{\textrm {i} k x} \,\text{d}k & = \frac {(k_j U_0 - \omega _n) \hat {h}(k_j)}{(1 + \textrm {i}) k_j \sqrt {r \big\vert c^g_{j k} \big\vert } D_\sigma (k_j, \omega _n)} J_n \left ( \frac {k_j U_T}{\omega _T} \right ) A_1(z; k_j, \omega _n) e^{\textrm {i} k^\pm _j x} \nonumber \\ \quad & + o \big(r^{-1/2} \big) . \end{align}

Accordingly, the stationary solution is divergent in the inviscid limit with the power law $\vert \psi \vert \propto r^{-1 /2}$ . Physically, if the group velocity is 0, wave energy accumulates locally above the topography. Interference between the accumulating wave and the topographic forcing results in resonant amplification (this situation does not correspond to a convective or absolute instability, as no exponentially growing modes exist when $Ri\gt 1/4$ (Miles–Howard theorem). Instead, the vanishing group velocity indicates a resonance condition, where the forced response is amplified due to the matching of the forcing frequency with the natural frequency of a stationary wave mode). Far from the topography, the wave amplitude is damped exponentially over the length scale $\sqrt {\vert c^g_{jk} \vert / r}$ according to (4.8). In this process, the group velocity dispersion represented by the curvature of the dispersion curve, $c^g_{jk}$ , plays distinctive roles. When $\vert c^g_{jk} \vert$ is small, a broad range of wavenumbers is involved in the resonance so that the wave is much more amplified. At the same time, the group velocity of these side-band components is still small, which results in shortening the decay length scale.

We next consider the contributions from the branch cut. For a prescribed $k$ and $\sigma$ , $A_1(z; k, \sigma )$ can have a branch point at $z = z'$ where $k U(z') = \sigma$ . In the same way, for a prescribed $\sigma$ and $z$ , a branch point $k = k'$ is defined where $k' U(z) = \sigma$ is satisfied. At a particular level where $U(z) = 0$ , however, this consideration does not make sense because $k U(z) = \sigma$ is identically satisfied for any $k$ if $\sigma = 0$ . This location is the critical level of the steady component, $n = 0$ , where an inviscid stationary solution is undefinable. In the following, we neglect this level and set $U(z) \neq 0$ and $n \gt 0$ .

Now, we write the branch point corresponding to the critical level as $k' (z, r) ={} (\omega _n + \textrm {i} r) / U(z)$ . Around this point, inserting $z = z^c (k', \omega _n + \textrm {i} r)$ and $z^c = z^c (k, \omega _n + \textrm {i} r)$ into (4.3), we expand $A_1$ as

(4.10) \begin{align} A_1(z; k, \omega _n + \textrm {i} r) & = \left ( \frac {U}{k' U_z} \right )^{1/2 - \textrm {i} \nu } \big(k - k' \big)^{1/2 - \textrm {i} \nu } \alpha _1(z) \nonumber \\ &\quad + \left ( \frac {U}{k' U_z} \right )^{1/2 + \textrm {i} \nu } \big(k - k' \big)^{1/2 + \textrm {i} \nu } \alpha _2(z) + \mathcal{R} , \end{align}

where $\partial _k z^c = - U(z^c) / [k U_z(z^c)]$ has been used. Note that $\nu$ was a function of $z^c(k, \omega _n + \textrm {i} r)$ in (4.3), but it is now a function of $z$ . This replacement yields an error $\propto (k - k')^{3/2 \pm \textrm {i} \nu } \log (k - k')$ , which is negligible in the following result and hence included in the residual term $\mathcal{R}$ . For $U(z) \gt 0$ , the branch point is located in the upper half of the complex $k$ plane. Its contribution to the integral at $x \gg 1$ is evaluated as

(4.11a) \begin{align} \lim _{r \to 0} \frac {1}{2 \pi } \int _{k \sim k'} \hat {\psi }_{1n} e^{\textrm {i} kx} \,\text{d}k & = \left [ \left ( \frac {U^2}{\omega _n U_z} \right )^{1/2 - \textrm {i} \nu } \frac {e^{- \pi \nu / 2} \alpha _1 x^{\textrm {i} \nu }}{\varGamma (- 1/2 + \textrm {i} \nu )}\right.\nonumber \\ &\quad+ \left.\left ( \frac {U^2}{\omega _n U_z} \right )^{1/2 + \textrm {i} \nu } \frac {e^{\pi \nu / 2} \alpha ^\dagger _1 x^{- \textrm {i} \nu }}{\varGamma (- 1/2 - \textrm {i} \nu )} \right ] \nonumber \\ &\quad\times (U_0 - U) \hat {h} \left ( \frac {\omega _n}{U} \right ) J_n \left ( \frac {n U_T}{U} \right ) \frac {\textrm {i}^{- 1/2} x^{- 3/2} e^{\textrm {i} \omega _n x / U}}{D(\omega _n / U, \omega _n)} + o\big ( x^{-3/2} \big ) ,\end{align}

(4.11b) \begin{align} \lim _{r \to 0} \frac {1}{2 \pi } \int _{k \sim k'} \hat {\rho }_{1n} e^{\textrm {i} kx} \,\text{d}k & = \left [ \left ( \frac {U^2}{\omega _n U_z} \right )^{1/2 - \textrm {i} \nu } \frac {e^{- \pi \nu / 2} \alpha _1 x^{\textrm {i} \nu }}{\varGamma (1/2 + \textrm {i} \nu )} \right.\nonumber \\ & \quad+ \left.\left ( \frac {U^2}{\omega _n U_z} \right )^{1/2 + \textrm {i} \nu } \frac {e^{\pi \nu / 2} \alpha _1^\dagger x^{- \textrm {i} \nu }}{\varGamma (1/2 - \textrm {i} \nu )} \right ] \nonumber \\ &\quad\times \frac { (U_0 - U) \omega _n N^2}{U^2} \hat {h} \left ( \frac {\omega _n}{U} \right ) J_n \left ( \frac {n U_T}{U} \right )\nonumber \\ &\quad\times\frac {\textrm {i}^{1/2} x^{- 1/2} e^{\textrm {i} \omega _n x / U}}{D(\omega _n / U, \omega _n)} + o\big ( x^{-1/2} \big ) ,\\[6pt]\nonumber \end{align}

where an asymptotic formula in Appendix C is used. Note that the exponents on $x$ in (4.11) match those reported by Camassa & Viotti (Reference Camassa and Viotti2013), who dealt with time-periodic but unforced disturbances in a stratified shear flow. The present study demonstrates that a forced problem solved with a different asymptotic approach yields consistent results.

A marked difference of (4.11) from the standing-mode wave solutions (4.6) is found in their spatial structure. The present solutions exhibit algebraic dependence on $x$ . Furthermore, owing to the factor $U(z)$ , the solution is not separable into the $x$ - and $z$ -dependent parts. As a result, $\psi$ , $\psi _z$ and $\psi _{zz}$ obey different power laws. The perturbations of velocity, density and vertical shear exhibit the dependence on $x \gg 1$ roughly as

(4.12) \begin{align} \vert w_1 \vert \propto x^{- 3/2} , \quad \vert u_1 \vert , \vert \rho _1 \vert \propto x^{- 1/2}, \quad \left \vert \frac {\partial u_1}{\partial z} \right \vert \propto x^{1/2} . \end{align}

While velocity and density decay, the velocity gradient grows far from the topography. Figure 6 illustrates the asymptotic formula (4.11a ). Notably, the structure of the solution differs from that in the single-wavenumber case. Now, the equiphase curves are almost straight and there does not exist a level where the solution is discontinuous. The angle of the equiphase lines varies in the horizontal direction, reflecting the increase of the vertical wavenumber and, consequently, the growth of the velocity gradient.

Figure 6.

Illustration of the asymptotic formula (4.11a ). Snapshots of $\psi _1$ , $\partial _z \psi _1$ and $\partial _z^2 \psi _1$ are computed for $N=1$ , $U=0.5z$ , $n=1$ , $\omega _T=1$ , $U_T = 0.1$ , $\hat {h}=1$ (i.e. $h = \delta (x)$ ) and $\alpha$ presented as (B5).

The algebraic dependence of the asymptotic solutions in (4.12) can be interpreted in terms of conservation of wave activity (pseudomomentum), a quadratic invariant associated with horizontal translation symmetry. A formal definition and the corresponding conservation law are given later in § 5.2; here, we use only the qualitative implication that, once the disturbance has left the forcing region, the pseudomomentum carried by the wave field is conserved during propagation. In the far field, pseudomomentum is dominated by the product of the density perturbation and the vorticity, $\rho _1 {\nabla} ^2 \psi _1/N^2$ (up to subdominant correction terms), so if $\rho _1$ decays, it must be compensated by the growth of ${\nabla} ^2 \psi _1$ . Because $w_1$ decays much faster than $u_1$ , the term $w_{1x}$ in ${\nabla} ^2 \psi _1=-u_{1z}+w_{1x}$ is subdominant and, hence, ${\nabla} ^2 \psi _1 \approx -u_{1z}$ . Consequently, pseudomomentum conservation suggests $\rho _1 u_{1z}=\mathcal{O}(1)$ as $x\to \infty$ , consistent with (4.12).

By contrast, wave energy $E$ and enstrophy $({\nabla} ^2 \psi _1)^2$ are not conserved in a stratified shear flow. Indeed, (4.12) tells that the energy decreases while the enstrophy increases with $x$ . Here, a difference from a single wavenumber solution is worth noting; for a singular solution (4.3), the power laws around the critical level are $\vert w_1 \vert \propto \vert z - z^c \vert ^{1/2}$ , $\vert u_1 \vert , \vert \rho _1 \vert \propto \vert z - z^c \vert ^{-1/2}$ and $\vert u_{1z} \vert \propto \vert z - z^c \vert ^{-3/2}$ . Accordingly, the wave energy density algebraically grows towards the critical level. On the other hand, for the present regular solution (4.11) composed of a continuous range of wavenumbers, the wave energy density decays during the horizontal propagation. In either case, because the velocity shear grows during propagation, waves will eventually break, resulting in energy dissipation.

In addition to $k = (\omega _n + \textrm {i} r) / U(z)$ , the integrand of the Fourier integral involves other branch points at $k = (\omega _n + \textrm {i} r) / U(0)$ and $k = (\omega _n + \textrm {i} r) / U(1)$ . Differing from the above-mentioned result, their contributions for a large $\vert x \vert$ are separable with respect to $x$ and $z$ . This study does not discuss the contributions to the solutions from these points.

5.1. Total energy budgets

We consider the fully nonlinear model (2.5) applied with the non-dimensionalisation in § 2.2. For this system, the total amounts of kinetic energy and potential energy in the frame moving with the oscillating barotropic flow are represented by

(5.1a) \begin{align} \mathcal{E}_K & = \iint _{\epsilon h}^1 \left [ \frac {1}{2}\left ( U - \frac {\partial \psi }{\partial z} \right )^2 + \frac {1}{2} \left ( \frac {\partial \psi }{\partial x} \right )^2 \right ] \,\text{d}z \,\text{d}x \end{align}

(5.1b) \begin{align} \text{and} \quad \mathcal{E}_{P} & = \iint _{\epsilon h}^1 z \rho _{tot} \,\text{d}z \,\text{d}x ,\\[6pt]\nonumber \end{align}

respectively, where the lower bound of the vertical integration $\epsilon h$ depends on $x$ and $t$ , and a sufficiently large extent is taken in the horizontal direction to cover the whole range in which the disturbance exists. The total density is constituted of the laminar reference part and the disturbances as $\rho _{tot} = \rho _0(z) + \rho (x,z,t)$ with $\rho _{0z} = - N^2$ understood.

We shall expand the total energy equations (5.1) in terms of $\epsilon$ . First, the integrand of the kinetic energy is expanded to obtain

(5.2) \begin{align} \mathcal{E}_K = \iint _{\epsilon h}^1 \frac {U^2}{2} \,\text{d}z \,\text{d}x + \epsilon \iint _{\epsilon h}^1 U u_1 \,\text{d}z \,\text{d}x + \epsilon ^2 \iint _{\epsilon h}^1 \left ( \frac {u_1^2 + w_1^2}{2} + U u_2 \right ) \,\text{d}z \,\text{d}x + \mathcal{O}(\epsilon ^3) , \end{align}

where $u_1 = - \partial _z \psi _1$ , $w_1 = \partial _x \psi _1$ and $u_2 = - \partial _z \psi _2$ are introduced. To remove the dependence on $\epsilon$ of the integration interval, we use a formula valid for an arbitrary analytic function $f(z)$ ,

(5.3) \begin{align} \int _{\epsilon h}^1 f \,\text{d}z = \int _0^1 \left [ f - \delta (z) \left ( \epsilon h f + \frac {\epsilon ^2 h^2 f_z}{2} + \cdots \right ) \right ] \,\text{d}z , \end{align}

where $\delta$ is the Dirac delta function. Consequently, the kinetic energy (5.2) is represented as the integration over a rectified domain as

(5.4) \begin{align} \mathcal{E}_K & = \int _0^1 \frac {U^2 \left \langle 1 - \epsilon \delta (z) h \right \rangle }{2} \,\text{d}z + \epsilon \int _0^1 U \left \langle u_1 \right \rangle \,\text{d}z \nonumber \\ &\quad + \epsilon ^2 \int _0^1 \left ( \left \langle \frac {u_1^2 + w_1^2}{2} \right \rangle + U \left \langle u_2^\ast \right \rangle \right ) \,\text{d}z + \mathcal{O}(\epsilon ^3) , \end{align}

in which the horizontal integration is denoted by $\left \langle \ \right \rangle \equiv \int \,\text{d}x$ and the second-order horizontal velocity is redefined as

(5.5) \begin{align} u_2^\ast \equiv u_2 - \delta (z) \left ( h u_1 + \frac {h ^2 U_z}{2} \right ) . \end{align}

Note that the additional terms represent the volume transport associated with the corrugation in the bottom boundary. If $\left \langle h \right \rangle$ were to vanish, the horizontal average of these terms would correspond to the so-called bolus velocity – a residual transport induced by the correlation between velocity and fluid thickness (this kind of transport for topography-generated internal waves is hinted in the footnote on p. 135 of Bühler (Reference Bühler2014)).

The first term on the right-hand side of (5.4) does not depend on time and is unimportant. To evaluate the second term, we shall integrate (2.11a ) over $x$ to obtain $\partial _t \partial _z \left \langle u_1 \right \rangle = 0$ , i.e. $\partial _z \left \langle u_1 \right \rangle$ is constant in time. It is natural to choose the initial condition such that $\partial _z \left \langle u_1 \right \rangle = 0$ , which yields, from (2.12), $\left \langle u_1 \right \rangle = [ U(0) - U_T \cos \omega _T t ] \left \langle h \right \rangle$ . We thus understand that the second term on the right-hand side of (5.4) is also determined by the reference state and irrelevant to the wave dynamics. It is enough to concentrate our attention on the remaining $\mathcal{O}(\epsilon ^2)$ terms.

In physical oceanography, it is widely known (e.g. Vallis Reference Vallis2017) that the potential energy in (5.1b ) can be separated into a constant part and a varying part, referred to as the background potential energy and the available potential energy, respectively. We present its detailed derivation in Appendix D and write the final result here as

(5.6) \begin{align} \mathcal{E}_{P} = \mathcal{E}_{BP} + \epsilon ^2 \int _0^1 \left \langle \frac {\rho ^2_1}{2 N^2} \right \rangle \,\text{d}z + \mathcal{O}( \epsilon ^3 ) , \end{align}

where $\mathcal{E}_{BP}$ is a constant.

Combining (5.2), (5.4) and (5.6), we represent the total amount of energy in a single expression as

(5.7) \begin{align} \mathcal{E}_K + \mathcal{E}_P = \mathcal{E}_B + \epsilon ^2 \int _0^1 \left ( \left \langle E \right \rangle + U \left \langle u_2^\ast \right \rangle \right ) \,\text{d}z + \mathcal{O}(\epsilon ^3) , \end{align}

where $\mathcal{E}_B$ is an unimportant part and can be ignored. In the remaining terms, we have defined

(5.8) \begin{align} E \equiv \frac {u_1^2 + w_1^2}{2} + \frac {\rho _1^2}{2N^2} , \end{align}

which is a quadratic function of the $\epsilon$ -order disturbance amplitude and the common definition of the energy density of internal gravity waves. The other part in (5.7), $\int _0^1 U \big \langle u_2^\ast \big \rangle \,\text{d}z$ , belongs to the energy contained in the horizontal mean flow. Variation in this energy can be regarded as feedback to the mean flow from the deviations. Importantly, wave energy $E$ is not conserved in shear flow, but, as follows from its definition, the summation of the two terms of disturbance energy, $\left \langle E \right \rangle + U \big \langle u_2^\ast \big \rangle$ , should be conserved in the absence of an external force.

We shall investigate the energy variations in the system induced by the topographic forcing. First, multiplying $-\psi _1$ to (2.11a ) and $\rho _1 / N^2$ to (2.11b ), combining them, and performing the horizontal integration, we derive

(5.9) \begin{align} \frac {\partial \left \langle E \right \rangle }{\partial t} + \frac {\partial }{\partial z} \left \langle - \psi _1 \frac {\partial ^2 \psi _1}{\partial t \partial z} \right \rangle = U \frac {\partial \left \langle u_1 w_1 \right \rangle }{\partial z} . \end{align}

In this expression, $\left \langle u_1 w_1 \right \rangle$ appearing on the right-hand side stands for the Reynolds stress, i.e. the vertical flux of the horizontal momentum. Through this term, waves exchange energy with the horizontal mean flow. Next, performing the horizontal integration to the horizontal part of equations in (2.5a ) and collecting the $\epsilon ^2$ -order terms, we derive an equation governing the variations in the mean flow as

(5.10) \begin{align} \frac {\partial \left \langle u_2 \right \rangle }{\partial t} + \frac {\partial \left \langle u_1 w_1 \right \rangle }{\partial z} = P . \end{align}

Here, $P$ denotes the difference in the pressure $p_2$ between $x = - \infty$ and $x = \infty$ , and is uniform in the vertical direction; because the fluid is motionless sufficiently far from the localised topography, a hydrostatic balance is maintained there and the pressure deviation from the reference state does not depend on the depth. This pressure force is essential to fulfil the constraint of the net volume transport. That is, expanding (2.7) in terms of $\epsilon$ and collecting the second-order terms, we derive $\int u_2^\ast \,\text{d}z = 0$ , which combined with (5.10) yields

(5.11) \begin{align} \left ( \left \langle u_1 w_1 \right \rangle - \frac{\text{d}}{\text{d}t} \left \langle h u_1\right \rangle \right )_{z=0} + P = 0 , \end{align}

where we have used ${\rm d} \big \langle h^2 \big \rangle / \text{d}t = 0$ . This expression represents the horizontal momentum balance in the whole system; the first two terms represent the form stress exerted by the bottom topography, which is compensated by the ambient pressure gradient to keep the net volume transport constant.

Even though the barotropic transport of the mean flow is prescribed by the reference state, the baroclinic part of the mean flow can be accelerated/decelerated by the convergence of the Reynolds stress as (5.10), and its reaction performs work on wave motion through the right-hand side of (5.9). Combining the two equations in (5.9) and (5.10), we may describe the energy balance at each depth as

(5.12) \begin{align} \frac {\partial \left ( \left \langle E \right \rangle + U \left \langle u_2 \right \rangle \right )}{\partial t} + \frac {\partial }{\partial z} \left \langle - \psi _1 \frac {\partial ^2 \psi _1}{\partial t \partial z} \right \rangle & = PU . \end{align}

The physical interpretation of the second term on the left-hand side is drawn from the horizontal momentum equation, $- \partial _t \partial _z \psi _1 = \partial _t u_1 = - U \partial _x u_1 - U_z w_1 - \partial _x p_1$ , and the integration by parts; this term represents the divergence of energy flux composed of pressure work and kinetic energy advection, i.e.

(5.13) \begin{align} \left \langle - \psi _1 \frac {\partial ^2 \psi _1}{\partial t \partial z} \right \rangle = \left \langle p_1 w_1 \right \rangle + \left \langle U u_1 w_1 \right \rangle . \end{align}

Vertically integrating (5.12), using (5.11) and taking into account the correction term at the bottom in (5.5), we write down the net production rate of disturbance energy as

(5.14) \begin{align} \frac{\text{d}}{\text{d}t} \int _0^1 \left ( \left \langle E \right \rangle + U \left \langle u^\ast _2 \right \rangle \right ) \,\text{d}z = \mathcal{P}^{(T)} + \mathcal{P}^{(S)} + \frac{\text{d}}{\text{d}t} \left ( \left ( U - \overline {U} \right ) \left \langle h \frac {\partial \psi _1}{\partial z} \right \rangle \right )_{z = 0} \equiv \mathcal{P}(t) \end{align}

with

(5.15) \begin{align} \mathcal{P}^{(T)} \equiv \left \langle - \psi _1 \frac {\partial ^2 \psi _1}{\partial t \partial z} \right \rangle _{z = 0} \quad \text{and} \quad \mathcal{P}^{(S)} \equiv - \overline {U} \left \langle u_1 w_1 \right \rangle _{z=0} , \end{align}

where we have written the barotropic velocity as $\overline {U} \equiv \int _0^1 U \,\text{d}z$ . The third term in the middle expression of (5.14) represents the variations in the kinetic energy associated with the volume transport on the varying boundary. This term owes its existence to the vertical gradient of the steady flow; if $U$ is uniform, $U_{z=0} - \overline {U}$ identically vanishes. In this form, the difference between the barotropic velocity and the velocity at the bottom demands a minor modification in the energy budget of the whole system. In any case, a long-term average of this term inevitably vanishes for a solution whose amplitude is bounded all the time. Therefore, volume transport induced by the corrugated topography makes no contributions to the net energy conversion from barotropic to baroclinic components. In the remaining terms of (5.14), as elucidated in (5.13), $\mathcal{P}^{(T)}$ stands for the vertical energy flux from the bottom and accounts for the energy supply from the oscillatory barotropic tide. In contrast, $\mathcal{P}^{(S)}$ corresponds to the pressure work done by the steady flow.

5.2. Pseudomomentum and pseudoenergy

As we have seen, the second-order horizontal velocity, $u^\ast _2$ , plays a key role in the energetics. Unlike the other first-order quantities, this variable is difficult to directly obtain. Nonetheless, its horizontal integration, $\big \langle u^\ast _2 \big \rangle$ , can be diagnosed based on the other variables by use of the concept of pseudomomentum. Adding the bottom-correction term to that formulated by Shepherd (Reference Shepherd1990), let us define the pseudomomentum density (the present definition of $\textsf {p}$ involves only Eulerian variables and it differs from the pseudomomentum density in the generalised Lagrangian mean theory. For a plane wave limit, averaged over a phase, the two definitions will coincide (Bühler Reference Bühler2014)) as

(5.16) \begin{align} {\textsf{p}} = \frac {\rho _1 {\nabla} ^2 \psi _1}{N^2} - \frac {U_{zz} \rho _1^2}{2 N^4} - \delta (z) \left ( h u_1 + \frac {h ^2 U_z}{2} \right ) , \end{align}

which is a quadratic quantity in terms of $(h, \psi _1, \rho _1)$ . The horizontal integration of $\textsf {p}$ follows a similar form of the equation as for $\big \langle u_2^\ast \big \rangle$ ; specifically, from (2.11),

(5.17) \begin{align} \frac {\partial\!\left \langle \textsf {p} \right \rangle }{\partial t} + \frac {\partial\!\left \langle u_1 w_1 \right \rangle }{\partial z} = 0 \end{align}

for $z \neq 0$ is established. Comparing this to (5.10), we learn that the two equations differ only by the vertically uniform pressure force $P$ . Given that $P$ is determined from a constraint $\int _0^1 u_2^\ast \,\text{d}z = 0$ , we come up with a relevant quantity for the pseudomomentum by projecting $\textsf {p}$ to the baroclinic component, i.e. $\textsf {p} - \int _0^1 \textsf {p} \,\text{d}z \equiv \textsf {p}'$ . Consequently, $\big \langle \textsf {p}' \big \rangle$ is shown to obey the same equation as $\big \langle u_2^\ast \big \rangle$ . Therefore, neglecting the difference in their initial conditions, the $\epsilon ^2$ -order modification in the horizontal mean velocity, $\big \langle u_2^\ast \big \rangle$ , is identified as the baroclinic part of the pseudomomentum, $\big \langle \textsf {p}' \big \rangle$ .

In addition to the pseudomomentum density, it is also informative to define the two kinds of pseudoenergy density as $\textsf {e} = E + U \textsf {p}$ and $\textsf {e}' = E + U \textsf {p}'$ . Note that the original definitions $\textsf {p}$ and $\textsf {e}$ are the locally conserved quantities – obeying (5.17) and

(5.18) \begin{align} \frac {\partial\!\left \langle \textsf {e} \right \rangle }{\partial t} + \frac {\partial }{\partial z} \left \langle - \psi _1 \frac {\partial ^2 \psi _1}{\partial t \partial z} \right \rangle = 0 \end{align}

for $z \neq 0$ . However, to discuss the production rate of the disturbance energy, the modified form of pseudoenergy $\textsf {e}'$ is more relevant because its volume integral coincides with the $\epsilon ^2$ -order disturbance energy defined in (5.7). It is enlightening to write the equivalence of three forms of energy representation,

(5.19) \begin{align} \int _0^1 \big \langle \textsf {e}' \big \rangle \,\text{d}z = \int _0^1 \left \langle E \right \rangle \,\text{d}z + \int _0^1 U \left \langle u_2^\ast \right \rangle \,\text{d}z = \int _0^1 \left \langle \textsf {e} \right \rangle \,\text{d}z - \overline {U} \int _0^1 \left \langle \textsf {p} \right \rangle \,\text{d}z . \end{align}

The final expression as well as (5.17) and (5.18) provide another physical interpretation for (5.14) and (5.15); $\mathcal{P}^{(T)}$ and $\mathcal{P}^{(S)}$ stand for the vertical flux of pseudoenergy and pseudomomentum multiplied by minus the barotropic velocity, both of which are injected at the bottom, respectively. The disturbance energy generated on the bottom at the rate of $\mathcal{P}$ is partitioned into the wave part, $\int _0^1 \left \langle E \right \rangle \,\text{d}z$ , and the mean flow part, $\int _0^1 U \big \langle \textsf {p}' \big \rangle \,\text{d}z$ . The mean flow is then accelerated at a rate of $\partial _t \big \langle u_2^\ast \big \rangle = \partial _t \big \langle \textsf {p}' \big \rangle$ at each depth via the convergence of the Reynolds stress and the barotropic pressure force. Even though the volume integral of $\textsf {p}'$ is identically $0$ , this quantity plays the role of vertical redistribution of momentum through the Reynolds stresses.

Before concluding this subsection, we highlight a marked contrast in energetics between two typical scenarios of topographic wave generation with distinct horizontal boundary conditions: the unbounded system analysed in this study versus a system with periodic boundary conditions. In the atmospheric context, the periodic conditions are more relevant, as the primary focus is often on the interactions between bottom-generated waves and a planetary-scale zonal mean flow. Under such conditions, the pressure gradient integrated across the zonal direction vanishes, leading to variations in the total volume transport in response to momentum injection at the bottom boundary. The acceleration rate of the mean flow at each altitude can be evaluated using the temporal evolution of the zonally integrated pseudomomentum, $\left \langle \textsf {p} \right \rangle$ . Similarly, variations in the total energy correspond directly to changes in the pseudoenergy, $\int \left \langle \textsf {e} \right \rangle \,\text{d}z$ . However, if the bottom topography is stationary (i.e. in the absence of oscillatory barotropic forcing), the amount of pseudoenergy for bottom-generated waves is null. It is because, in such a case, the energy conversion from the steady current to wave motion via topographic interactions represents an internal redistribution of energy within the coupled wave-mean flow system. Consequently, this conversion is not reflected in the pseudoenergy production rate. This outcome is a well-established result in lee wave generation theory (e.g. Bühler Reference Bühler2014, p. 136). In contrast, in a horizontally unbounded system, as considered in this study, the energetic framework fundamentally differs due to the presence of a net pressure force acting far from the bottom topography. That is to say, a steady flow can sustain the generation of internal waves without losing its energy, supported by an external pressure force applied to the system.

5.3. Energy conversion rate for a stationary state

We now examine the disturbance energy production rate, as specified by (5.14) and (5.15), for a stationary state considered in § 4. Given that all variables are periodic with period $T \equiv 2 \pi / \omega _T$ , the net energy production rate is quantified by averaging over time. The disturbance motions conceptually involve both the barotropic and baroclinic components. However, since the barotropic velocity is prescribed by the reference state, its energy remains bounded. Consequently, the time-averaged energy production rate essentially represents the conversion rate from the barotropic energy in the reference state to the baroclinic energy in disturbances.

To inspect the contribution from the steady flow part $\mathcal{P}^{(S)}$ , we first focus on the topographic form stress, which is equivalent to the pseudomomentum input or minus the ambient pressure force, as shown in (5.11). Using Parseval’s theorem, we may write the time-averaged form stress as

(5.20) \begin{align} \frac {1}{T}\int _0^T (-P) \,\text{d}t = \frac {1}{T} \int _0^T \left \langle u_1 w_1 \right \rangle _{z=0}\, \text{d}t = \sum _{n = -\infty }^\infty\!\tau _n . \end{align}

An explicit expression for the form stress of the $n$ th harmonic component is found by inserting the modal expansion (3.1) and (4.1) in $\left \langle u_1 w_1 \right \rangle$ with a condition $\hat {\psi }_{1n}(z; k, r) = \hat {\psi }^\dagger _{1-n}(z; - k, r)$ , leading to

(5.21) \begin{align} \tau _n & = \lim _{r \to 0} \Re \left [ \frac {\textrm {i}}{2\pi } \int k \hat {\psi }^\dagger _{1n} \frac {\partial \hat {\psi }_{1n}}{\partial z}\, {\rm d} k \right ]_{z=0} \nonumber \\ & = \lim _{r \to 0} \frac {1}{2 \pi } \int k \big \vert \hat {\psi }_{1n}(0; k) \big \vert ^2 \Im \left [ \frac {- A_{1z}(0; k, \omega _n + \textrm {i} r)}{D(k, \omega _n + \textrm {i} r)} \right ] \,\text{d}k , \end{align}

where the second equality stems from the definition of $\hat {\psi }_{1n}$ in (4.2), its bottom condition is $\hat {\psi }_{1n}(0; k) = (k U_0 - \omega _n) \hat {h}(k) J_n(k U_T / \omega _T) / k$ , $A_1$ is defined in (A4) and $D(k, \sigma )= A_1(0; k, \sigma )$ is understood. The integrand is finite when $\omega _n$ belongs to the continuous part of the spectrum $\varSigma ^c_k$ and we shall invert this condition to write $k \in K^c_{\omega _n}$ . For $k \notin K^c_{\omega _n}$ , since $A_1$ is real in the limit of $r \to 0$ , the integrand vanishes almost everywhere. The exceptional points are the zeros of $D$ , i.e. the discrete parts of the spectrum, $K^d_{\omega _n}$ . Using (4.5a ), contributions from these points can be evaluated by using the formula

(5.22) \begin{align} \lim _{r \to 0} \Im \frac {- 1}{2 \pi D(k, \sigma + \textrm {i} r)} = \sum _j \frac {\delta (k - k_j)}{2 \big \vert c_j^g (\sigma ) \big \vert D_\sigma (k, \sigma )} , \end{align}

which is valid for $k \notin K^c_\sigma$ . Here, we have used

(5.23) \begin{align} D(k, \sigma + \textrm {i} r) & = D(k, \sigma ) + \textrm {i} r D_\sigma (k, \sigma ) + \mathcal{O}(r^2) \nonumber \\ & = (k - k_j) F(k, \sigma ) + \textrm {i} r D_\sigma (k, \sigma ) + \mathcal{O}(r^2) \end{align}

with $F = D_k$ at $k = k_j$ understood and a formula of a nascent delta function

(5.24) \begin{align} \lim _{r\rightarrow \pm 0} \frac {r}{\pi (x^2+r^2)} = \pm \delta (x) . \end{align}

The form stress for the $n$ th harmonic component (5.21) is consequently represented as

(5.25) \begin{align} \tau _n & = \sum _{k_j \in K^d_{\omega _n}} k_j \chi _j (\omega _n) \big \vert \hat {\psi }_{1n}(0; k_j) \big \vert ^2 + \frac {1}{2 \pi } \int _{K^c_{\omega _n}} k R(k, \omega _n) \big \vert \hat {\psi }_{1n}(0; k) \big \vert ^2 \,\text{d}k \end{align}

with

(5.26) \begin{align} \chi _j (\sigma ) \equiv \frac {A_{1z}(0; k_j, \sigma )}{2 \big \vert c_j^g(\sigma ) \big \vert D_\sigma (k_j, \sigma )} \quad \text{and} \quad R(k, \sigma ) \equiv - \Im \left [ \frac { A_{1z}(0; k, \sigma )}{D(k, \sigma )} \right ] . \end{align}

The two functions $\chi _j$ and $R$ play the roles of transforming a prescribed boundary condition of $\psi$ at $z=0$ to the vertical wave action flux at the bottom. Here, wave action is regarded as the pseudomomentum divided by the wavenumber. A more common definition of the wave action would be the wave energy divided by the intrinsic frequency – wave frequency observed in a frame moving with a background current. Indeed, this classical definition also applies to the present case, but care should be taken that the background current is not uniform. This point will be discussed in § 6. We shall call $\chi _j$ and $R$ the action production coefficients for the discrete and continuous spectra, respectively.

With the harmonic decomposition of the form stress in hand, the time-averaged energy conversion rate from the steady background flow to the baroclinic motion follows directly as $(1/T) \int _0^T \mathcal{P}^{(S)} \, \text{d}t = -\overline {U} \sum _{n = -\infty }^\infty \tau _n$ . Similarly, the time-averaged energy conversion rate from the barotropic tide to the baroclinic motion, or equivalently the pseudoenergy production rate, is written as $(1/T) \int _0^T \mathcal{P}^{(T)} \, \text{d}t = \sum _{n = -\infty }^\infty \mathcal{P}^{(T)}_n$ with

(5.27) \begin{align} \mathcal{P}^{(T)}_n = \sum _{k_j \in K^d_{\omega _n}} \omega _n \chi _j (\omega _n) \big \vert \hat {\psi }_{1n}(0; k_j) \big \vert ^2 + \frac {1}{2 \pi } \int _{K^c_{\omega _n}} \omega _n R(k, \omega _n) \big \vert \hat {\psi }_{1n}(0; k) \big \vert ^2 \,\text{d}k . \end{align}

This expression makes clear that the pseudoenergy flux is given by the wave frequency multiplied by the wave-action flux, consistent with the standard pseudoenergy–action relation. In particular, the pseudoenergy production rate of the zero-frequency component vanishes identically, $\mathcal{P}^{(T)}_0 = 0$ ; for this mode, the baroclinic motion is driven solely by the steady background flow and the oscillatory tidal forcing does not contribute energetically, a general consequence of the pseudoenergy-based framework (Bühler Reference Bühler2014). Note that the form stress can remain finite even for the zero-frequency component.

Combining the steady and tidal contributions, we obtain the total time-averaged barotropic-to-baroclinic energy conversion rate, $(1/T) \int _0^T \mathcal{P} \, \text{d}t = \sum _{n = -\infty }^\infty \mathcal{P}_n$ , with

(5.28) \begin{align} \mathcal{P}_n & = \sum _{k_j \in K^d_{\omega _n}} \left ( \omega _n - \overline {U} k_j \right ) \chi _j (\omega _n) \big \vert \hat {\psi }_{1n}(0; k_j) \big \vert ^2 \nonumber \\ & + \frac {1}{2 \pi } \int _{K^c_{\omega _n}} \left ( \omega _n - \overline {U} k \right ) R(k, \omega _n) \big \vert \hat {\psi }_{1n}(0; k) \big \vert ^2 \,\text{d}k . \end{align}

By analogy with the pseudomomentum or pseudoenergy flux, one might be tempted to interpret (5.28) as a wave-energy flux written as an intrinsic frequency times an action flux. However, this interpretation is generally incorrect because $\omega _n - \overline {U} k$ is not an intrinsic frequency unless $U(z)$ is uniform. It is important to remember that $(1/T)\int _0^T \mathcal{P} \, \text{d}t$ represents the barotropic-to-baroclinic energy conversion rate, not the wave-energy production rate.

Returning to (5.26), the discrete-spectrum coefficient $\chi _j$ depends on the group velocity $c_j^g$ . If $c_j^g$ happens to vanish, $\chi _j$ diverges, and the form stress and the conversion rate become unbounded within inviscid linear theory. Aside from this special case, the time-averaged conversion rate is finite and time independent. In an unbounded domain, the injected disturbance energy is carried away by radiating motions towards $x=\pm \infty$ ; consequently, the total disturbance energy integrated over the domain grows linearly in time.

This constant-rate growth can be interpreted as a quasi-resonant response of a spatially unbounded system to spatially localised forcing. As discussed in § 3.4, for each fixed horizontal wavenumber $k$ , the Taylor–Goldstein equation admits a discrete set of regular eigenfrequencies, say $\sigma =\sigma _j(k)$ . As $k$ varies, these eigenfrequencies depend continuously on $k$ , forming dispersion curves in the $(k,\sigma )$ plane. Because localised topography excites a continuum of horizontal wavenumbers, a given forcing frequency $\omega _n$ can interact with these dispersion curves at isolated wavenumbers $k_j$ satisfying $D(k_j,\omega _n)=0$ (i.e. $k_j\in K^d_{\omega _n}$ ). When the intersection is regular (non-vanishing group velocity), each contribution remains finite and the forcing supplies energy at a constant rate. Colin de Verdiëre & Saint-Raymond (Reference Colin De Verdiëre and Saint-Raymond2020) referred to this mechanism as quasi-resonance and contrasted it with a genuine resonance, in which the energy grows faster than linearly so that a time-independent conversion rate is not defined. The divergence of $\chi _j$ at $c_j^g=0$ corresponds to the onset of this resonant behaviour.

To illustrate the conversion-rate formula (5.28), figure 7 shows numerical evaluations for a simple background state. We take uniform stratification $N=1$ and a linear shear flow $U(z)=0.05+0.15z$ , together with a barotropic tide of amplitude $U_T=0.2$ and frequency $\omega _T=0.2$ . The topography is chosen as an extremely narrow ridge; accordingly, its Fourier transform is nearly flat and, for simplicity, we set $\hat {h}(k)=1$ in the computation. We evaluate separately the contributions associated with the tidal and steady-current parts, and further split each into the discrete- and continuous-spectrum contributions. Harmonics up to $n=2$ are retained ( $n=0,1,2$ ). For the discrete spectrum, we also separate left- and right-propagating contributions according to the sign of the horizontal group velocity $c^g$ . In this example, a substantial contribution comes from the discrete modes of the zero-frequency component ( $n=0$ ), i.e. lee waves. For the tide-frequency component ( $n=1$ ), the upstream-propagating branch ( $c^g\lt 0$ ) is stronger than the downstream-propagating branch ( $c^g\gt 0$ ), consistent with the numerical simulations of Lamb & Dunphy (Reference Lamb and Dunphy2018).

Figure 7.

Numerically computed energy conversion rates for the parameter set described in the text ( $N=1$ , $U(z)=0.05+0.15z$ , $U_T=0.2$ , $\omega _T=0.2$ and $\hat {h}(k)=1$ ). (a) Spectrum of the operator $\unicode{x1D648}_k$ : red and blue curves show the discrete spectrum (dispersion curves of regular eigenmodes), while the magenta shading indicates the continuous spectrum. The dashed line corresponds to $\sigma /k=\overline {U}=0.125$ (phase speed equal to the barotropic velocity). Horizontal dotted lines indicate the forcing frequencies $\sigma = \omega _n$ for $n=1,2$ ; the blue dot marks an upstream-propagating intersection with the discrete spectrum ( $c^g\lt 0$ ). (b) Integrand of the continuous-spectrum contribution in (5.28) for $n=1$ , separated into the tidal part, the steady-current part and their sum. The vertical dashed line indicates $k=\omega _1/\overline {U}$ . (c) Contributions to the conversion rates from the discrete and continuous spectra for $n=0,1,2$ . Discrete-spectrum contributions are further divided into left- ( $c^g\lt 0$ ) and right-propagating ( $c^g\gt 0$ ) components. For $n=1,2$ , the plotted values are multiplied by 2 to include the equal contributions from the negative harmonics $n=-1, -2$ (i.e. $\mathcal{P}_{-n}=\mathcal{P}_n$ ).

Although the net contribution from the continuous spectrum is small in this example, this does not imply that the critical-level (singular) response is negligible. Rather, the tidal and steady-current parts contribute comparably, but with opposite signs, leading to a strong cancellation. This is illustrated in figure 7(b), which plots the integrand of the continuous-spectrum contribution to (5.28) for $n=1$ . Over the integration interval, the tidal part (proportional to $\omega _n$ ) is positive whereas the steady-current part (proportional to $-\overline {U}k$ ) is negative, because $R(k,\omega _n)\gt 0$ . Their magnitudes cross at $k=\omega _n/\overline {U}$ (vertical dashed line), where the sum changes sign. For the second harmonic ( $n=2$ ), the discrete-spectrum contribution is smaller and the continuous-spectrum contribution becomes more prominent; although the two parts still oppose each other, their individual magnitudes are larger than in the $n=1$ case. This trend suggests that continuous-spectrum effects become increasingly important for higher-frequency forcing. Notably, in classical no-shear formulations (e.g. Khatiwala Reference Khatiwala2003), the spectrum is purely discrete and no continuous-spectrum (critical-layer) contribution arises. The present results suggest that this discrete-mode-only approximation can miss an important part of the conversion budget when background shear is present, especially at higher harmonics.

In the previous computations, $A_1$ and $D$ in (A4) and (A7) were evaluated using the closed-form expressions derived in Appendix B. For the continuous spectrum, the calculation can be streamlined by expressing $A_1(z;k,\sigma )$ in terms of a pair of linearly independent solutions $\gamma ^\pm (z;k,\sigma )$ of the Taylor–Goldstein equation (Appendix A). When the Richardson number at the critical level is large, $N^2/U_z^2\gg 1/4$ at $z=z^c$ , the connection rule (A10) implies $|\gamma ^+(0)|\gg |\gamma ^-(0)|$ , whereas $|\gamma ^+(1)|=|\gamma ^-(1)|$ . Accordingly, the action production coefficient can be approximated by

(5.29) \begin{align} R(k,\sigma ) \approx -\Im \left [\frac {\gamma ^+_z(0)}{\gamma ^+(0)}\right ] \equiv R_\infty (k,\sigma ), \end{align}

which is essentially insensitive to the upper boundary condition. Physically, when the Richardson number at the critical level is large, incident waves are mostly absorbed and do not reflect back towards the bottom, so the problem approaches the classical unbounded setting. Indeed, setting $U_0=0$ and applying a WKB approximation near $z=0$ yields $\gamma ^+ \propto e^{\textrm {i} m z}$ with $m=-|k|\sqrt {N^2-\sigma ^2}/\sigma$ . One then finds $\sigma R_\infty \approx |k|\sqrt {N^2-\sigma ^2}$ for $|\sigma |\lt N$ , so that the continuous-spectrum term in (5.27) reduces to (5.3) of Bell (Reference Bell1975). A similar insensitivity to the upper boundary has been noted in finite-depth viscous problems when strong dissipation prevents reflection (Shakespeare & Hogg Reference Shakespeare and Hogg2017; Shakespeare, Arbic & Hogg Reference Shakespeare, Arbic and Hogg2021). In oceanographic applications, Bell’s unbounded formula is often used to estimate the generation of short-wavelength lee waves (Legg Reference Legg2021); the present analysis provides a rationale for this practice. Nevertheless, whether turbulent mixing can be parametrised based solely on such elementary lee-wave theory remains debated, as discussed in the next section.