2.1. Base state

We specify that the background velocity is in the $x$-direction, and both background velocity and stratification are steady and vary only in the vertical. The background velocity is given by $(U(z),0,0)$, pressure by $\bar {p}(y,z)$ and buoyancy by $\bar {b}(y,z)$, with the bar notation indicating a background field. Assuming that the impact of perturbations on the mean flow is not leading order, from (2.1) it must satisfy both geostrophic and hydrostatic balance

(2.4)\begin{gather} -fU ={-}\rho_0^{{-}1}\bar{p}_y, \end{gather}

(2.5)\begin{gather}0 ={-}\rho_0^{{-}1}\bar{p}_z + \bar{b}. \end{gather}

Eliminating $\bar {p}$ from (2.4)–(2.5) gives the thermal wind balance

(2.6)\begin{equation} -fU_z = \bar{b}_y. \end{equation}

Requiring that the stratification $N^{2} = \bar {b}_z$ is a function of $z$ only, (2.6) gives that $fU_{zz} = 0$. We therefore only consider base states such that $fU_{zz} = 0$, but continue the derivation for general $U(z)$ for use when $f=0$. This ensures that although $\bar {p}$ and $\bar {b}$ are functions of $y$, $\bar {p}_y$ and $\bar {b}_y$ are not, and the problem remains effectively two-dimensional so that all coefficients of the linearised problem to be derived in § 2.3 are independent of $y$.

2.2. Energy loss

A representation of lee wave energy loss is crucial to understanding the structure of the lee wave field in the vertical. Lee wave energy must either be reabsorbed by the mean flow, or lost to dissipation and mixing. The latter is a result of energy transfer to smaller scales through instabilities of the waves themselves, or through nonlinear interactions with other waves and the background flow. In our idealised linear model, we cannot properly represent either the dynamics of the waves which can lead to instabilities and breaking, or small scales from other sources of turbulence that act to eventually dissipate even linear waves. The effect of this energy lost from the lee wave field must therefore be parametrised.

Parametrisation of dissipation and mixing at the sub-grid scale in models is generally implemented through Laplacian (or higher-order) viscous and diffusive terms in the momentum and buoyancy equations – as shown in (2.1)–(2.2). Shakespeare & Hogg (Reference Shakespeare and Hogg2017) provide a comprehensive overview of the role of Laplacian viscosity and diffusivity in the linear lee wave problem, with a focus on preventing excessive dissipation in wave resolving models. Here, we do not represent the processes that drain energy from the lee wave field, so aim to model them diffusively with this parametrisation. However, unlike Shakespeare & Hogg (Reference Shakespeare and Hogg2017), we are dealing with background flows that vary in the vertical, and thus including the vertical components $\mathcal {A}_v\boldsymbol {u}_{zz}^{{{\dagger}} }$ and $\mathcal {D}_vb_{zz}^{{{\dagger}} }$ of the Laplacian terms in our study significantly complicates the solution.

For mathematical convenience, we therefore represent the total viscous and diffusive terms by the horizontal components only. This allows some scale selection for energy loss (improving on, say, a simple Rayleigh friction), without overly complicating the problem. Using only the horizontal component as a proxy for the total dissipation and mixing has certain drawbacks, including invalidating any solutions where the vertical wavelength changes drastically or becomes very small, e.g. at critical levels. It is important to keep in mind the simplifications made here when analysing the model mixing and dissipation in § 4. Direct comparisons between our horizontal turbulent viscosity $\mathcal {A}_h$ and diffusivity $\mathcal {D}_h$ parameters and other studies or models should also be made with care, since they represent both horizontal and vertical viscosity and diffusivity. Furthermore, since $\mathcal {A}_h$ and $\mathcal {D}_h$ represent both background turbulent processes and breaking of the lee wave field itself, their ‘real’ values should depend on nonlinearity of the wave field and properties of the background flow, among other things. Although the simplifications made with this parametrisation are likely to modify our solutions somewhat, we believe that the key results of this study are unaffected.

2.3. Linearisation

For $Fr_L \ll 1$, we consider small perturbations to the base state described in § 2.1. We consider a quasi-2-D flow such that the topography and waves are uniform in $y$. This choice is sufficient to demonstrate our main findings and greatly simplifies the problem, although the theory can be extended to 2-D topography without significant difficulties. The coefficients of the linearised equations are independent of $y$ due to the constraints on the base state described in § 2.1, thus the perturbation variables can be taken to be independent of $y$. We also assume here that the perturbations are steady, although this need not be imposed at this point and follows from the application of the steady boundary conditions to be described in § 2.4.

Letting $\boldsymbol {u}^{{{\dagger}} } = (U(z) + u(x,z),v(x,z),w(x,z))$, $b^{{{\dagger}} } = \bar {b}(y,z) + b(x,z)$, $p^{{{\dagger}} } = \bar {p}(y,z) + p(x,z)$ and linearising (2.1)–(2.3) gives

(2.7)\begin{gather} wU_z + Uu_x -fv ={-}\rho_0^{{-}1}p_x + \mathcal{A}_h u_{xx}, \end{gather}

(2.8)\begin{gather}Uv_x + fu = \mathcal{A}_h v_{xx}, \end{gather}

(2.9)\begin{gather}\alpha Uw_x ={-}\rho_0^{{-}1}p_z + b + \alpha \mathcal{A}_h w_{xx}, \end{gather}

(2.10)\begin{gather}Ub_x - fvU_z +wN^{2} = \mathcal{D}_h b_{xx}, \end{gather}

(2.11)\begin{gather}u_{x} + w_{z} = 0, \end{gather}

where $\alpha \in \{0,1\}$, so that when $\alpha = 0$ the equations are hydrostatic. The hydrostatic assumption is made when the ratio of vertical to horizontal scales is small, as is often the case for lee waves. We introduce a perturbation streamfunction $\psi$ such that $u = -\psi _z$, $w = \psi _x$, with Fourier transform $\hat {\psi }(k,z)$ defined such that

(2.12)\begin{equation} \psi(x,z) =\frac{1}{2{\rm \pi}} \int_{-\infty}^{\infty} \hat{\psi}(k,z)\textrm{e}^{\textrm{i}kx} \,\textrm{d}k . \end{equation}

Taking the Fourier transform in $x$ of (2.7)–(2.11) and solving for the transformed streamfunction $\hat {\psi }(k,z)$ gives a second-order ordinary differential equation

(2.13)\begin{equation} \hat{\psi}_{zz} + P(k,z)\hat{\psi}_z + Q(k,z)\hat{\psi} = 0, \end{equation}

where

(2.14)\begin{gather} P(k,z) = \frac{f^{2}U_z\left(2U - \textrm{i}k(\mathcal{A}_h + \mathcal{D}_h)\right)}{\left(k^{2}(U - \textrm{i}k\mathcal{A}_h)^{2} - f^{2}\right)\left(U -\textrm{i}k\mathcal{A}_h\right)\left(U-\textrm{i}k\mathcal{D}_h\right)}, \end{gather}

(2.15)\begin{gather}Q(k,z) = \frac{k^{2}\left(U-\textrm{i}k\mathcal{A}_h\right)\left(N^{2} - \alpha k^{2}(U-\textrm{i}k\mathcal{A}_h)(U-\textrm{i}k\mathcal{D}_h)\right)}{\left(U - \textrm{i}k\mathcal{D}_h\right)\left(k^{2}(U-\textrm{i}k\mathcal{A}_h)^{2} - f^{2}\right)} - \frac{k^{2}U_{zz}(U-\textrm{i}k\mathcal{A}_h)}{k^{2}(U-\textrm{i}k\mathcal{A}_h)^{2} - f^{2}}. \end{gather}

