2.1 QG model equations in cylindrical coordinates

In this section we present the QG potential vorticity (PV) and energy equations for axially symmetric mean flow and bathymetry. Quasi-geostrophy is an approximation to fluid flow in a rotating reference frame, which is often a good approximation for synoptic scale oceanic flows (oceanic mesoscale), i.e. with characteristic length scales comparable with the Rossby radius of deformation, defined below (Pedlosky Reference Pedlosky1987). The necessarily small parameter in the approximation is the Rossby number $Ro=U/fX\ll 1$ , where $f$ is the Coriolis parameter, $U$ the velocity scale and $X$ the horizontal length scale. In these cases the Coriolis force approximately balances the pressure gradient, and to first order in $Ro$ , the evolution of the flow field is given by the QG PV equations. The QG approximation also requires the bathymetry and isopycnals to exhibit small variations relative to their respective domain-wide averages. While these conditions are not necessarily satisfied over continental slopes, previous studies suggest that QG captures the essential features of large-scale flows over topographic steepnesses typical of the ocean’s continental slopes (Williams, Read & Haine Reference Williams, Read and Haine2010; Stewart, Dellar & Johnson Reference Stewart, Dellar and Johnson2011, Reference Stewart, Dellar, Johnson, von Larcher and Williams2014; Poulin et al. Reference Poulin, Stegner, Hernández-Arencibia, Marrero-Díaz and Sangrà2014; Stern, Nadeau & Holland Reference Stern, Nadeau and Holland2015).

We use the $f$ -plane approximation (Pedlosky Reference Pedlosky1987), in which the reference frame revolves around the vertical axis with the same rate everywhere in the domain, neglecting the effect of the Earth’s curvature on the Coriolis acceleration. This isolates the effect of continental slope curvature, and thereby simplifies our analysis. This is partially justified by the fact that, dynamically, a topographic gradient induces a similar dynamical effect on rotating flow as does the latitudinal gradient of the rotation rate. This so-called topographic ${\it\beta}$ effect is usually much larger than the planetary ${\it\beta}$ effect in the local dynamics of slope-trapped currents.

We write the two-layer QG PV equations (Pedlosky Reference Pedlosky1964) in cylindrical coordinates,

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\partial }{\partial \tilde{t}}+\tilde{u} _{1r}\frac{\partial }{\partial \tilde{r}}+\frac{\tilde{u} _{1{\it\phi}}}{\tilde{r}}\frac{\partial }{\partial \tilde{{\it\phi}}}\right]\left[\tilde{{\rm\nabla}}^{2}\tilde{{\it\psi}}_{1}-\frac{1}{L_{1}^{2}}(\tilde{{\it\psi}}_{1}-\tilde{{\it\psi}}_{2})\right]=0, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\partial }{\partial \tilde{t}}+\tilde{u} _{2r}\frac{\partial }{\partial \tilde{r}}+\frac{\tilde{u} _{2{\it\phi}}}{\tilde{r}}\frac{\partial }{\partial \tilde{{\it\phi}}}\right]\left[\tilde{{\rm\nabla}}^{2}\tilde{{\it\psi}}_{2}-\frac{1}{L_{2}^{2}}\left(\tilde{{\it\psi}}_{2}-\tilde{{\it\psi}}_{1}-\frac{g^{\prime }}{f_{0}}\tilde{{\it\eta}}_{b}\right)\right]=0. & \displaystyle\end{eqnarray}$$

$\tilde{R}_{i}$

$\tilde{R}_{e}$

$\tilde{r}$

$\tilde{{\it\phi}}$

$\tilde{{\it\psi}}_{j}$

$(\tilde{u} _{jr},~\tilde{u} _{j{\it\phi}})=(-\tilde{r}^{-1}\partial \tilde{{\it\psi}}_{j}/\partial \tilde{{\it\phi}},~\partial \tilde{{\it\psi}}_{j}/\partial \tilde{r})$

$\tilde{{\it\zeta}}_{j}$

$\tilde{{\it\zeta}}_{j}=\tilde{{\rm\nabla}}^{2}\tilde{{\it\psi}}_{j}$

$g$

${\it\rho}_{j}$

$g^{\prime }=g(({\it\rho}_{2}-{\it\rho}_{1})/{\it\rho}_{1})$

$H_{j}$

$f_{0}$

$L_{j}=\sqrt{g^{\prime }H_{j}}/f_{0}$

$\tilde{{\it\eta}}_{b}(\tilde{r})$

$\partial \tilde{{\it\psi}}_{j}/\partial \tilde{{\it\phi}}=0$

$|_{\tilde{r}=\tilde{R}_{i},\tilde{R}_{e}}$

To study the instability of currents flowing parallel to the bathymetric isobaths, we assume a geostrophic, axially symmetric, azimuthal mean flow $\tilde{U} _{j{\it\phi}}(\tilde{r})=\partial \overline{\tilde{{\it\psi}}}_{j}/\partial \tilde{r}$ . This is an exact steady solution of (2.1a )–(2.1b ). We partition the streamfunction into mean and perturbation components, $\overline{\tilde{{\it\psi}}}_{j}$ and $\tilde{{\it\psi}}_{j}^{\prime }$ respectively. Linearizing the QG PV equations (2.1a )–(2.1b ) yields a linear system of equations for the perturbation streamfunctions,

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\partial }{\partial \tilde{t}}+\frac{\tilde{U} _{1{\it\phi}}}{\tilde{r}}\frac{\partial }{\partial \tilde{{\it\phi}}}\right]\left[\tilde{{\rm\nabla}}^{2}\tilde{{\it\psi}}_{1}^{\prime }-\frac{1}{L_{1}^{2}}(\tilde{{\it\psi}}_{1}^{\prime }-\tilde{{\it\psi}}_{2}^{\prime })\right]-\frac{1}{\tilde{r}}\frac{\partial \tilde{{\it\psi}}_{1}^{\prime }}{\partial \tilde{{\it\phi}}}\frac{\partial \tilde{Q}_{1}}{\partial \tilde{r}}=0, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\partial }{\partial \tilde{t}}+\frac{\tilde{U} _{2{\it\phi}}}{\tilde{r}}\frac{\partial }{\partial \tilde{{\it\phi}}}\right]\left[\tilde{{\rm\nabla}}^{2}\tilde{{\it\psi}}_{2}^{\prime }-\frac{1}{L_{2}^{2}}(\tilde{{\it\psi}}_{2}^{\prime }-\tilde{{\it\psi}}_{1}^{\prime })\right]-\frac{1}{\tilde{r}}\frac{\partial \tilde{{\it\psi}}_{2}^{\prime }}{\partial \tilde{{\it\phi}}}\frac{\partial \tilde{Q}_{2}}{\partial \tilde{r}}=0, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{Q}_{j}=\tilde{{\rm\nabla}}^{2}\overline{\tilde{{\it\psi}_{j}}}-\frac{1}{L_{j}^{2}}\left[(-1)^{j}(\overline{\tilde{{\it\psi}}_{2}}-\overline{\tilde{{\it\psi}}_{1}})-{\it\Delta}_{j2}\frac{g^{\prime }}{f_{0}}\tilde{{\it\eta}}_{b}\right]. & \displaystyle\end{eqnarray}$$