With constant background velocity and stratification and in the absence of viscosity and diffusivity, this reduces to the familiar equation for the steady lee wave problem (Bell Reference Bell1975)

(2.16)\begin{equation} \hat{\psi}_{zz}(k,z) + k^{2}\frac{N^{2} - \alpha U^{2}k^{2}}{U^{2}k^{2} - f^{2}}\hat{\psi}(k,z) = 0, \end{equation}

with solution

(2.17)\begin{equation} \hat{\psi}(k,z) = A(k)\textrm{e}^{\textrm{i}m(k)z} + B(k)\textrm{e}^{-\textrm{i}m(k)z}, \end{equation}

for some functions $A$ and $B$ to be specified by the boundary conditions, where

(2.18)\begin{equation} m^{2}(k) = k^{2}\frac{N^{2} - \alpha U^{2}k^{2}}{U^{2}k^{2} - f^{2}}. \end{equation}

It is clear from (2.17) and (2.18) that there are radiating solutions (lee waves) only when $m$ is real, that is when the topographic wavelength $k$ satisfies

(2.19)\begin{equation} |\,f| < |Uk| < |N|. \end{equation}

For wavenumbers $k$ in this radiating range, rotation can be neglected when $|\,f| \ll |Uk|$, and the hydrostatic assumption ($\alpha = 0$) can be made when $|Uk| \ll |N|$, since in this case the vertical wavenumber $m \sim {N}/{U}$ (from (2.18)), so $|Uk|/|N|$ represents the ratio of vertical to horizontal wavelengths.

2.4. Boundary conditions

2.4.1. Bottom boundary condition

For a given $k$, (2.13) requires two boundary conditions. A free slip condition to ensure that the flow is parallel to the 2-D topography $h(x)$ is given by

(2.20)\begin{equation} w^{{{{\dagger}}}}(x,h) = u^{{{{\dagger}}}}(x,h)h_x. \end{equation}

Linearising about the base state then gives

(2.21)\begin{equation} w(x,0) = U(0)h_x, \end{equation}

or equivalently, defining the Fourier transform of the topography $\hat {h}(k)$ similarly to (2.12)

(2.22)\begin{equation} \hat{\psi}(k,0) = U(0)\hat{h}(k). \end{equation}

Given this requirement, we write $\hat {\psi }(k,z) = U(0)\hat {h}(k)\hat {\zeta }(k,z)$, where $\hat {\zeta }(k,z)$ is the normalised vertical structure function for a wavenumber $k$, so that

(2.23)\begin{equation} \psi(x,z) =\frac{U(0)}{2{\rm \pi}} \int_{-\infty}^{\infty} \hat{\zeta}(k,z)\hat{h}(k)\textrm{e}^{\textrm{i}kx}\, \textrm{d}k, \end{equation}

and $\hat {\zeta }(k,z)$ satisfies

(2.24)\begin{gather} \hat{\zeta}_{zz} + P(k,z)\hat{\zeta}_z + Q(k,z)\hat{\zeta} = 0, \end{gather}

(2.25)\begin{gather}\hat{\zeta}(k,0) = 1. \end{gather}

2.4.2. Upper boundary condition

For the second condition, first consider the classical unbounded lee wave problem, which requires that waves propagate freely through the upper boundary. When the background flow is uniform in $z$, $P$ vanishes, $Q$ is constant in $z$ and a vertical wavenumber $m(k) = \pm \sqrt {Q(k)}$ can be found. If the flow is also inviscid, $m(k)$ is given up to a sign by (2.18), and for the wave-like solutions with $k$ in the radiating range (2.19), there is then a well-defined vertical group velocity (to be discussed further in § 2.5). The vertical group velocity must be positive when the solutions are wave like so that energy radiates away from topography, and this is ensured by choosing $m(k)$ to have the same sign as $Uk$. When $m$ is imaginary (non-wave-like solutions), physical intuition necessitates that the positive root is taken so that disturbances decay away from topography rather than increase exponentially (see (2.17)).

If viscosity and diffusivity are non-zero, there is still a well-defined complex vertical wavenumber $m(k) = \pm \sqrt {Q(k)}$. However, since $m$ is now complex, the correct choice for the sign of $m$ must always be that with positive imaginary part so that the solution decays away from the topography. Note that the previously non-radiating wavenumbers gain a small radiating component, although this is insignificant in the realistic limit of weak viscosity and diffusivity – see Shakespeare & Hogg (Reference Shakespeare and Hogg2017) for a detailed discussion. Here, we consider radiating lee waves from topography such that $|\,f| < |U(0)k| < |N(0)|$.

When the coefficients $P$ and $Q$ are not constant in $z$, this radiating upper boundary condition is in general poorly defined, since for each $k$ there is not a well-defined vertical wavelength and group velocity. Waves can internally reflect and refract from changes in background density or velocity, so the solution cannot be restricted to upward propagating components. However, in some cases WKB solutions can be found for slowly varying background flows – see § 2.11.

If lee waves reach the upper ocean, the radiating upper boundary condition is inappropriate, and the air–sea interface may instead be better represented by a free surface boundary condition. A simpler condition is the rigid lid – we will show that for this problem, these are essentially equivalent.

At a free surface given by $z = H + \eta (x)$, where $\eta \ll H$, the linearised kinematic boundary condition (cf. (2.21)) is

(2.26)\begin{equation} \psi(x,H) = U(H)\eta(x). \end{equation}

A further boundary condition is required to close the problem and is given by the dynamic condition that the pressure at the surface is equal to the atmospheric pressure $p_A$ (assumed constant). The full pressure at the surface is $p^{{{\dagger}} }(x,z) = \bar {p}(z) + p(x,z)$, plus the linear term $-\rho _0 g (z-H)$ that was removed in the definition of $p^{{{\dagger}} }$ (see (2.1)). Expanding $p^{{{\dagger}} }(x,H + \eta (x))$ to first order in the perturbation variables and $\eta$, and using $\bar {p}(H)= p_A$, gives

(2.27)\begin{align} p^{{{{\dagger}}}}(x,H + \eta(x)) &\simeq p^{{{{\dagger}}}}(x,H) + \eta(x)\left.\frac{\partial {p^{{{{\dagger}}}}}}{\partial {z}}\right|_{z=H}, \end{align}

(2.28)\begin{align} &\simeq p_A + p(x,H) + \eta(x)\left.\frac{\partial {\bar{p}}}{\partial {z}}\right|_{z=H}. \end{align}

Invoking the boundary condition then gives

(2.29)\begin{equation} p_A = p_A + p(x,H) + \eta(x)\left.\frac{\partial {\bar{p}}}{\partial {z}}\right|_{z=H} - \rho_0 g \eta(x). \end{equation}

Using hydrostatic balance of the base state (2.5) then gives the dynamic boundary condition

(2.30)\begin{equation} p(x,H) = \rho_0 (g -\bar{b}(H))\eta(x) = \bar{\rho}(H)g\eta (x) = \rho_0 g \eta(x) , \end{equation}

where the reference density $\rho _0$ is taken to be the base state surface density.

Eliminating the unknown $\eta$ from the surface boundary conditions (2.26) and (2.30) gives the boundary condition

(2.31)\begin{equation} \psi(x,H) = \frac{U(H)p(x,H)}{\rho_0 g}. \end{equation}