${\it\Delta}_{j2}=0,1$

$j=1,2$

The model describes a concave (convex) continental slope if $\tilde{{\it\eta}}_{b}(\tilde{r})$ is monotonically increasing (decreasing) with radius. A given convex (concave) along-slope flow can be transformed to the analogous concave (convex) along-slope flow by a radial reflection $P(\tilde{r}-\tilde{R}_{i})\rightarrow P(\tilde{R}_{e}-\tilde{r})$ , for any scalar radial property $P(\tilde{r}-\tilde{R}_{i})$ of the mean state, such as bathymetry $\tilde{{\it\eta}}_{b}(r)$ or isopycnal profile $\tilde{Z}_{I}(\tilde{r})$ .

The baroclinic growth rate in uniform rectilinear flow (Mechoso Reference Mechoso1980) peaks close to the wavenumber corresponding to the first baroclinic Rossby radius of deformation. Therefore we non-dimensionalize the equations by scaling $\tilde{r}\sim L$ , where

(2.3) $$\begin{eqnarray}\displaystyle L=\sqrt{\frac{g^{\prime }H_{1}H_{2}}{f_{0}^{2}(H_{1}+H_{2})}}. & & \displaystyle\end{eqnarray}$$

We denote the velocity scale (to be specified later) by $U$ . The non-dimensional variables are defined by

(2.4a-d ) $$\begin{eqnarray}\displaystyle t=(L/U)^{-1}\tilde{t},\quad U_{j{\it\phi}}=U^{-1}\tilde{U} _{j{\it\phi}},\quad {\it\eta}_{b}=(ULf_{0}/g^{\prime })^{-1}\tilde{{\it\eta}}_{b},\quad r=L^{-1}\tilde{r}. & & \displaystyle\end{eqnarray}$$

For notational convenience we also define $F_{j}=L^{2}/L_{j}^{2}=1-H_{j}/(H_{1}+H_{2})$ , which measures the fraction of the total depth that is not occupied by layer $j$ . Although $F_{1}$ and $F_{2}$ are not independent, we shall keep both parameters to preserve some symmetry in the presentation of the equations.

2.2 Method of solution

In what follows we drop the prime notation from the perturbation streamfunction for ease of presentation. The eigenvalue problem is derived by decomposing the perturbation streamfunction into normal azimuthal and temporal modes,

(2.5) $$\begin{eqnarray}\displaystyle {\it\psi}_{j}=\text{Re}\{{\it\Psi}_{j}(r)\exp (\text{i}(m{\it\phi}-{\it\sigma}t))\}. & & \displaystyle\end{eqnarray}$$

The notation $\text{Re}\{\}$ indicates the real part of the expression in the curly braces, and $i\equiv \sqrt{-1}$ . The azimuthal wavenumber is denoted as $m$ , and ${\it\sigma}$ is the non-dimensional complex frequency (dimensional $\tilde{{\it\sigma}}$ scales like $(U/L)$ by (2.4)). The real and imaginary parts of ${\it\sigma}$ are the frequency and the growth rate, respectively. The no normal flow boundary condition (stated above) simplifies to ${\it\Psi}_{j}|_{r=R_{i},R_{e}}=0$ . Writing ${\rm\nabla}_{r}^{2}=(\partial /\partial r+1/r)\partial /\partial r$ , the linear vorticity equations (2.2a )–(2.2b ) may be simplified as

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{U_{1{\it\phi}}}{r}m-{\it\sigma}\right]\left[{\rm\nabla}_{r}^{2}{\it\Psi}_{1}-\frac{m^{2}}{r^{2}}{\it\Psi}_{1}-F_{1}({\it\Psi}_{1}-{\it\Psi}_{2})\right]-\frac{m}{r}{\it\Psi}_{1}\frac{\partial Q_{1}}{\partial r}=0, & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{U_{2{\it\phi}}}{r}m-{\it\sigma}\right]\left[{\rm\nabla}_{r}^{2}{\it\Psi}_{2}-\frac{m^{2}}{r^{2}}{\it\Psi}_{2}-F_{2}({\it\Psi}_{2}-{\it\Psi}_{1})\right]-\frac{m}{r}{\it\Psi}_{2}\frac{\partial Q_{2}}{\partial r}=0, & \displaystyle\end{eqnarray}$$

(2.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial Q_{j}}{\partial r}=\frac{\partial }{\partial r}\left(\frac{\partial }{\partial r}+\frac{1}{r}\right)U_{j{\it\phi}}-F_{j}(-1)^{j}\left[\left(U_{2{\it\phi}}-U_{1{\it\phi}}\right)-{\it\Delta}_{j2}\frac{\partial {\it\eta}_{b}}{\partial r}\right]. & \displaystyle\end{eqnarray}$$

In most cases presented below we solve the eigenproblem posed by (2.6a )–(2.6c ) numerically. We discretize equations (2.6a )–(2.6b ) using second-order centred finite differences and solve the resulting matrix eigenvalue problem using the ‘eig’ function in Matlab, which uses the QZ algorithm (Moler & Stewart Reference Moler and Stewart1973). The grid resolution is $\text{d}r=0.025$ , giving 40 grid points per Rossby radii, thus resolving well the spatial scales normally associated with QG dynamics. Verification of the numerical set-up including convergence tests and comparison with some analytic results are presented in appendix A. The standard experiment parameters are: $F_{1}=F_{2}=1/2$ , $R_{i}=3$ , $R_{e}=10$ . The chosen channel width $(R_{e}-R_{i})$ is motivated by the widths of DWBCs, which are typically at least a few Rossby radii (Stommel & Arons Reference Stommel and Arons1972; Xu et al. Reference Xu, Rhines, Chassignet and Schmitz2015). Similar bathymetric curvature radii (in the range of 3–10 Rossby radii) are found around the GB and FC. Other parameter ranges and sensitivity tests are discussed in appendix A.

2.3 Energy equation

To study the modes of energy conversion from the mean state to perturbations, we derive the volume-integrated energy equation. The general method is standard (Pedlosky Reference Pedlosky1987): multiplying equations (2.2a ) and (2.2b ) by $D_{1}{\it\psi}_{1}$ and $D_{2}{\it\psi}_{2}$ respectively, adding the two resulting equations together, integrating in the entire domain, and using several integrations by parts and the no normal flow boundary conditions. We defined the relative layer thicknesses $D_{i}$ by $D_{1}=F_{2}=H_{1}/(H_{1}+H_{2})$ and $D_{2}=F_{1}=H_{2}/(H_{1}+H_{2})$ . In addition, one line integral over the domain boundaries, $\sum _{j=1}^{2}D_{j}\oint {\it\psi}_{j}^{\prime }(\partial ^{2}/\partial t\partial n){\it\psi}_{j}^{\prime }\,\text{d}s$ (where $n$ is the normal to the boundary), is required to vanish (McWilliams (Reference McWilliams1977), specifically equation (13)), for consistency with the analogous asymptotic expansion (in Rossby number) of the primitive equations’ energy balance. The derived energy equation in non-dimensional variables is

(2.7a ) $$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}E=\frac{\partial }{\partial t}\left\{\mathop{\sum }_{j=1}^{2}\text{EKE}_{j}+\text{EPE}\right\}=\mathop{\sum }_{j=1}^{2}\text{RS}_{j}+\text{PEC}, & & \displaystyle\end{eqnarray}$$

(2.7b,c ) $$\begin{eqnarray}\displaystyle \text{EKE}_{j}=\frac{1}{2}D_{j}\iint \left(\boldsymbol{{\rm\nabla}}{\it\psi}_{j}\right)^{2}r\,\text{d}r\,\text{d}{\it\theta},\quad \text{EPE}=\frac{1}{2}D_{1}D_{2}\iint ({\it\psi}_{1}-{\it\psi}_{2})^{2}r\,\text{d}r\,\text{d}{\it\phi},\quad & & \displaystyle\end{eqnarray}$$

(2.7d ) $$\begin{eqnarray}\displaystyle \text{RS}_{j}=D_{j}\iint \left(r\frac{\partial }{\partial r}\frac{U_{j{\it\phi}}}{r}\right)\left(\frac{1}{r}\frac{\partial {\it\psi}_{j}}{\partial {\it\phi}}\right)\left(\frac{\partial {\it\psi}_{j}}{\partial r}\right)r\,\text{d}r\,\text{d}{\it\phi}, & & \displaystyle\end{eqnarray}$$

(2.7e ) $$\begin{eqnarray}\displaystyle \text{PEC}=D_{1}D_{2}\iint (U_{1{\it\phi}}-U_{2{\it\phi}}){\it\psi}_{1}\frac{1}{r}\frac{\partial {\it\psi}_{2}}{\partial {\it\phi}}r\,\text{d}r\,\text{d}{\it\phi}. & & \displaystyle\end{eqnarray}$$

The energy of perturbations to the mean flow ( $E$ ) is a sum of the so-called eddy kinetic energy ( $\sum \text{EKE}_{j}$ ) and eddy potential energy (EPE). Thus energy tendency $\partial _{t}E$ is balanced by the volume-integrated Reynolds stress work ( $\sum \text{RS}_{j}$ ) and by potential energy conversion (PEC), i.e. conversion rates from mean kinetic and mean potential energy, respectively (Pedlosky Reference Pedlosky1987). When the net perturbation energy tendency (i.e. left-hand side of (2.7a )) is positive (i.e. perturbations grow), we may define a purely baroclinic instability as one where the Reynolds stress volume-integrated work is zero, as occurs in uniform rectilinear flow (Pedlosky Reference Pedlosky1987). We later show that when $\text{RS}_{j}$ do not vanish, they in fact are negative, i.e. they decrease the perturbation growth rate in all cases we study here.

Bathymetry does not enter the energy equation explicitly: it does contribute to energy exchange locally, but integrates to zero over the entire domain. The energy equation has zero energy tendency for an azimuthally constant perturbation, and therefore such perturbations are necessarily neutral. $\text{RS}_{j}$ are identically zero when the radial strain,

(2.8) $$\begin{eqnarray}\displaystyle S_{r}\equiv r\frac{\partial }{\partial r}\frac{U_{j{\it\phi}}}{r}, & & \displaystyle\end{eqnarray}$$

is identically zero, which in an annular channel occurs everywhere only for flow in solid-body rotation. Therefore solid-body rotation is the only annular flow that has zero Reynolds stress volume-integrated work for any infinitesimal perturbation. If $S_{r}$ is non-zero anywhere then there exist many particular ${\it\psi}(r,{\it\phi})$ perturbation shapes that make $\text{RS}_{j}$ non-zero.