This surface boundary condition could be used with the bottom boundary condition (2.21) to solve (2.7)–(2.11) , then the surface height recovered from (2.30) or (2.26). However, in practice this is unnecessary, as (2.31) can be well approximated by the rigid lid condition $\psi (x,H) = 0$, equivalent to imposing $\eta (x) = 0$ (and thereby not satisfying the dynamic boundary condition). To see why, first notice from (2.7) that for negligible rotation, shear and viscosity, $p \sim -\rho _0 U u \sim \rho _0 U \psi _z$. For slowly varying background conditions, we expect $\psi (x,z)$ to locally have a sinusoidal structure in the vertical, so let (for fixed $x$)

(2.32)\begin{gather} \psi(z) \sim A\sin (mz + \theta), \end{gather}

(2.33)\begin{gather}\psi_z(z) \sim Am\cos (mz + \theta), \end{gather}

for some amplitude $A$, wavenumber $m$ and phase $\theta$. Then, using the boundary relation (2.31):

(2.34)\begin{equation} \sin( mH + \theta) \sim \frac{mU^{2}}{g}\cos (mH + \theta). \end{equation}

Assuming that $m \sim N/U$ (the hydrostatic, non-rotating limit of (2.18))

(2.35)\begin{equation} \tan (mH + \theta) \sim \frac{NU}{g} \ll 1, \end{equation}

even for large upper ocean values of $U$ and $N$, and this scaling still holds for realistic conditions with rotation and non-hydrostatic waves. Therefore, the phase $\theta$ is such that $\psi (H) \simeq 0$, and it is clear that a rigid lid approximation is sufficient for determining the interior structure of the lee waves. The full free surface boundary condition could be implemented to determine exactly the (linear) surface height $\eta (x)$, but hereafter we only consider the rigid lid boundary condition. Since the interior flow is relatively unaffected by this approximation, we could still estimate the surface height without explicitly solving for it, using (2.30)

(2.36)\begin{equation} \eta(x) \sim \frac{p(x,H)}{\rho_0 g}, \end{equation}

where $p(x,H)$ is found from the rigid lid solution.

With the rigid lid condition $\psi (x,H) = 0$, the solution to the bounded problem is then given by (2.24)–(2.25), with the upper boundary condition

(2.37)\begin{equation} \hat{\zeta}(k,H) = 0. \end{equation}

2.5. Group velocities

The behaviour of lee waves in a bounded domain depends strongly on the direction of their group velocity. Consider the inviscid and unbounded problem with uniform background stratification and velocity, so that the vertical wavenumber $m$ is independent of $z$. Re-deriving the governing equation (2.16) with time dependence by considering plane wave solutions ${\sim } \exp ({\textrm {i}(kx + mz -\omega t)})$ gives the dispersion relation (cf. (2.18))

(2.38)\begin{equation} (\omega - Uk)^{2} = \frac{N^{2}k^{2} + f^{2}m^{2}}{\alpha k^{2} + m^{2}}, \end{equation}

where $\omega = 0$ for steady lee waves satisfying the boundary condition (2.21). The phase velocity is zero as a result, but the group velocity is non-zero and can be found by differentiating (2.38)

(2.39)\begin{align} \boldsymbol{c_g} &= \left(\frac{\partial {\omega}}{\partial {k}}, \frac{\partial {\omega}}{\partial {m}}\right), \end{align}

(2.40)\begin{align} &= \left( \frac{f^{2}(N^{2} - \alpha U^{2}k^{2}) + \alpha U^{2}k^{2}(U^{2}k^{2} - f^{2})}{Uk^{2}(N^{2} - \alpha f^{2})} , \frac{(U^{2}k^{2} - f^{2})^{3/2}(N^{2} - \alpha U^{2}k^{2})^{1/2}}{Uk^{2}(N^{2}-\alpha f^{2})}\right), \end{align}

where the sign of the vertical group velocity is taken to be positive when $m$ is real, as is appropriate for the unbounded case, and it is assumed that $|\,f| < |Uk| < |N|$ so that the waves are radiating. It is clear from (2.40) that, in the non-rotating and hydrostatic case ($f = \alpha = 0$), the horizontal component of group velocity is zero, and waves propagate vertically upwards. Supposing now that they encounter the surface, the waves will reflect and propagate directly downwards – still with zero horizontal group velocity and now with negative vertical group velocity – superimposing exactly on the upward propagating wave field. This scenario is illustrated for monochromatic topography in figure 2(a) (blue lines). The reflected waves can then be expected to directly increase or decrease the topographic wave drag and energy conversion by constructive or destructive interference with the upwards propagating wave field at the topography. The extent to which this occurs is determined by the energy lost to mixing and dissipation during propagation, to be discussed in § 4.1. When the horizontal group velocity is zero, no energy propagates downstream, so without dissipative energy loss there can be no energy conversion into lee waves at the topography and also no wave drag (see § 2.8). However, there may be resonance (to be discussed in § 2.7).

Figure 2. (a) Diagram showing monochromatic topography with indicative ray paths for several values of the overlap parameter $\gamma$, demonstrating some possible idealised paths of lee waves with different directions of group velocity. (b) Diagram showing the vertical structure function $\hat {\zeta }(k,z)$ for the analytic solution (2.42), for some vertical wavenumbers $m$ such that the solution is near resonant (blue) and at its minimum amplitude (magenta).

If $|Uk|$ is of comparable magnitude to the Coriolis or buoyancy frequency, the waves will have a positive horizontal component of group velocity and will propagate both upwards and downstream, reflecting at the surface downstream of the topography. Without dissipation and mixing this could continue indefinitely and allow the lee wave energy to propagate far downstream, although in reality it seems unlikely that a significant amount of wave energy would undergo multiple reflections due to nonlinear interactions near the bottom boundary. For flow moving over the top $\sim U/N$ of an isolated topographic peak and generating a continuous range of wavenumbers $k$, if the angle of propagation (the angle of the group velocity vector to the vertical) is large enough for all radiating components, the reflected wave will not significantly interact with the generation process and the wave drag will be unchanged from the unbounded case. If the bump is not isolated, the reflected wave could be incident on the generation of a lee wave at a different topographic feature, and the drag (and energy flux) modification would be more complex.

To determine the likelihood of a lee wave superimposing on itself at the topography, we can determine the angle of propagation using (2.40), assuming for simplicity that $U$ and $N$ are constant and viscosity and diffusivity are negligible. An ‘overlap parameter’ $\gamma (k)$ can be defined as

(2.41a,b)\begin{equation} \gamma (k) = \left|\frac{kH}{\rm \pi}\right|\tan \theta ,\quad \tan \theta = \left.\frac{\partial {\omega}}{\partial {k}} \right/ \frac{\partial {\omega}}{\partial {m}}, \end{equation}

where $\theta$ is the angle of wave propagation to the vertical, and for each $k$, $\gamma$ is the horizontal distance travelled by a wave whilst propagating to the surface at $z=H$ and back to the topography at $z=0$, normalised by the horizontal wavelength. This is illustrated for $\gamma = 0, 0.5$, and $1$ in figure 2(a). Figure 3 shows the variation in $\gamma$ with horizontal wavelength $2{\rm \pi} /k$ for various $U$ and $f$. Each curve tends to infinity (not shown) at $k = |N/U|$ and $k = |\,f/U|$, at which point the vertical group velocity reaches zero and the solutions become evanescent. The increase in $\gamma$ for both large and small horizontal wavelengths is due to the increasing downstream component of group velocity.

Figure 3. Overlap parameter $\gamma$ defined in (2.41a,b) for the non-hydrostatic case and (a) fixed $U = 0.1\ \textrm {m}\,\textrm {s}^{-1}$, varying $f$ and (b) fixed $f = -1\times 10^{-4}\ \textrm {s}^{-1}$, varying $U$, with $H = 3\ \textrm {km}$, $N = 1\times 10^{-3}\ \textrm {s}^{-1}$. The black dashed line shows $\gamma = 1$, below which a reflected lee wave could be expected to modify its own generation mechanism.

For smaller values of $f$ (black dashed and orange lines in figure 3a) and larger values of $U$ (blue and magenta lines in figure 3b), there exists a range of scales at which $\gamma \lesssim 1$, indicating that reflected lee waves could impact on the generation mechanism by direct superposition. For $f = -1\times 10^{-4}\ \textrm {s}^{-1}$, characteristic of the SO, there exist horizontal scales at which this may be the case for $U \gtrsim 0.2\ \textrm {m}\,\textrm {s}^{-1}$. However, for $f = -1\times 10^{-4}\ \textrm {s}^{-1}$ and $U = 0.1\ \textrm {m}\,\textrm {s}^{-1}$ (orange line in figure 3b), $\gamma > 2$ for all radiating wavelengths, and all reflected waves return to topography at least $4\ \textrm {km}$ downstream of the generating topographic feature. Of course, this argument does not cover the more likely scenario of varying velocity and stratification with height. We conclude that for lee waves in shallow areas, low latitudes or high background flows it is possible for lee waves generated by isolated topography to reflect at the surface and modify the original wave drag and energy conversion, but that this is unlikely for deep generation, low background velocities and high latitudes.

When the topography is not isolated, and in particular when an artificially discrete topographic spectrum is used as in this study, the wave drag modification can be significant even when the overlap parameter is larger than one. For monochromatic topography, when $\gamma (k) = n \in \mathbb {N}$, a wave generated at a topographic peak reflects at the surface and is incident on the topographic peak $n$ wavelengths downstream from the original, as shown in figure 2(a) for $n=1$, and impacts the wave field at that generation site in a similar way to the case $\gamma = 0$.

2.6. Analytic solution

When $U$ and $N$ are constant with height, so that $P$ vanishes and $Q$ is a function of $k$ only, the solution to (2.24) for $\hat {\zeta }(k,z)$ is (extended from Baines Reference Baines1995)

(2.42)\begin{equation} \hat{\zeta}(k,z) = \frac{\sin(m(k)(H-z))}{\sin ( m(k)H )}, \end{equation}

where $m(k)$ is the complex vertical wavenumber defined by $m^{2}(k) = Q(k)$, and the choice of sign does not matter. The solution can then be found numerically for general topography via (2.23), or analytically for monochromatic topography $h(x) = h_0\cos k_0 x$ to be

(2.43)\begin{equation} \psi(x,z) = Uh_0 \mathrm{Re} \left(\frac{\sin(m(k_0)(H-z))}{\sin ( m(k_0)H )}\textrm{e}^{\textrm{i}k_0 x}\right). \end{equation}

The above solutions are valid only when $|m(k)H| \neq n{\rm \pi}$, $n \in \mathbb {N}$. At such points, resonances of the system occur.

2.7. Resonance

For uniform $U$ and $N$, under the assumption that lee waves are hydrostatic ($|Uk| \ll |N|$), rotation is unimportant ($|Uk| \gg |\,f|$), and the system is inviscid, the vertical wavenumber is simply $m(k) = N/U$. The resonances of (2.42) are then independent of $k$, and occur when $|NH/U| = n{\rm \pi}$ for some $n \in \mathbb {N}$. There are no steady solutions to (2.24), (2.25) and (2.37) in this inviscid limit. Physically, this occurs when a whole number of half-wavelengths fits in the vertical domain and there is constructive interference of the upwards and downwards propagating waves. Figure 2(b) shows the vertical structure function $\hat {\zeta }$, defined in (2.42), for two real values of $m$. When $mH = 5.05{\rm \pi}$ (blue) the system is near resonance, as the half-wavelength nearly divides the depth $H$ (true resonance is at $mH = 5{\rm \pi}$). Thus, $\hat {\zeta }(z) = 0$ near $z=0$, so in order to satisfy the boundary condition $\hat {\zeta }(z=0) = 1$, the amplitude of the wave must be very large. At true resonance, this boundary condition cannot be met. In the opposite case, (shown for $mH = 5.5{\rm \pi}$ in magenta), there is destructive interference and the amplitude is at a minimum.

Under the above assumptions, the horizontal group velocity is zero, therefore energy cannot escape downstream and the wave generation at resonance continually reinforces the wave field. If this were to happen in practice, the wave amplitude would become large enough to invalidate the linearity of the wave field, perhaps causing nonlinear wave breaking or modifying the wave field or the boundary condition so as to move the system away from resonance.

When non-hydrostaticity is included, the horizontal group velocity is non-zero and the nature of the resonance changes slightly. The vertical wavenumber $m(k) = \sqrt {N^{2}/U^{2} - k^{2}}$, thus the solution (2.42) has singularities at

(2.44)\begin{equation} k^{2} = \frac{N^{2}}{U^{2}} - \frac{n^{2}{\rm \pi}^{2}}{H^{2}}, \quad n \in \mathbb{N}\,. \end{equation}

Physically, these singularities still represent modes where an exact number of half vertical wavelengths fit in the domain, but now this happens at different values of $U$, $N$ and $H$ for each component $k$ of the wave field.

Mathematically, the resulting singularities of (2.42) are simple poles, so when the topographic spectrum $\hat {h}(k)$ is continuous (as for isolated topography), the integral (2.23) along the real line can be moved to a contour of integration in complex $k$ space that avoids the poles. To ensure that there is no disturbance at upstream infinity, the contour must be taken below rather than above the poles (McIntyre Reference McIntyre1972). The solution can then be expressed as the Cauchy principal value of (2.23) plus half the residues of the simple poles, which represent the non-hydrostatic resonant modes (Baines Reference Baines1995). The solutions, when $\hat {h}(k)$ is continuous, could be found numerically from (2.23) by choosing some contour of integration sufficiently far from the poles to avoid numerical difficulties. However, this becomes more difficult once rotation is included since the poles no longer all lie on the real axis. The numerical solution is also problematic since periodicity in the horizontal is assumed by default when taking a discretised Fourier transform, leading to spurious waves upstream of the isolated topography. If the topographic spectrum is discrete and includes one of the singular wavelengths defined by (2.44), then true resonance occurs and no steady solution exists.

The inclusion of energy loss through viscosity and diffusivity aids the numerical solution by moving all poles off of the real line so that the integral (2.23) can be found numerically with a simple fast Fourier transform (FFT). Although true resonance is avoided, states can still be near resonant, as will be shown in § 4.1. The topographic representation used here (see § 3.3) consists of a spectrum of topographic wavenumbers, which numerically becomes a sum of discrete components. This is likely to enhance the resonance effect compared with a more realistic and inhomogeneous topography.

2.8. Energy and momentum

The vertical linear lee wave energy flux at a given height is given by $\overline {pw}$, where an overbar here represents a horizontal average. At the topography ($z=0$), this is equal to the product of the bottom mean flow velocity and the horizontally averaged form drag exerted by the topography on the mean flow, since using (2.21)

(2.45)\begin{equation} \overline{pw}|_{z=0} = U(0)\overline{ph_x}|_{z=0}. \end{equation}

Taking the inner product of (2.7)–(2.9) with the perturbation velocity and multiplying (2.10) by the perturbation buoyancy gives the energy equation for the wave field. Taking a horizontal average and assuming a periodic domain in the horizontal then gives an expression for the divergence of the energy flux

(2.46)\begin{equation} \overline{pw}_z ={-}\rho_0(U_zF + \bar{D}), \end{equation}