2.4 Mean flow profiles

Throughout this paper we compare our results against the case of uniform flow in a straight channel over linear bathymetry (Pedlosky Reference Pedlosky1964; Mechoso Reference Mechoso1980), which is described in appendix B. We hereafter refer to this case as uniform rectilinear flow for short. In the annular channel, we investigate in detail two specific configurations of the bathymetry and the mean azimuthal flow. Since the mean flows we prescribe are geostrophic, the isopycnal profile $Z_{I}(r)$ is determined by the Margules relation (Cushman-Roisin Reference Cushman-Roisin1994). In dimensional variables,

(2.9) $$\begin{eqnarray}\displaystyle \tilde{U} _{1{\it\phi}}-\tilde{U} _{2{\it\phi}}=-\frac{g^{\prime }}{f}\frac{\partial \tilde{Z}_{I}}{\partial \tilde{r}}. & & \displaystyle\end{eqnarray}$$

The first case, solid-body rotation, is motivated by the fact that both uniform rectilinear flow and solid-body rotation have zero strain rate, defined for solid-body rotation by (2.8), and thus it is a simple starting point from which to study the effect of horizontal curvature. We assume parabolic bathymetry to simplify the analysis, though we later briefly explore linear bathymetry too (see § 6). Formally, we define our solid-body rotation case as

(2.10a-c ) $$\begin{eqnarray}\displaystyle U_{j{\it\phi}}={\it\Omega}_{j}r,\quad Z_{I}\sim -({\it\Omega}_{1}-{\it\Omega}_{2})r^{2},\quad {\it\eta}_{b}={\textstyle \frac{1}{2}}pr^{2}, & & \displaystyle\end{eqnarray}$$

where ${\it\Omega}_{j}$ are the constant angular velocities of the flow in each layer, and $p$ is a quadratic coefficient for the bathymetry.

The second case is uniform azimuthal flow, where we assume constant mean azimuthal velocity everywhere. This is similar to uniform rectilinear flow in that the speed is uniform, and the isopycnals are linear in the cross-flow coordinate ( $r$ ). It is different in that the velocity direction varies, i.e. speed is azimuthal everywhere but the azimuthal direction varies with the azimuthal angle ${\it\phi}$ . We take the bathymetry to be linear as well (as in the uniform rectilinear flow case), though we later briefly explore parabolic bathymetry as well (see § 6). Formally, we define our uniform azimuthal flow case as

(2.11a-c ) $$\begin{eqnarray}\displaystyle U_{j{\it\phi}}=\text{const.},\quad Z_{I}\sim -(U_{1{\it\phi}}-U_{2{\it\phi}})r,\quad {\it\eta}_{b}=br, & & \displaystyle\end{eqnarray}$$

where $U_{1{\it\phi}}$ and $U_{2{\it\phi}}$ are the azimuthal velocities and $b$ is the linear coefficient for the bathymetry.

We note that in uniform azimuthal velocity $R_{i}$ cannot be chosen to approach $r=0$ , both because the azimuthal velocity must be zero in the $r\rightarrow 0$ limit, and because even before the actual limit, the centrifugal force becomes larger than the Coriolis force, in violation of the QG conditions. The balance between the two forces results in a local Rossby number, $Ro=U/f\tilde{r}=U/(r\sqrt{g^{\prime }H_{1}H_{2}/(H_{1}+H_{2})})$ . For example, taking $H_{1}=H_{2}\approx 500~\text{m}$ , $g^{\prime }\approx 10^{-3}~\text{g}$ , and $U\approx 0.1~\text{m}~\text{s}^{-1}$ , we have $r>1$ (and $R_{i}>1$ ) as an approximate condition for $Ro=o(1)$ . Therefore our choice of $R_{i}=3$ (§ 2.2) is also consistent with the QG approximation.

We define the mean vertical rotation rate shear and velocity shear for solid-body rotation and for uniform azimuthal flow as follows: ${\it\Omega}_{s}={\it\Omega}_{1}-{\it\Omega}_{2}$ , and $U_{s}=U_{1{\it\phi}}-U_{2{\it\phi}}$ , respectively. Motivated by the fact that the baroclinic instability growth rate in uniform rectilinear flow varies linearly in the vertical velocity shear (Mechoso Reference Mechoso1980), we choose the velocity scale $U=L\tilde{{\it\Omega}}_{s}$ for solid-body rotation, and $U=\tilde{U} _{s}$ for uniform azimuthal flow. We assume everywhere that $U>0$ (and hence $Z_{I}(r)$ is monotonously decreasing). This assumption is general since we explore both positive and negative ${\it\delta}$ and can then deduce corresponding results for $U<0$ results by symmetry (see § 5.3).

Similarly we define for solid-body rotation and uniform azimuthal flow, the barotropic mean rotation rate and velocity as follows: ${\it\Omega}_{bt}=({\it\Omega}_{1}+{\it\Omega}_{2})/2$ , and $U_{bt}=(U_{1{\it\phi}}+U_{2{\it\phi}})/2$ , respectively. In fact ${\it\Omega}_{bt}$ (or $U_{bt}$ ) is exactly the barotropic component only if mean layer thicknesses are equal, but for ease of notation we refer to it as the barotropic component in what follows.

3.1 Derivation of the Rayleigh theorem

Our starting point is the modal PV equations (2.6a )–(2.6b ). We define the complex phase speed $c={\it\sigma}/m$ , and its real ( $c_{r}$ ) and imaginary ( $c_{i}$ ) parts. By (2.5), only unstable eigenmodes have a non-zero $c_{i}$ , so for unstable eigenmodes we may divide the equation for the layer $j$ by $m(U_{j{\it\phi}}/r-c)$ .