where $\bar {D} = \overline {\varepsilon } + \bar {\varPhi }$ is the horizontally averaged energy loss from the flow, consisting of the dissipation rate $\varepsilon = \mathcal {A}_h|\boldsymbol {u}_x|^{2}$ and irreversible mixing $\varPhi = \mathcal {D}_hb_x^{2}/N^{2}$, and

(2.47)\begin{equation} F = \overline{uw} - \frac{f\overline{vb}}{N^{2}} \end{equation}

is the wave pseudomomentum flux, or the Eliassen–Palm (E–P) flux (Eliassen & Palm Reference Eliassen and Palm1960). If there are no critical levels ($U \neq 0$) it can be shown from (2.7)–(2.10) that the E–P flux $F$ is related to the energy flux as (extended from Eliassen & Palm Reference Eliassen and Palm1960)

(2.48)\begin{align} \overline{pw} &={-}\rho_0U\left[ \overline{\left(u - u_x\mathcal{A}_h/U\right)w} - \frac{f}{N^{2}}\overline{\left(b - \mathcal{D}_hb_x/U\right)v}\right] \end{align}

(2.49)\begin{align} &={-}\rho_0UF\left(1 + {O}\left(\mathcal{A}_hk/U\right)\right), \end{align}

where $k$ is the characteristic wavenumber of the topography, and $\mathcal {A}_h \sim \mathcal {D}_h$. Taking typical values considered here, $\mathcal {A}_h \sim 1\ \textrm {m}^{2}\,\textrm {s}^{-1}$, $k \sim 0.005\ \textrm {m}^{-1}$, and $U \sim 0.1\ \textrm {m}\,\textrm {s}^{-1}$, giving $\mathcal {A}_h k/U \sim 0.05 \ll 1$. Thus the energy flux is approximately equal to the local velocity multiplied by the E–P flux even when there is energy loss.

Since the source of the waves is at the bottom of the domain, the energy flux is always positive (or zero). (2.49) then gives that the E–P flux must have the opposite sign to $U$, thus $F \leq 0$ in the cases considered here. This can be seen from the form of $F$ (2.47) in the non-rotating case, since a positive vertical velocity perturbation of the flow ($w > 0$) will correspond to the deceleration of a fluid parcel in the horizontal ($u < 0$), thus the momentum flux $\overline {uw} \leq 0$.

In the inviscid problem, (2.46) and (2.49) together give (Eliassen & Palm Reference Eliassen and Palm1960)

(2.50)\begin{equation} F_z =0. \end{equation}

Therefore, the E–P flux is conserved when there is no energy lost to dissipation and mixing. When there is also no vertical shear of the mean flow ($U_z = 0$), (2.46) gives that the energy flux $\overline {pw}$ is also conserved. When the mean velocity increases or decreases with height, the energy flux increases or decreases correspondingly, but the E–P flux is still conserved. Any divergence of the E–P flux thus corresponds to the force exerted on the flow by the waves as they dissipate (Andrews & McIntyre Reference Andrews and McIntyre1976). It is the divergence of $F$ rather than the Reynolds stress (or momentum flux) $\overline {uw}$ that gives the relevant lee wave forcing on the mean flow, since $\overline {uw}$ is in general not conserved – a paradox explained by Bretherton (Reference Bretherton1969).

The total wave drag on the mean flow (defined here to be a positive quantity, although acting in the negative $x$ direction) is therefore given by the integral of $\rho _0 F_z$ over the depth of the ocean. Since there cannot be any energy or momentum flux through the upper boundary, $\overline {pw}|_{z = H} = F(H) = 0$, thus the wave drag is given by $-\rho _0 F(0)$. Comparison of (2.45) and (2.49) then shows that up to ${O}(\mathcal {A}_h k/U)$ the wave drag is equal to the form drag.

Since $\overline {pw}|_{z = H} = F(H) = 0$, if energy loss $\bar {D}$ is zero, $F = 0$ everywhere (from (2.50)) and $\overline {pw} = 0$ everywhere (from (2.49)), thus there is no net topographic wave drag on the flow or energy conversion to lee waves in steady state. This is true for a periodic domain – in a non-periodic domain the waves would propagate infinitely far downstream, and boundary fluxes would become important in the (2.46) and (2.49). Energy loss is therefore a key component of the bounded study, as there can be no topographic wave drag without it. This is a realisation of the ‘non-acceleration theorem’, first discussed by Charney & Drazin (Reference Charney and Drazin1961). Of course, in the unbounded problem there must also be energy loss in order to have wave drag at the topography – but that energy loss can implicitly occur by allowing the lee waves to exit the given domain (such that $\overline {pw}|_{z =H} > 0$) and dissipate ‘elsewhere’.

From (2.46), the wave energy flux can change both by exchange with a mean flow through the E–P flux (Kunze & Lien Reference Kunze and Lien2019), and by mixing and dissipation. Integrating (2.46) over the entire height of the domain gives

(2.51)\begin{equation} \overline{pw}|_{z=0} -\overline{pw}|_{z=H} = \rho_0\int_0^{H} U_zF + D\, \textrm{d}z. \end{equation}

If there is an upper boundary and no background shear ($U_z = 0$), then $\overline {pw}|_{z=H} = 0$, and topographic energy conversion and wave drag are directly proportional to the total mixing and dissipation in the water column.

2.9. Time dependence

When calculating lee wave fluxes, it is usually assumed that the background fields and lee waves are steady. In reality, the geostrophic flows that generate lee waves vary on time scales of days to weeks. After a change in the background flow, the time taken for the lee wave field to equilibrate to the new steady state could be long compared with the typical time scales of the flow.

The relevant time scale in this study is the time taken for the lee wave to propagate from the topography to the surface. Figure 4 shows the vertical group velocity (defined in (2.40)) for various values of $f$ and $U$. Lee waves generated at smaller horizontal scales (larger $k$) propagate faster, although they are also more likely to dissipate along the way due to sharper horizontal gradients. The effect of rotation on larger scales significantly slows the vertical group velocity, so that larger horizontal scale waves will take significantly longer to develop. The group velocity increases with $U$, so larger background velocities allow faster lee wave propagation. For a wave with wavelength $3\ \textrm {km}$ in a background flow $U = 0.1\ \textrm {m}\,\textrm {s}^{-1}$, $N = 1 \times 10^{-3}\ \textrm {s}^{-1}$ and $f = -1 \times 10^{-4}\ \textrm {s}^{-1}$, the vertical group velocity is approximately $1\ \textrm {km}$ per day, suggesting that a wave would take 3 days to propagate to the surface and a further 3 days to reflect back to the topography in an ocean of depth $3000\ \textrm {m}$. For a vertically varying background flow as is typical of the ocean, the vertical group velocity and wavelength will change during propagation, modifying this result. The time scale separation of full water-column lee wave formation and the mesoscale eddy field is therefore not clear, and depends on the scale of the waves. However, in energetic regions of the ocean such as the Drake Passage shown in figure 1, large velocities can enable high vertical group velocities and the steady approximation for lee waves at certain scales is expected to be valid.

Figure 4. Vertical group velocity defined in (2.40) for the non-hydrostatic case and (a) fixed $U = 0.1\ \textrm {m}\,\textrm {s}^{-1}$, varying $f$ and (b) fixed $f = -1 \times 10^{-4}\ \textrm {s}^{-1}$, varying $U$, with $H = 3\ \textrm {km}$, $N = 1 \times 10^{-3}\ \textrm {s}^{-1}$.

2.10. Critical levels