(3.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\rm\nabla}_{r}^{2}{\it\Psi}_{1}-\frac{m^{2}}{r^{2}}{\it\Psi}_{1}-F_{1}({\it\Psi}_{1}-{\it\Psi}_{2})-\frac{1}{U_{1{\it\phi}}-cr}{\it\Psi}_{1}\frac{\partial Q_{1}}{\partial r}=0, & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\rm\nabla}_{r}^{2}{\it\Psi}_{2}-\frac{m^{2}}{r^{2}}{\it\Psi}_{2}-F_{2}({\it\Psi}_{2}-{\it\Psi}_{1})-\frac{1}{U_{2{\it\phi}}-cr}{\it\Psi}_{2}\frac{\partial Q_{2}}{\partial r}=0. & \displaystyle\end{eqnarray}$$

$D_{1}{\it\Psi}_{1}^{\ast }$

$D_{2}{\it\Psi}_{2}^{\ast }$

$^{\ast }$

$(r\,\text{d}r)$

$R_{i}$

$R_{e}$

${\it\Psi}_{j}(R_{i})={\it\Psi}_{j}(R_{e})=0$

(3.2) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{R_{i}}^{R_{e}}\left[\mathop{\sum }_{j=1}^{2}D_{j}|\frac{\partial }{\partial r}{\it\Psi}_{j}|^{2}+\mathop{\sum }_{j=1}^{2}D_{j}\frac{m^{2}}{r^{2}}|{\it\Psi}_{j}|^{2}+D_{1}D_{2}|{\it\Psi}_{1}-{\it\Psi}_{2}|^{2}\right]r\,\text{d}r\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{R_{i}}^{R_{e}}\left[\mathop{\sum }_{j=1}^{2}\frac{D_{j}}{U_{j{\it\phi}}-cr}\frac{\partial Q_{j}}{\partial r}|{\it\Psi}_{j}|^{2}\right]r\,\text{d}r=0.\end{eqnarray}$$

The imaginary part of this expression is

(3.3) $$\begin{eqnarray}\displaystyle c_{i}\int _{R_{i}}^{R_{e}}\mathop{\sum }_{j=1}^{2}\frac{D_{j}}{|U_{j{\it\phi}}-cr|^{2}}\frac{\partial Q_{j}}{\partial r}|{\it\Psi}_{j}|^{2}r^{2}\,\text{d}r=0. & & \displaystyle\end{eqnarray}$$

For an unstable mode $c_{i}$ is non-zero, and so the last integral must vanish. Therefore a necessary condition for instability (hereafter, the Rayleigh criterion) is that the mean PV gradient must be somewhere negative and somewhere positive in the domain interior.

Another necessary condition for instability can be found using the real part of (3.2). Substituting (3.3) into (3.2) eliminates the terms proportional to $c_{r}$ and $c_{i}$ , leaving

(3.4) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{R_{i}}^{R_{e}}\mathop{\sum }_{j=1}^{2}\frac{D_{j}}{|U_{j{\it\phi}}-cr|^{2}}|{\it\Psi}_{j}|^{2}\left(U_{j{\it\phi}}\frac{\partial Q_{j}}{\partial r}\right)r\,\text{d}r\nonumber\\ \displaystyle & & \displaystyle \quad =-\int _{R_{i}}^{R_{e}}\left[\mathop{\sum }_{j=1}^{2}D_{j}\left|\frac{\partial }{\partial r}{\it\Psi}_{j}\right|^{2}+\mathop{\sum }_{j=1}^{2}D_{j}\frac{m^{2}}{r^{2}}|{\it\Psi}_{j}|^{2}+D_{1}D_{2}|{\it\Psi}_{1}-{\it\Psi}_{2}|^{2}\right]r\,\text{d}r<0.\end{eqnarray}$$

Therefore, another necessary instability condition (hereafter, the Fjortoft criterion) can be stated as: at least one of the products $U_{1{\it\phi}}(\partial Q_{1}/\partial r)$ and $U_{2{\it\phi}}(\partial Q_{2}/\partial r)$ must be negative inside at least part of the domain $(R_{i},R_{e})$ .

Both the first and second conditions as phrased here are the same as found in a straight channel (Pedlosky Reference Pedlosky1964). These (straight channel) conditions are used frequently to identify unstable flow regimes in boundary currents, as well as other ocean and atmosphere flow regimes.

3.2 Stability bounds for solid-body rotation and uniform azimuthal flow

For solid-body rotation, using the same notation as (Mechoso Reference Mechoso1980), the ratio between bathymetric slope and isopycnal slope is the bathymetric parameter ${\it\delta}=pr/{\it\Omega}_{s}r=p/{\it\Omega}_{s}$ . Vorticity is constant and hence the PV gradients are simply

(3.5a,b ) $$\begin{eqnarray}\displaystyle \frac{\partial Q_{1}}{\partial r}=-F_{1}r,\quad \frac{\partial Q_{2}}{\partial r}=F_{2}(1-{\it\delta})r. & & \displaystyle\end{eqnarray}$$

By the Rayleigh criterion instability is possible only if the PV gradient changes sign, which is seen from (3.5) to occur only if ${\it\delta}<1$ , exactly as in uniform rectilinear flow.

For uniform azimuthal flow, in non-dimensional variables, $U_{1{\it\phi}}=U_{2{\it\phi}}+1$ . The bathymetric parameter is now ${\it\delta}=-b$ , and from (2.6c ),

(3.6a,b ) $$\begin{eqnarray}\displaystyle \frac{\partial Q_{1}}{\partial r}=-\frac{1}{r^{2}}\left(U_{2{\it\phi}}+1\right)-F_{1},\quad \frac{\partial Q_{2}}{\partial r}=-\frac{1}{r^{2}}U_{2{\it\phi}}+F_{2}(1-{\it\delta}). & & \displaystyle\end{eqnarray}$$

Using (3.6) in the instability criteria (§ 3.1), it follows that a necessary condition for instability is that, at least somewhere inside the domain,

(3.7) $$\begin{eqnarray}\displaystyle {\it\delta}<{\it\delta}_{0}(r)\equiv 1-\frac{1}{F_{2}}\frac{U_{2{\it\phi}}}{r^{2}}. & & \displaystyle\end{eqnarray}$$