In the inviscid and non-rotating problem, there are singularities of (2.13)–(2.15) at levels where $U = 0$ (Booker & Bretherton Reference Booker and Bretherton1967; Maslowe Reference Maslowe1986). These are known as critical levels, where the horizontal phase speed of the wave (here equal to zero) equals the mean flow speed. At these levels, the vertical wavelength and group velocity vanish. No energy or momentum flux at the original wavenumber can propagate any further vertically, and the perturbation velocities become very large, invalidating the linear solution. In reality, instabilities and energy loss can lead to wave breaking and reflection at this level (Wurtele Reference Wurtele1996), thus critical levels may be a sink of lee wave energy in the ocean (Bell Reference Bell1975). However, this requires that the mean flow speed reaches zero somewhere in the water column. This is not a ubiquitous feature of the geostrophic eddies of the ACC, although critical levels may exist. This mechanism may be more important in regions of layered currents such as near the equator or in western boundary currents.

When rotation is included, there exist two further singularities of (2.13) at $U= \pm |\,f/k|$, above and below the critical level $U = 0$ (Jones Reference Jones1967). These act to prevent the vertical propagation of the wavenumber $k$, in a similar way to the critical level of the non-rotating problem at $U = 0$. However, since each critical level is specific to the wavenumber $k$ (unlike for the non-rotating problem), if the spectrum of the topography is continuous it can be shown that there need not be singularities of the linear problem at these critical levels since the relevant solutions of (2.13) are logarithmic and thus integrable over a spectrum (Wurtele, Datta & Sharman Reference Wurtele, Datta and Sharman1996). Therefore, in reality there is not a single well-defined critical level for lee waves with rotation and a continuous spectrum of wavenumbers. However, the solutions may still become nonlinear so as to invalidate the linear solution and cause breaking. It can also be shown that when $f \neq 0$, the solution at $U = 0$ is no longer singular (Grimshaw Reference Grimshaw1975) – physically this is because all components have already reached their first critical level and stopped propagating.

When the flow is sheared such that $|U|$ decreases with height, energy transfers from the lee waves to the mean flow via the E–P flux (see (2.46)), leaving less energy to be dissipated at the critical level for a particular wavenumber $k$. Kunze & Lien (Reference Kunze and Lien2019) examine this mechanism as a possible sink for lee wave energy in regions of intensified bottom flow. In particular, for lee wave energy generated at wavenumbers far from the inertial limit $|Uk| = |\,f|$, a greater proportion of the initial energy is available to be reabsorbed by a mean flow decreasing with height, allowing a smaller percentage to be dissipated at the critical level at $|Uk| = |\,f|$ or elsewhere. For waves generated close to the inertial limit (from large-scale topography), little energy is available for transfer to the mean flow, as it will instead soon reach its critical level and dissipate.

The inclusion of viscosity and diffusivity allows non-singular solutions to be found at critical levels where $|Uk| = |\,f|$. However, near these levels the wave fields can become nonlinear, invalidating the linear approach. Furthermore, on the approach to these levels the vertical wavelength tends to zero, creating sharp vertical gradients and enhancing energy loss. Having neglected vertical viscosity and diffusivity in our solution, this energy loss does not take place. When the horizontal viscosity and diffusivity are large enough and shear small enough for the solutions to stay appropriately linear as a critical level is approached, the linear solution is valid, but may be unrealistic due to the lack of vertical dissipation and mixing. Velocity profiles that decrease with height are therefore not considered hereafter.

For positively sheared background flows where the flow speed increases with height, energy instead transfers from the mean flow to the lee waves during propagation (see (2.46)). Since wind driven oceanic currents tend to be surface intensified, this may be a common occurrence. In this case (or if stratification $N$ decreases with height), the waves may reach ‘turning levels’, whereby their intrinsic frequency $|Uk|$ reaches the buoyancy frequency $N$ (Scorer Reference Scorer1949). At such levels the vertical wavenumber $m$ tends to zero (see (2.18)), and the wave is reflected downward. Scorer (Reference Scorer1949) showed that wave amplitudes in the resulting ‘trapped’ wave field could be increased by the superposition of reflected waves, much like in the current study due to the upper boundary. These turning levels are not the focus of our study, but may occur in the solutions for certain wavenumbers.

2.11. WKB solutions

Before solving the full equation (2.24) numerically, we first consider how much progress can be made with WKB theory (Gill Reference Gill1982). Under the assumption that waves are propagating through a slowly varying medium such that the wavelength is small compared with the scale of changes in the background flow, WKB theory can often be used to find closed form solutions. In addition to the necessary assumption for linear theory that there are no critical levels (where $|U(z)k| = |\,f|$), WKB theory requires the further assumption that there are no turning levels (where $|U(z)k| = |N|$), since the vertical wavelength tends to infinity at these levels (see (2.18)). However, under the assumption that waves are generated and stay within the propagating range (2.19), solutions with both the rigid lid and freely radiating boundary condition can be found. For the latter, the component of the solution with positive local vertical group velocity is taken as in the uniform background case. de Marez et al. (Reference de Marez, Lahaye and Gula2020) find inviscid and non-rotating lee wave solutions under this approximation with a freely radiating boundary condition. Here, we extend their analysis to include rotation and viscous terms.

First, we write (for each $k$)

(2.52)\begin{equation} \hat{\zeta}(z) = A(z)\textrm{e}^{\textrm{i}\varphi(z)}, \end{equation}

where $A(z)$ and $\varphi (z)$ are some real amplitude and phase such that $\varphi '(z) \sim m \sim 1/\lambda _z$, where $\lambda _z$ is the vertical wavelength of the waves, itself varying over some length scale $L \gg \lambda _z$; $A(z)$ also varies on the length scale $L$ due to changes in the background flow and to viscosity. At this point we must therefore also take the weakly viscous approximation $\mathcal {A}_h^{2} k^{2}/U^{2} \ll 1$ to ensure that the decay scale due to viscosity is larger than the vertical wavelength of the waves (see (2.56) and discussion for justification). This is likely to be realistic in the ocean, as discussed after (2.49). Substituting (2.52) into (2.24), taking real and imaginary parts, and neglecting second-order terms in $\lambda _z /L$ allows two independent solutions to be found

(2.53)\begin{equation} \hat{\zeta}_\pm(k,z) = Q_r^{-{1}/{4}} \exp\left({\int_0^{z} -\tfrac{1}{2}(P_r \pm Q_i/\sqrt{Q_r}) \pm \textrm{i} \sqrt{Q_r} \, \textrm{d}\hat{z}}\right), \end{equation}

where $\hat {z}$ is a dummy integration variable, and $Q_r, P_r$ and $Q_i, P_i$ are the real and imaginary parts of $P$ and $Q$, defined in (2.14) and (2.15). Setting $\mathcal {A}_h = \mathcal {D}_h$, and using the weakly viscous approximation, these become

(2.54a,b)\begin{gather} Q_r = \frac{k^{2}(N^{2} - \alpha U^{2} k^{2} - U U_{zz})}{U^{2}k^{2} - f^{2}}, \quad Q_i = \frac{\mathcal{A}_h k^{3}\left[ 2Uk^{2}(N^{2} - \alpha f^{2}) - U_{zz}(U^{2}k^{2} + f^{2})\right]}{(U^{2}k^{2} - f^{2})^{2}}, \end{gather}

(2.55a,b)\begin{gather} P_r = \frac{2f^{2}U_z}{U(U^{2}k^{2} - f^{2})}, \quad P_i = \frac{2\mathcal{A}_h f^{2} U_z k (3U^{2} k^{2} - f^{2})}{U^{2}(U^{2} k^{2} - f^{2})^{2}} . \end{gather}