Equivalently, a sufficient condition for stability is ${\it\delta}>\max ({\it\delta}_{0})$ , similarly to the uniform rectilinear flow case which is stable for ${\it\delta}>1$ . Note that ${\it\delta}_{0}=1$ exactly if $U_{2{\it\phi}}=0$ , and ${\it\delta}_{0}\approx 1$ if $(1/F_{2})|U_{2{\it\phi}}(r)|/r^{2}=o(1)$ . Unlike the case in uniform rectilinear flow, the stability threshold depends on the mean velocity magnitude, i.e. the cutoff bathymetric parameter, ${\it\delta}_{0}$ , increases (decreases) for negative (positive) $U_{2{\it\phi}}$ . As evident from the Rayleigh criterion and from (3.6), the difference is due to the non-zero mean relative vorticity caused by curved streamlines.

We add that for $U_{2{\it\phi}}<-1-F_{1}R_{i}^{2}$ , instability is not prohibited, irrespective of the ${\it\delta}$ value, since a PV gradient sign change occurs within a single layer. However, since the mean flow becomes almost barotropic, that regime is less relevant to this study.

5.1 Normal modes

In figure 3, we plot growth rate (GR) as a function of normalized wavenumber $\hat{m}$ and of ${\it\delta}$ for uniform azimuthal flow (each panel for a different $U_{bt}$ value), and for reference also the GR of the uniform rectilinear flow and solid-body rotation cases. Note that at each point in $(\hat{m},{\it\delta})$ space there may be multiple unstable modes, so we have plotted the growth rate of the most unstable mode in each case. While for zero barotropic flow $U_{bt}=0$ the growth rate is similar to uniform rectilinear flow, non-zero barotropic velocity results in very different $\text{GR}(\hat{m},{\it\delta})$ dependence. In contrast, uniform rectilinear flow and solid-body rotation have no barotropic velocity dependence. Additional local maxima in $\text{GR}(\hat{m},{\it\delta})$ appear in uniform azimuthal flow for non-zero barotropic velocity.

Figure 4.

Mean uniform azimuthal flow, selected unstable eigenmodes. Upper (lower) layer streamfunctions are shown in the upper (lower) panels. The bathymetric slope parameter is ${\it\delta}=-0.2$ , and the azimuthal wavenumber is $m=4$ . (a,b) $U_{bt}=0$ , fastest growing eigenmode. (c,d) $U_{bt}=-1$ , fastest growing eigenmode. (e,f) $U_{bt}=-1$ , third-fastest growing eigenmode. The inner and outer circles mark the domain boundaries at $r=R_{i}$ and $R_{e}$ , respectively. The lines intercepting the boundaries are the zero contours of the streamfunctions, while positive (negative) streamfunction contours are denoted by full (dashed) closed curves. The absolute value of the contours is not given since eigenmode amplitudes are arbitrary unless specified by initial conditions.

The streamfunctions for several unstable uniform azimuthal flow modes are presented in figure 4. Two geometrical differences from uniform rectilinear flow and solid-body rotation (see examples in figure 2) are evident: (i) while in solid-body rotation streamfunctions are always centred in the channel, the streamfunctions in uniform azimuthal flow cases with non-zero $U_{bt}$ are shifted off the centre of the channel; and (ii) while in uniform rectilinear flow the streamfunction circulation cell axes are aligned with the radial direction, those in uniform azimuthal flow cases with non-zero $U_{bt}$ are often tilted. Reynolds stress work varies linearly with the strain and the tilts of the circulation axes, and vanishes when the tilt is zero (Pedlosky Reference Pedlosky1987). In polar coordinates,

(5.1) $$\begin{eqnarray}\displaystyle \text{RS}_{j}\sim -(S_{r})_{j}\left(\frac{\partial r}{\partial {\it\phi}}\right)_{{\it\psi}_{j}}. & & \displaystyle\end{eqnarray}$$

The uniform azimuthal flow streamfunction for the first mode at $U_{bt}=0$ (figure 4 a,b) has zero or very small tilt, implying insignificant Reynolds stress work. Figure 4(c–f) (for the first and third eigenmodes with $U_{bt}=-1$ ) show progressively higher positive tilts, implying higher negative $\text{RS}_{j}$ , since by (2.8) strain rate is positive for constant negative velocity.

The local maxima in figure 3(d–f) are due to changes in the number and character of growing eigenmodes with $U_{bt}$ value. This can be seen in figure 5, where we plot several properties for all uniform azimuthal flow growing eigenmodes at ${\it\delta}=-0.2$ and $U_{bt}=1$ (compare with solid-body rotation and uniform rectilinear flow, figure 2). We find that up to four unstable eigenmodes can co-exist at a given wavenumber, whereas no more than two co-existed for uniform rectilinear flow and solid-body rotation. While in uniform rectilinear flow the second mode has considerably lower growth rate than the first, in uniform azimuthal flow they have similar maximum values but still peak at different wavenumbers, thus explaining the multiple maxima observed in figure 3(d–f). The growth rates of the third- and fourth-most unstable modes are considerably smaller. In figure 5(c) we plot the eigenmodes’ phase speeds, Doppler-corrected and normalized via ${\hat{c}}_{r}=c_{r}R-U_{bt}$ to compare approximately with equivalent values in uniform rectilinear flow. While uniform rectilinear flow has waves propagating with shallow water to their right (prograde), in uniform azimuthal flow the second most rapidly growing mode is retrograde, and eigenmodes 3 and 4 have much smaller propagation speeds ${\hat{c}}_{r}$ .

Figure 5(b) shows the ratios of volume-integrated Reynolds stress work ( ${\it\Sigma}_{j}\text{RS}_{j}$ ) to potential energy conversion (PEC), which are negative in all cases. In contrast, uniform rectilinear flow has zero $\text{RS}_{j}$ values in all cases. Thus the uniform azimuthal flow unstable eigenmodes are largely baroclinic modes whose growth rates are somewhat diminished by Reynolds stress work. The two extra eigenmodes that appear with non-zero barotropic velocity have much higher ${\it\Sigma}_{j}\text{RS}_{j}$ to PEC ratio magnitudes, consistent with their very low growth rates. These results are also consistent with the tilts of streamfunctions shown in figure 4, and are qualitatively similar for other values of ${\it\delta}$ , $m$ and $U_{bt}$ . The general reduction in growth rate with $|U_{bt}|$ is thus partially attributed to an increase in Reynolds stress work.

Figure 5.

Mean uniform azimuthal flow with barotropic velocity $U_{bt}=1$ and bathymetric to isopycnal slopes ratio ${\it\delta}=-0.2$ . (a) Growth rate, (b) ratio of Reynolds stress volume-integrated work to potential energy conversion, and (c) Doppler-corrected Cartesian phase speed versus normalized azimuthal wavenumber $\hat{m}=m/R$ , for all growing eigenmodes. In (c) the (real) phase speed is Doppler-corrected and normalized to Cartesian values by ${\hat{c}}_{r}=c_{r}R-U_{bt}$ , to facilitate comparison with the other mean flow cases (figure 2).

Figure 6.

(a) Maximum growth rate (filled contours) in the uniform azimuthal flow (UAF) case as a function of the barotropic velocity, $U_{bt}$ , and the ratio of the bathymetric to isopycnal slopes, ${\it\delta}$ . The dashed line marks the barotropic velocity corresponding to the largest growth rates at each ${\it\delta}$ . The dotted line marks ${\it\delta}=1$ , above which straight channel uniform flow is stable. (b) Maximum uniform azimuthal flow growth rate (full line) and the barotropic velocity at which it is achieved (full line with circles), as a function of ${\it\delta}$ . The dashed line is maximum growth rate for uniform rectilinear flow (URF). (c) (Half the value of) potential energy conversion (PEC) and Reynolds stress work ( $\text{RS}={\it\Sigma}_{j}\text{RS}_{j}$ ), for the fastest growing eigenmode, in three different ${\it\delta}$ values. Discontinuities (as a function of $U_{bt}$ ) are expected since $\text{PEC}{-}\text{RS}$ distribution changes with $m$ , and since up to four different eigenmodes exist per $m$ .

In figure 6 (a,b) we plot the maximum growth rate over all unstable modes and over all wavenumbers as a function of $U_{bt}$ and ${\it\delta}$ . Unless ${\it\delta}\sim 1$ , the growth rate peaks at or close to $U_{bt}=0$ , and is close to peak growth rate for the uniform rectilinear flow case, while lower growth rates are found for non-zero $U_{bt}$ . However, the uniform rectilinear flow instability has a cutoff at ${\it\delta}=1$ , whereas the uniform azimuthal flow cutoff depends on $U_{bt}$ and can occur for ${\it\delta}$ larger than 1, as predicted by the instability criteria derived in § 3.2. Therefore horizontal curvature decreases the eigenmodes’ growth rates when $U_{bt}$ is non-zero, unless the isopycnals are approximately parallel to the bathymetry ( ${\it\delta}\sim 1$ ) and $U_{bt}<0$ , in which case the curvature destabilizes the flow. Although we report above that Reynolds stress work is partially responsible for the reduction in growth rate (both relative to uniform rectilinear flow and between different uniform azimuthal flow modes), figure 6(c) demonstrates that a reduction in potential energy conversion is responsible for a ${\sim}2{-}4$ times larger fraction of the growth rate reduction than is the $|{\it\Sigma}_{j}\text{RS}_{j}|$ increase. For an eigenmode, $GR=(\text{PEC}+{\it\Sigma}_{j}\text{RS}_{j})/2$ and hence the changes in growth rate are proportional to changes in PEC and $\text{RS}_{j}$ . Generally PEC decreases monotonously with $|U_{bt}|$ , thus supporting the barotropic governor effect interpretation given below.

Reduction in baroclinic growth rate in the presence of lateral barotropic shear is a somewhat general phenomenon, often called the barotropic governor effect (James & Gray Reference James and Gray1986; James Reference James1987). James (Reference James1987) attributes the effect to the horizontal shear of advection. To remain in phase in the cross-flow direction in the presence of advective shear, the unstable eigenmodes are tilted in the horizontal plane and have a reduced cross-flow extent. These circulation features make the unstable eigenmodes less ideally suited for extracting mean potential energy and hence their growth rates are smaller.

Though the usual barotropic governor effect interpretation is due to barotropic shear, we find that in curved flow geometry, it may be more general to refer to barotropic strain rather than shear. In uniform azimuthal flow, the eigenmodes are azimuthally travelling waves, with constant angular phase velocity $(\text{rad}~\text{s}^{-1})$ . If they were to have no tilt, azimuthal advection would need to be radially constant. The azimuthally advecting quantity is the (mean) angular velocity, ${\it\Omega}_{j}(r)\equiv U_{j{\it\phi}}/r$ , and its radial gradient is the strain rate (2.8). Thus the barotropic governor effect can generally occur in azimuthal flow with non-zero strain, i.e. flow not in solid-body rotation (§ 4). In accordance, we find (figure 4) that some uniform azimuthal flow eigenmodes have substantial horizontal tilts, often with much lower growth rates. Similar results were obtained in a primitive-equations two-layer vortex instability model (Dewar & Killworth Reference Dewar and Killworth1995). The authors found that Gaussian vortices with co-rotating lower layers had reduced PEC and growth rates relative to vortices with counter-rotating lower layers, and attributed the result to the barotropic governor effect.

5.2 Non-normal growth

While the above diagnosis focuses on perturbation growth by individual eigenmodes (aka normal modes), non-orthogonality of eigenmodes, which is a common occurrence in sheared flow, means that what is known as non-normal growth is also possible (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Farrell & Ioannou Reference Farrell and Ioannou1996). Linear evolution of two or more non-orthogonal eigenmodes, even if they are all neutral or decaying, can result in transient (non-normal) growth before the eventual decay of the disturbance. For parameter values where growing eigenmodes exist, they do dominate the linearized dynamics, at long enough times. But transient non-normal growth may dominate at shorter times, as well as for parameter ranges where no growing eigenmodes exist.