Each term in (2.53) can be interpreted to understand the solutions; $Q_r^{-{1}/{4}} \exp ({-\frac {1}{2}\int _0^{z} P_r\, \textrm {d}\hat {z}})$ determines the change in amplitude of solution due to changes in the background flow, and is constant in a uniform flow; $\exp ({\mp \frac {1}{2} \int _0^{z} Q_i/\sqrt {Q_r} \, d\hat {z}})$ determines how the wave amplitude exponentially decays (or ‘grows’, for a downwards propagating component) due to viscosity; $\exp ({\pm \textrm {i} \int _0^{z} \sqrt {Q_r} \, \textrm {d}\hat {z} })$ is the oscillatory component, with wavenumber $\sim \sqrt {Q_r}$ (cf. (2.17)).

The solution to the bounded problem can then be found using bottom and upper boundary conditions (2.25) and (2.37), and the solution to the unbounded problem can be found by taking $\zeta _+(k,z)$ (see discussion in § 2.4.2) and using bottom boundary condition (2.25).

Simplifying further by considering the non-rotating and hydrostatic case with $f = \alpha = 0$, and assuming a linear velocity profile with $U_{zz} = 0$, the WKB solutions are given by

(2.56)\begin{equation} \hat{\zeta}_\pm(k,z) = \sqrt{\frac{U}{N}} \exp\left({\pm \int_0^{z} \frac{N}{U}\left(\textrm{i} - \frac{\mathcal{A}_h k}{U}\right) \,\textrm{d}\hat{z}}\right). \end{equation}

Aside from reflections and viscous dissipation, we expect that as $U$ and $N$ vary, $w \sim \hat {\zeta } \sim \sqrt {U/N}$, $u \sim \hat {\zeta }_z \sim \sqrt {N/U}$ and $\overline {uw} \sim \text {constant}$, as expected by conservation of the E–P flux (2.50). The length scale of viscous decay is given by $U^{2}/\mathcal {A}_h N k$, and the ratio of the vertical wavelength to the viscous decay scale is $\mathcal {A}_h k /U$, justifying the condition $\mathcal {A}_h^{2} k^{2} / U^{2} \ll 1$ for use of the weakly viscous approximation.

Figure 5 shows the WKB non-hydrostatic, rotating, bounded solution for $\hat {\zeta }(k,z)$ and the corresponding full numerical solution (to be explained in § 3) for three different values of $k$ and various background profiles of $U(z)$ and $N(z)$. In figure 5(a,c,d) each of the wavenumbers $k_1, k_2, k_3$ remains within the propagating range (2.19) throughout the vertical domain, and the solutions are oscillatory. In each case, as $k$ increases, the amplitude of the solutions decreases, since the length scale of viscous decay decreases (cf. (2.56)). The exponential decay when $k=k_3$ (blue lines) is large enough that oscillatory component is small in comparison.

Figure 5. Numerical and WKB solutions $\mathrm {Re} [\hat{\zeta}(k,z)]$ to the rotating, non-hydrostatic, bounded problem defined in (2.24), with boundary conditions (2.25) and (2.37). Background fields vary linearly from $U(0) = 0.1\ \textrm {m}\,\textrm {s}^{-1}$ and $N(0) = 1 \times 10^{-3}\ \textrm {s}^{-1}$; (a) $U(H) = U(0)$ and $N(H) = N(0)$, (b) $U(H) = 0.3\ \textrm {m}\,\textrm {s}^{-1}$ and $N(H) = N(0)$, (c) $U(H) = U(0)$ and $N(H) = 3 \times 10^{-3}\ \textrm {s}^{-1}$ and (d) $U(H) = 0.3\ \textrm {m}\,\textrm {s}^{-1}$ and $N(H) = 3 \times 10^{-3}\ \textrm {s}^{-1}$. Solutions are shown at $k = \{k_1 = k_{min} + 0.1(k_{max}-k_{min}), k_2 = 0.5*(k_{min} + k_{max}), k_3 = k_{max} - 0.1(k_{max} - k_{min})\}$, where $k_{min} = |f/U(0)|$, $k_{max} = |N(0)/U(0)|$, and $f = -1 \times 10^{-4}\ \textrm {s}^{-1}$.

In general, the numerical and WKB solutions agree well, although to a lesser extent as $k$ increases since $m$ decreases (cf. (2.18)), and thus $\lambda _z$ increases, invalidating the WKB approximation. In figure 5(b) $U(z)$ increases and the solutions do not all stay within in the propagating range (2.19). Components $k_2$ and $k_3$ reach their turning levels where $|Uk| = |N|$, and the WKB approximation is not valid for these solutions.

Comparing the $k=k_1$ component in figure 5(a) (constant background fields) with figures 5(b) ($U$ increasing) and 5(c) ($N$ increasing), we see that the scalings inferred from (2.53) hold, and the amplitude of $\hat {\zeta }$ (and $w$) increases with increasing $U$ and decreases with increasing $N$. The vertical wavelengths scale with $U/N$, although they are also dependent on $k$, as is expected from (2.18).

A further observation from figure 5(a,d) is the reduced impact of viscosity and larger amplitudes when $U$ and $N$ increase together, even though $U/N$ is fixed. This is because the vertical group velocity of the waves is larger in this case – from (2.40) we can see that for $f^{2} \ll U^{2}k^{2} \ll N^{2}$, the vertical group velocity ${\sim } U^{2}k/N$. The waves have therefore decayed less by the time they reach and reflect from the surface. Equivalently, the decay scale due to viscous dissipation $U^{2}/\mathcal {A}_h N k$ increases with $U$ even when $U/N$ is fixed.

The WKB solutions are insightful provided there are no turning levels, and offer closed form solutions and insights that are not available from the numerical solutions. Similar approximations could also be made to implement a vertical and/or vertically variable viscosity, although we do not pursue this here to avoid further complexity. WKB solutions may therefore be more useful than numerical solutions in any eventual parametrisation of this process. However, hereafter we employ the full numerical solution, with no assumptions on the length scale of variations of the background flow or the viscous decay scale in order to study a wider range of wavenumbers and background flows.

3.1. Unbounded solver

The solution can be found in the traditional way (Bell Reference Bell1975), with the requirement the solution decays away from the topography as discussed in § 2.4. The solution for $\psi$ is given by (2.23), where $\hat {\zeta }(k,z)$ satisfies (2.24)–(2.25), with the radiating upper boundary condition satisfied by taking the correct choice of branch for $m$. Note that $P(k,z) = 0$ and $Q(k,z) = Q(k)$, so the solution for $\hat {\zeta }$ is simply

(3.1)\begin{equation} \hat{\zeta}(k,z) = \textrm{e}^{\textrm{i}m(k)z}, \end{equation}

where $m^{2}(k) = Q(k)$ and $\mathrm {Im}(m) > 0$. This can be implemented numerically for general topography 0$h(x)$ by performing the Fourier transforms with a FFT. Once $\psi$ is found, all other wave fields can be recovered.

3.2. Bounded solver

When the background flow is uniform in $z$, the solutions can be found similarly to the unbounded case above, using the analytic solution (2.42) for $\hat {\zeta }(k,z)$. When $U$ and $N$ are not uniform, we use Galerkin methods to solve (2.24)–(2.37), an unforced second-order ordinary differential equation with inhomogeneous boundary conditions. First, we transform it into a forced problem with homogeneous boundary conditions. Let

(3.2)\begin{equation} \hat{\zeta}(k,z) = \hat{\phi}(k,z) + G(k,z), \end{equation}

where $G$ is some function such that $G(k,0) = 1$ and $G(k,H) = 0$. Then $\hat {\phi }$ satisfies

(3.3)\begin{gather} \hat{\phi}_{zz} + P(k,z)\hat{\phi}_z + Q(k,z)\hat{\phi} = R(k,z), \end{gather}