To calculate the non-normal growth, using the same numerical eigenvalue solver described above, we recast the PV equations (2.6) in the form: $M{\it\Psi}={\it\sigma}B{\it\Psi}$ , where $M$ and $B$ are differential operators, and ${\it\Psi}=({\it\Psi}_{1},{\it\Psi}_{2})^{\text{T}}$ ( $\text{T}$ for transpose). We refer the reader to Farrell & Ioannou (Reference Farrell and Ioannou1996) for details of the method. Assuming $B$ is invertible we can rewrite the differential equation as ${\it\sigma}{\it\Psi}=L{\it\Psi}$ , where $L=B^{-1}M$ . And since ${\it\Psi}\sim \exp (-\text{i}{\it\sigma}t)$ ,

(5.2) $$\begin{eqnarray}\displaystyle \partial _{t}{\it\Psi}=-\text{i}L{\it\Psi}. & & \displaystyle\end{eqnarray}$$

So the propagator to time $t$ is $\exp (-\text{i}Lt)$ . If we define $\hat{L}=N^{1/2}LN^{-1/2}$ , where $N$ is the energy norm operator (Farrell & Ioannou Reference Farrell and Ioannou1996), then the maximal instantaneous growth-rate of disturbances is given by the eigenvalues (and eigenstates) of the operator $H=\text{i}(\hat{L}^{\dagger }-\hat{L})/2$ . The energy-norm operator in the annulus case is, prior to the performed discretization of the differential operators and of $r$ ,

(5.3) $$\begin{eqnarray}\displaystyle N=\frac{1}{2}\left(\begin{array}{@{}cc@{}}rD_{1} & 0\\ 0 & rD_{2}\end{array}\right)\left(-{\rm\nabla}_{r}^{2}+\frac{m^{2}}{r^{2}}\right)+\frac{1}{2}D_{1}D_{2}\left(\begin{array}{@{}cc@{}}1 & -1\\ -1 & 1\end{array}\right). & & \displaystyle\end{eqnarray}$$

In figure 7 the maximal instantaneous non-normal growth in energy norm is shown for uniform azimuthal flow and uniform rectilinear flow. Interestingly, the result is independent of ${\it\delta}$ , as the bathymetry does not appear in the energy equation (2.7). Bathymetry affects local energy conversion but not its domain integral. However, in finite time bathymetry certainly effects energy growth or decay since it affects the streamfunction evolution, and would likely render sub-optimal the fastest-growing non-normal perturbations calculated using (5.3). For $U_{bt}=0$ , uniform azimuthal flow growth is very similar to uniform rectilinear flow (which is independent of barotropic velocity), and both have non-normal growth just slightly higher than peak normal growth rate (compare with figure 6). However, non-normal growth occurs for a wider range of wavenumbers than the range in which unstable normal modes exist. The decay of the non-normal growth rate with wavenumber is slower in uniform azimuthal flow relative to uniform rectilinear flow. In addition, at non-zero $U_{bt}$ , the uniform azimuthal flow growth rate is higher everywhere, and decays even slower with $\hat{m}$ , or even (not shown) grows and oscillates in $\hat{m}$ before decaying again. Growth at very high wavenumbers is probably not physical, and would likely not appear if some form of scale-selective ‘eddy’ viscosity were included.

Figure 7.

Maximal instantaneous (non-normal) growth rates. The blue curve corresponds to mean uniform rectilinear flow. The other curves correspond to mean uniform azimuthal flow with three different barotropic velocities $U_{bt}$ . Linear bathymetry was used in all cases.

5.3 Convex and concave cases

The convex to concave transformation (by reflection of ${\it\eta}_{b}(r)$ and $Z_{I}(r)$ , as explained in § 2.1), results for uniform azimuthal flow in $b\rightarrow -b$ and $U_{j{\it\phi}}\rightarrow -U_{j{\it\phi}}$ . Therefore (2.6) are unaltered if in addition ${\it\sigma}\rightarrow -{\it\sigma}$ . Also, if $\{{\it\sigma},{\it\Psi}_{j}(r)\}$ are an eigenvalue–eigenfunction pair of equations (2.6) then so are their complex conjugates $\{{\it\sigma}^{\ast },{\it\Psi}_{j}^{\ast }(r)\}$ , as can be verified by taking the complex conjugate of (2.6). Combining the last two observations, if $\{{\it\sigma},{\it\Psi}_{j}(r)\}$ is an eigenvalue–eigenfunction pair in a convex geometry, then $\{-{\it\sigma}^{\ast },{\it\Psi}_{j}^{\ast }(r)\}$ is an eigenvalue–eigenfunction pair in a concave geometry, and vice versa. Thus the growth rate and the real (physical) part of the streamfunction are unaltered, and the phase speed is reversed. The reversal of phase speed, along with reflection of ${\it\eta}_{b}(r)$ (shallow water at the other side of the channel) results in the same phase speed direction relative to shallow water, and therefore the physical propagation direction is also unaltered.

Different eddy growth rate at convex versus concave sections is not accounted for in the linear uniform azimuthal flow model, contrasting with the pronounced instabilities observed around convex bends in the real ocean’s continental slopes (see § 1). Despite this symmetry of our QG model, there is still a potentially important difference between unstable modes growing over convex and concave bends in the continental slope. If the perturbation streamfunction is displaced off the centre of the channel, say toward shallower water, then switching between concave and convex continental slopes will result in the streamfunction being displaced toward deep water instead. As noted in § 5.1, the uniform azimuthal flow perturbation streamfunctions are indeed typically displaced from the centre of the channel (see figure 4). The significance of this difference between the structure of growing modes over convex and concave continental slopes cannot be determined from our linear instability analysis, and warrants further investigation using a nonlinear model.