(3.4)\begin{gather}\hat{\phi}(k,0) = 0, \end{gather}

(3.5)\begin{gather}\hat{\phi}(k,H) = 0 , \end{gather}

and $R$ satisfies

(3.6)\begin{equation} G_{zz} + P(k,z)G_z + Q(k,z)G ={-} R(k,z). \end{equation}

$G$ can be chosen to be any function satisfying $G(k,0) = 1$ and $G(k,H) = 0$. We choose it so that $R(k,0) = R(k,H) = 0$ by taking $G$ to be a cubic polynomial in $z$, and solving for the coefficients; $R$ can then be found via (3.6). The problem (3.3)–(3.5) can now be solved numerically using Galerkin methods. Specifically, for each $k$ we decompose $\hat {\phi }$, $P$, $Q$ and $R$ into finite Fourier sums with some truncation limit $M$

(3.7)\begin{equation} \left. \begin{aligned} \hat{\phi}(k,z) & = \sum_{s=1}^{M}a_s(k)\sin \frac{s{\rm \pi} z}{H}, \quad P(k,z) = \sum_{j=1}^{M}p_j(k)\sin \frac{(j-1){\rm \pi} z}{H}, \\ Q(k,z) & = \sum_{i=1}^{M}q_i(k)\cos \frac{(i-1){\rm \pi} z}{H}, \quad R(k,z) = \sum_{n=1}^{M}r_n(k)\sin \frac{(n-1){\rm \pi} z}{H} , \end{aligned} \right\} \end{equation}

where the $q_i$, $p_j$ and $r_n$ are known and found via the relevant sine or cosine transform, and the coefficients $a_s$ are to be found. Notice that the sine expansion of $\hat {\phi }$ and $R$ ensures that their boundary conditions are satisfied. However, if $P \neq 0$ or $Q_z \neq 0$ at $z = 0,H$, the sine and cosine expansions of $P$ and $Q$ respectively must represent one or more discontinuities in $P$ or $Q_z$ at endpoints. The numerical solution is therefore an approximation that is valid only in the interior, although (3.3) is satisfied everywhere by the series expansions. There can also be noise at the frequency of the truncation limit near the endpoints of the series representations due to the Gibbs phenomenon. With increasing truncation limit and vertical resolution, the interior series solution approaches the actual solution at all interior points – (2.46) can be used to validate this. As a consequence, quantities should be evaluated with care at $z=0,H$, and (2.51) is used to find the wave drag rather than direct evaluation at $z=0$.

Substituting (3.7) into (3.3), integrating over $z \in [0,H]$ and using the orthogonality properties of sine and cosine gives a matrix equation for the coefficients $a_s$

(3.8)\begin{equation} A_{sn}a_n = B_{sn}r_n, \end{equation}

where

(3.9)\begin{gather} A_{sn} ={-}\left(\frac{n{\rm \pi}}{H}\right)^{2}\delta_{s,n} \!+\! \frac{n{\rm \pi}}{2H}(p_{s-n+1} \!+\! p_{s+n+1} - p_{n-s+1}) + \frac{1}{2}(q_{s-n+1} + q_{n-s+1} - q_{s+n+1}), \end{gather}

(3.10)\begin{gather}B_{sn} = \delta_{n,s+1}. \end{gather}

The coefficients $a_s$ can now be found from (3.8) by inverting the matrix $A$; $\hat {\phi }$ can then be recovered from the coefficients $a_s$, $\hat {\zeta }$ found from (3.2) and $\hat {\psi }$ found from (2.23).

3.3. Topography

The topography is found from the theoretical abyssal hill topographic spectrum proposed by Goff & Jordan (Reference Goff and Jordan1988), and is similar to that used in previous lee wave modelling studies (Nikurashin & Ferrari Reference Nikurashin and Ferrari2010a, Reference Nikurashin and Ferrari2011; Nikurashin et al. Reference Nikurashin, Ferrari, Grisouard and Polzin2014; Klymak Reference Klymak2018; Zheng & Nikurashin Reference Zheng and Nikurashin2019). This choice allows both comparison to the results of these studies, and consideration of the behaviour of a range of different wavenumbers without expanding the parameter space of investigation. The 1-D topography is found from the theoretical topographic spectrum $P_{2D}(k,l)$

(3.11)\begin{equation} P_{2D}(k,l) = \frac{2{\rm \pi} h_0^{2} (\mu - 2)}{k_0 l_0}\left(1 + \frac{k^{2}}{k_0^{2}} + \frac{l^{2}}{l_0^{2}}\right)^{-({\mu}/{2})}, \end{equation}

by integrating over wavenumbers $l$; $k_0$ and $l_0$ are the characteristic horizontal wavenumbers, $\mu$ is the high wavenumber spectral slope and $h_0$ is the r.m.s. abyssal hill height. For comparison with other recent lee wave studies, we set $k_0 = 2.3\times 10^{-4}\ \textrm {m}^{-1}$, $l_0 = 1.3\times 10^{-4}\ \textrm {m}^{-1}$ and $\mu = 3.5$, in line with representative parameters of the Drake Passage region used in Nikurashin & Ferrari (Reference Nikurashin and Ferrari2010b) and Zheng & Nikurashin (Reference Zheng and Nikurashin2019). Next, $P_{1D}(k)$ is set to zero for wavenumbers $k$ such that $|U(0)k| < |\,f|$ or $|U(0)k| > |N(0)|$, since solutions in these ranges are non-propagating. The typical values used are $N(0) =1\times 10^{-3}\ \textrm {s}^{-1}$, $U(0) = 0.1\ \textrm {m}\,\textrm {s}^{-1}$ and $f = -1\times 10^{-4}\ \textrm {s}^{-1}$, corresponding to a topography with wavelengths between $\sim 630\ \textrm {m}$ and $\sim 6300\ \textrm {m}$. Note that the same topography is used throughout, even when $f =0$.

The topographic height differs from that used in the aforementioned studies, since the solver is linear and the solutions must therefore remain approximately linear to be valid. We normalise the topography resulting from the above steps so that the r.m.s. of the final topography $h_{rms} = 25\ \textrm {m}$. This gives a Froude number $Fr_L = Nh_{rms}/U = 0.25$ and is sufficient to keep the solution near linear such that the perturbation horizontal velocity $u$ is less than the background velocity $U$, with the exception of resonant cases. This is an unrealistically low Froude number for the rough topography of many parts of the SO (Nikurashin & Ferrari Reference Nikurashin and Ferrari2010b), but since the perturbation quantities are linear in $\hat {h}(k)$ (e.g. (2.23)), simple scaling arguments can recover the dependence on $h_{rms}$. The goal of this study is not to make predictions of the actual magnitude of the lee wave field, but its structure in the vertical and dependence on viscosity and diffusivity, background fields and boundary condition.

3.4. Numerical set-up

In the following section, the numerical solver is used to solve for the wave fields in a 2-D domain of width $40\ \textrm {km}$, and height $3\ \textrm {km}$. The number of grid points in $x$ and $k$ is 800, and in $z$ is 257. The truncation limit $M$ (see (3.7)) is 200. Sensitivity tests were performed to ensure that increasing these resolutions does not impact the results.

The horizontal Prandtl number $Pr_h = \mathcal {A}_h/\mathcal {D}_h$ is assumed to be equal to one throughout – Shakespeare & Hogg (Reference Shakespeare and Hogg2017) discuss the effect of non-zero Prandtl number on lee waves. Hereafter, we refer only to the viscosity $\mathcal {A}_h$, with the understanding that the diffusivity $\mathcal {D}_h$ varies similarly.