2.1 General linear theory à la Solberg

In this exploratory study of near-equatorial inertial instability, we shall only consider the symmetric inertial instability. As noted above, this implies conservation of momentum. Much insight is obtained with the heuristic argument of Solberg (Reference Solberg1936), which is easily adapted to our system (the ‘parcel method’ of Solberg in meteorology is further discussed in Holton (Reference Holton1992)). It is convenient to write the equations (2.5b,c ) for the meridional motions $\boldsymbol{v}=v\boldsymbol{j}+w\boldsymbol{k}$ as

(a constant term has been absorbed by the pressure). For a steady zonal flow $\boldsymbol{u}=U\boldsymbol{i}$ , we will indicate the associated momentum with

The balance is described by

(with $P$ the pressure field and $\bar{\unicode[STIX]{x1D70C}}$ the density field), and the thermal wind balance is

Solberg’s argument is a Lagrangian argument with the velocity $\boldsymbol{v}$ seen as the time derivative of a meridional displacement vector $\unicode[STIX]{x1D6FF}\boldsymbol{x}(t)$ , which we write as $\boldsymbol{v}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6FF}\boldsymbol{x}$ . Under small displacements, to lowest order, conservation of momentum and density are

with $\unicode[STIX]{x1D6FF}m$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}$ the local perturbations. If we denote the pressure perturbation with respect to the balance pressure $P$ with $p^{\prime }$ , the equations of motion (2.10) are

Subtracting the balanced part (2.12) and substituting from (2.14), one finds the linearized equations

or in vector–matrix notation

with

Here, $N^{2}$ is the usual square of the buoyancy frequency. In view of the thermal wind balance (2.13), $\unicode[STIX]{x1D648}=\left(\!\begin{smallmatrix}A & -B\\ -B & C\end{smallmatrix}\!\right)$ is a symmetric $2\times 2$ matrix. Thus, the eigenvalues $\{\unicode[STIX]{x1D706}_{1},\unicode[STIX]{x1D706}_{2}\}$ are real and the corresponding eigenvectors $\{\boldsymbol{e}_{1},\boldsymbol{e}_{2}\}$ are orthogonal.

Solberg (Reference Solberg1936) argued that the flow is unstable if the displaced element accelerates away from its original position; if ‘pushed back’, it is stable. In this sense, there is stability if both eigenvalues of $\unicode[STIX]{x1D648}$ are positive because then (2.17) with $p^{\prime }=0$ is a 2D harmonic oscillator equation. Instability follows if an eigenvalue $\unicode[STIX]{x1D706}_{i}<0$ because then a displacement $\unicode[STIX]{x1D6FF}x_{i}$ in the $\boldsymbol{e}_{i}$ -direction would amplify in time with the exponential rate $\exp (\sqrt{-\unicode[STIX]{x1D706}_{i}}t)$ . Excluding convective instability due to unstable stratification ( $N^{2}<0$ ), inertial instability is then associated with just one negative eigenvalue, say $\unicode[STIX]{x1D706}_{1}<0$ , and the preferred direction of the instability is along $\boldsymbol{e}_{1}$ , i.e. displacements in this direction exhibit maximal growth.

Refinements of this simple argument were developed by Ooyama (Reference Ooyama1966), who took into account the pressure perturbations. He proved stability/instability with an energy norm, as follows with the signs of $\unicode[STIX]{x1D706}_{1}(\boldsymbol{x})$ and $\unicode[STIX]{x1D706}_{2}(\boldsymbol{x})$ . Earlier, Fjørtoft (Reference Fjørtoft1950) arrived at the stability conditions with normal-mode analysis, but Ooyama pointed out that often normal modes cannot exist in bounded domains, as was first noted by Høiland (Reference Høiland1962). Related approaches, i.e. stability proofs through the use of invariant quadratic functionals, are found in Kloosterziel & Carnevale (Reference Kloosterziel and Carnevale2007) and in Fruman & Shepherd (Reference Fruman and Shepherd2008) for compressible flows.

Since much depends on the (spatial) properties of the eigenvalues of $\unicode[STIX]{x1D648}$ , it should be noted that

where $\boldsymbol{Q}$ is the absolute vorticity vector and $Q_{E}$ is the Ertel potential vorticity, i.e.

With the traditional approximation (TA) in the lower row of $\unicode[STIX]{x1D648}$ , the terms multiplied by $2\unicode[STIX]{x1D6FA}$ need to be discarded and $M=M_{TA}=U-(1/2)\unicode[STIX]{x1D6FD}y^{2}$ .

If Rayleigh’s argument mentioned in the introduction is extended to include stratification, i.e. if potential energy changes are included, the coefficients of $\unicode[STIX]{x1D648}$ appear as follows. If fluid elements located at $\boldsymbol{x}-(\unicode[STIX]{x1D6FF}\boldsymbol{x})/2$ , $\boldsymbol{x}+(\unicode[STIX]{x1D6FF}\boldsymbol{x})/2$ exchange positions, the energy change is $\unicode[STIX]{x0394}E=\unicode[STIX]{x1D6FF}\boldsymbol{x}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{x}$ . Diagonalized, this is $\unicode[STIX]{x0394}E=\unicode[STIX]{x1D706}_{1}\unicode[STIX]{x1D6FF}x_{1}^{2}+\unicode[STIX]{x1D706}_{2}\unicode[STIX]{x1D6FF}x_{2}^{2}$ . For infinitesimal $\unicode[STIX]{x1D6FF}\boldsymbol{x}$ , $\unicode[STIX]{x0394}E<0$ when, for example, $\unicode[STIX]{x1D706}_{1}<0$ and $\unicode[STIX]{x1D6FF}x_{1}\neq 0$ . If, additionally, $\unicode[STIX]{x1D706}_{2}\geqslant 0$ , maximal energy release is found for an exchange along the direction defined by $\unicode[STIX]{x1D706}_{1}$ , i.e. along $\boldsymbol{e}_{1}$ . The orientation of the secondary motions sketched in figure 1(a) in circular flows without stratification is explained with the far more sophisticated approach of Ooyama (Reference Ooyama1966). He showed that if, in some region, $\unicode[STIX]{x1D706}_{1}(\boldsymbol{x})<0$ , an initial field of displacements $\unicode[STIX]{x1D6FF}\boldsymbol{x}(\boldsymbol{x})$ can be chosen that leads to unbounded growth. For this, the displacements have to be concentrated within sufficiently elongated elliptically shaped subdomains, with the major axis of the ellipses along the ‘instability’ direction locally defined by $\boldsymbol{e}_{1}$ wherever $\unicode[STIX]{x1D706}_{1}(\boldsymbol{x})<0$ . Early numerical experiments by Yanai & Tokiaka (Reference Yanai and Tokiaka1969) proved the validity of Ooyama’s prediction.

Figure 3. Side view of a slice of the spherical planet rotating with angular velocity $\unicode[STIX]{x1D6FA}$ . The distance from the rotation axis is $R$ . The local Cartesian coordinates $y,z$ are indicated. With the complete Coriolis force (NTA), Taylor columns $U=U(R)$ are invariant along lines of constant distance $R\approx R_{0}+z-y^{2}/2R_{0}$ .

2.2 Homogeneous fluid

In a homogeneous fluid (constant density $\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{0}$ ), steady zonal flows $\boldsymbol{u}=U\boldsymbol{i}$ have to be Taylor columns according to the Taylor–Proudman theorem; that is, flows that do not vary parallel to the axis of rotation: $(\unicode[STIX]{x1D734}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=0$ (Pedlosky Reference Pedlosky1987). This means with (2.6) that only $U=U(R)=U(z-y^{2}/2R_{0})$ can be steady. This is sketched in figure 3. The stability tensor $\unicode[STIX]{x1D648}$ reduces to

with eigenvalues and eigenvectors

having left out a common normalizing factor in $\{\boldsymbol{e}_{1},\boldsymbol{e}_{2}\}$ . Specifically,

with $R$ approximated by the parabolas (2.6). The direction $\boldsymbol{e}_{2}$ associated with $\unicode[STIX]{x1D706}_{2}=0$ is everywhere tangent to $R=\text{const.}$ , i.e. along the direction of $\unicode[STIX]{x1D734}$ . This ‘neutrally stable’ direction is parallel to the orientation of the Taylor column sketched in figure 3 and is the direction in which in Rayleigh’s argument an exchange does not liberate energy. The direction $\boldsymbol{e}_{1}$ is perpendicular to the lines of constant $R$ . The criterion for instability ( $\unicode[STIX]{x1D706}_{1}<0$ ) is nothing but Rayleigh’s criterion (1.1) with the approximations we have made: as noted above after (2.9), $L\approx R_{0}M$ and in (1.1)

because $2L/R^{2}\approx 2\unicode[STIX]{x1D6FA}$ if $|U|\ll \unicode[STIX]{x1D6FA}R_{0}$ (on the Earth, $\unicode[STIX]{x1D6FA}R_{0}\approx 500~\text{m}~\text{s}^{-1}$ ).

When $\unicode[STIX]{x1D706}_{1}(\boldsymbol{x})<0$ in some region, the flow is unstable and the preferred direction of the instability is along $\boldsymbol{e}_{1}$ , which near the equator is almost perpendicular to the surface, i.e. upward/downward, parallel to lines of increasing depth $z$ (as sketched in figure 2 b). We believe that this has rarely been investigated in any detail.

2.3 Model flow

Within a homogeneous fluid, with the NTA, only flows $U(z-y^{2}/2R_{0})$ are steady (see § 2.2). These are flows $U(R)$ with $R$ approximated by (2.6). We shall use the simple example

From here on, $L$ defines the latitudinal width of the jet – not the angular momentum – and $U_{0}$ is the peak surface velocity at the equator. At the surface ( $z=0$ ), this flow matches the Gaussian jet we recently studied in the context of the TA approximation (see Kloosterziel, Orlandi & Carnevale Reference Kloosterziel, Orlandi and Carnevale2015). However, then a balanced flow has to be a Taylor column $U(y)$ and there is no $z$ -dependence. With the NTA, the amplitude has to decay with depth below the surface, $z<0$ , with a vertical scale $L^{2}/R_{0}$ . If $R_{0}\approx 6380$  km (Earth’s radius), for $L=25$  km, the decay scale is a mere 100 m, for $L=50$  km, approximately 400 m, etc.

Only the westward flowing current can be unstable ( $U_{0}<0$ ). With non-dimensional $\tilde{z}=zR_{0}/L^{2}$ , ${\tilde{y}}=y/L$ , it follows with (2.23) that

For instability ( $\unicode[STIX]{x1D706}_{1}<0$ ), it is required that $A>1$ ( $A$ is an equatorial Rossby number). With $\unicode[STIX]{x1D6FA}=7.26\times 10^{-5}~\text{rad}~\text{s}^{-1}$ (Earth’s radius angular velocity), this requires a westward peak velocity $|U_{0}|\gtrsim 1.5~\text{cm}~\text{s}^{-1}$ if $L=25$  km, if $L=50$  km, $|U_{0}|\gtrsim 6~\text{cm}~\text{s}^{-1}$ , and so on.

Figure 4. Schematic illustrating the prediction of absolute momentum mixing. In (a), a Taylor-column flow as sketched in figure 3 is assumed, which is unstable for $R>R_{ins}$ where absolute momentum decreases: $\text{d}M/\text{d}R<0$ (dark grey area). For smaller $R$ , the flow is stable (uncoloured area). In the unstable (dark grey) area, the overturning motions as sketched in figure 2 are expected to emerge. At high Reynolds numbers, through nonlinear interactions, these secondary motions are expected to propagate towards the rotation axis (indicated by fat horizontal black arrows). By doing so, they advect and mix momentum. In (b), we show the expected final state. A new Taylor columnar flow is established with constant absolute momentum $m=M_{c}$ between the axial radius $R=R_{l}$ and the surface (light grey area), while the momentum is unchanged for $R<R_{l}$ (uncoloured area). The flow is stable because everywhere $\text{d}M/\text{d}R\geqslant 0$ . Continuity at $R=R_{l}$ is established by setting $M_{c}=M(R_{l})$ , and $M_{c}$ , and then also $R_{l}$ , is determined by conservation of total (area integrated) momentum according to (2.32).

If exponential growth of small perturbations is assumed, i.e. $\exp (\unicode[STIX]{x1D6FE}t)$ , the initial growth rate $\unicode[STIX]{x1D6FE}$ of a displacement $\unicode[STIX]{x1D6FF}x_{1}$ is via Solberg’s approach locally $\unicode[STIX]{x1D6FE}=\sqrt{-\unicode[STIX]{x1D706}_{1}}$ , where $\unicode[STIX]{x1D706}_{1}<0$ . The instability region is contained in between the parabola $\tilde{z}-{\tilde{y}}^{2}/2=-\text{ln}\,A$ and the surface $\tilde{z}=0$ . The lowest point on the instability boundary parabola is at depth $\tilde{z}_{ins}=-\text{ln}\,A$ and it meets the surface at ${\tilde{y}}_{\pm }=\pm \sqrt{2\ln A}$ . With the approximation (2.6), the boundary is the line of constant distance $R=R_{0}-(L^{2}/R_{0})\ln A$ from the rotation axis. For smaller $R$ , one has $\text{d}M/\text{d}R>0$ , and for larger $R$ , the instability region is where $\text{d}M/\text{d}R<0$ (see (2.23) again and figure 4 a). The fastest rate of growth, $\unicode[STIX]{x1D6FE}_{\infty }$ , is at the equator at the surface, i.e. at $\{{\tilde{y}},\tilde{z}\}=\{0,0\}$ . Everywhere, the instability direction is along $\boldsymbol{e}_{1}$ , i.e. radial or normal to lines of constant $R$ (see the end of § 2.1). If time is scaled with $T=1/2\unicode[STIX]{x1D6FA}$ , non-dimensionally,

With the TA, the flow that has the same surface signature as (2.25) is

which is also only unstable when flowing westward with $A>1$ . The eigenvalues and eigenvectors of $\unicode[STIX]{x1D648}$ are

with $M_{TA}=U-(1/2)\unicode[STIX]{x1D6FD}y^{2}$ . Analogous to the NTA example, the neutral stability direction $\boldsymbol{e}_{2}$ is parallel to the axis of the column, but this is now aligned with the traditional rotation component $\unicode[STIX]{x1D734}_{\bot }=(\unicode[STIX]{x1D6FD}y/2)\boldsymbol{k}$ , i.e. upwards/downwards. If $\unicode[STIX]{x1D706}_{1}<0$ for some latitudinal range, the instability direction is horizontally/laterally. If $A>1$ , two instability regions appear on either side of the equator where $\unicode[STIX]{x1D706}_{1}<0$ (see Kloosterziel et al. (Reference Kloosterziel, Orlandi and Carnevale2015) and figure 5 below). The outer edges of the instability region are again at ${\tilde{y}}={\tilde{y}}_{\pm }=\pm \sqrt{2\ln A}$ . Maximal growth is at ${\tilde{y}}=\pm {\tilde{y}}_{m}$ inside the two regions:

with $W_{0}$ the Lambert $W$ -function. For example, ${\tilde{y}}_{m}\approx 1.12$ when $A=5$ . The fastest non-dimensional rate of growth is

These simple considerations suggest that narrow flows (small $L/R_{0}$ ) will exhibit a far slower development of the instability when the traditional approximation (TA) is adopted than with the full Coriolis dynamics (NTA).

Figure 5. Contour plots of zonal vorticity $\unicode[STIX]{x1D714}_{x}=\unicode[STIX]{x2202}_{y}w-\unicode[STIX]{x2202}_{z}v$ from linearized simulations of the barotropic Gaussian jet in the TA with $A=5$ , $L/R_{0}=0.25$ and (a)  $Re=10^{5}$ , (b)  $Re=10^{6}$ . The initial condition is a random pointwise perturbation. The view is from positive to negative $x$ . The grey/black contours correspond to positive/negative  $\unicode[STIX]{x1D714}_{x}$ . The boundaries of the linear instability region are ${\tilde{y}}_{\pm }=\sqrt{2\ln A}\approx \pm 1.79$ . The latitudes of maximum instability are $\pm {\tilde{y}}_{m}\approx \pm 1.12$ . The equator ( ${\tilde{y}}=0$ ) is marked by the dotted vertical line. Only a portion of the computational domain is shown.

2.4 Viscous modifications of linear theory: scale selection

In the above simple theory, applied to the Gaussian jets in both the TA and NTA cases, no viscosity is taken into account. The subscript ‘ $\infty$ ’ in (2.27) and (2.31) indicates that this is the expected maximal rate of growth in the limit of vanishing viscosity ( $\unicode[STIX]{x1D708}=0$ ) or Reynolds number $Re\rightarrow \infty$ . In the case of the traditional approximation dynamics (TA) with $U=U(y)$ , and assuming $N(z)=N$ , a constant, a simple normal-mode analysis is possible with perturbations proportional to $\exp (\text{i}k_{z}z)$ . For inviscid flow, it is found that the overturning motions associated with inertial instability are amplified more rapidly the higher the value of $k_{z}$ ; that is, the smaller their vertical scale (e.g. Dunkerton Reference Dunkerton1981; Smyth & McWilliams Reference Smyth and McWilliams1998; Griffiths Reference Griffiths2003a ; Kloosterziel & Carnevale Reference Kloosterziel and Carnevale2008). This can be understood with the energy analysis of Ooyama (Reference Ooyama1966). Maximal energy is extracted from the background flow by motions predominantly in the unstable direction predicted by Solberg’s argument (see the end of § 2.1). Constant stratification $N$ introduces a lower-wavenumber cutoff for instability: only perturbations with $k_{z}>k_{c}^{N}$ can exhibit inertial instability. Increasing $N$ has the effect of increasing $k_{c}^{N}$ , thus requiring smaller vertical scales for the instability to proceed. This was demonstrated with numerical simulations by Kloosterziel et al. (Reference Kloosterziel, Carnevale and Orlandi2007a ).

The addition of viscosity, which preferentially damps small-scale motions, introduces a viscous upper wavenumber cutoff $k_{c}^{\unicode[STIX]{x1D708}}$ , which is such that $k_{c}^{\unicode[STIX]{x1D708}}\rightarrow \infty$ for $Re\rightarrow \infty$ . Thus, the instability can only grow if $k_{c}^{N}<k_{z}<k_{c}^{\unicode[STIX]{x1D708}}$ . That is, for finite but large enough Reynolds numbers, only perturbations with vertical scales within a finite range will be amplified. Within this range, a maximum growth rate is found at a specific vertical scale. Hence, if a flow is subjected to small perturbations and the ‘fastest’ mode is excited, one expects meridional motions with the vertical scale of the fastest growing mode to emerge. For the TA, the characteristic behaviour is that as $Re$ increases, the vertical scale of the fastest growing mode diminishes, and interestingly its horizontal scale also diminishes, with the mode becoming ever more concentrated in $y$ at the locations predicted by (2.30).

This scale selection in inertial instability for near-equatorial flows with the full Coriolis dynamics has to the best of our knowledge never before been investigated in any detail. For the NTA, we shall show in § 4 that, for the surface-confined jet (2.25) and $N=0$ , the upper limit $\unicode[STIX]{x1D6FE}_{\infty }$ (2.27) is approached via surface-confined modes that with increasing $Re$ have ever smaller scales and contract towards the equator (see figure 6 below). Asymptotic analysis detailed in appendix A explains much of the measured behaviour that we present in § 4.

Figure 6. Contour plots of $\unicode[STIX]{x1D714}_{x}$ from linearized simulations of the Gaussian jet in the NTA with $A=5$ , $L/R_{0}=0.25$ and (a)  $Re=10^{5}$ , (b)  $Re=10^{6}$ . The coordinate system of lines parallel and perpendicular to the planet rotation axis are drawn as thin dashed lines in the background. The aspect ratio of the plots is such that the units of length in the unscaled physical variables $y$ and $z$ are the same in each direction, ensuring that the aspect ratios of the rolls as shown are the physical aspect ratios. Only a portion of the computational domain is shown.

2.5 Nonlinear equilibration

In our first detailed study of the unfolding of inertial instability of barotropic vortices on the midlatitude $f$ -plane (Kloosterziel et al. Reference Kloosterziel, Carnevale and Orlandi2007a ), we discovered that for large $Re$ , the outcome of the fully nonlinear turbulent evolution can be predicted. This was extended to parallel shear flows on the $f$ -plane (Kloosterziel, Orlandi & Carnevale Reference Kloosterziel, Orlandi and Carnevale2007b ) and recently to barotropic flows like the Gaussian jet on the traditional equatorial $\unicode[STIX]{x1D6FD}$ -plane (Kloosterziel et al. Reference Kloosterziel, Orlandi and Carnevale2015). The prediction is that the instability mixes angular momentum or absolute momentum to such an extent that a new barotropic (depth-independent) inertially stable flow emerges. Stability is accomplished by mixing momentum to a homogeneous state with no gradients. Over some spatial range, this means that the flow has become neutrally stable (both eigenvalues of $\unicode[STIX]{x1D648}$ become zero). The range over which the mixing alters the original flow is determined by the constraint that total momentum is conserved.

For the full Coriolis dynamics, we predict that the instability will mix the absolute momentum from the surface downwards. The basic idea is sketched in figure 4. A thoroughly mixed region is expected to form between distance $R=R_{l}$ from the rotation axis and the surface. In our Cartesian formulation, this will be some parabola $\tilde{z}-{\tilde{y}}^{2}/2=\text{const}$ . Between this parabola and the surface, the new flow will have some constant momentum $m=M_{c}$ . The value of $M_{c}$ is determined by the ‘boundary condition’ that at the hypothesized location $R=R_{l}$ , the new momentum distribution is continuous, i.e. $M_{c}=M(R_{l})$ , and that the new flow is stable, which means that the mixed region (light grey area in figure 4 b) encompasses the original instability region that started at $R=R_{ins}$ , as sketched in figure 4(a) (dark grey area). With this information, we can determine both $R_{l}$ and $M_{c}$ by the constraint that the domain integrated momentum is conserved (it should be noted that, even at finite $Re$ , free-slip boundaries mean that the bulk momentum cannot change). Thus, one needs to solve

At the equator, the parabola describing $R_{l}$ is at non-dimensional depth $\tilde{z}_{l}$ . For the Gaussian jet with $A>1$ , we find with (2.32) how $\tilde{z}_{l}$ varies with $A$ by solving a transcendental relation between $\tilde{z}$ and $A$ (see appendix B, where details are provided). We verify these predictions in § 5.

4.1 Spatial structure of the fastest growing modes

For the TA, both the linear and nonlinear evolutions have been discussed in great detail in Kloosterziel et al. (Reference Kloosterziel, Orlandi and Carnevale2015). Here, we will present only one representative example of linearized evolution in the TA case for comparison with the NTA results. Figure 5 shows the fastest growing mode in the TA from the linearized numerical simulations. They are visualized by plotting the zonal vorticity $\unicode[STIX]{x1D714}_{x}=\unicode[STIX]{x2202}_{y}w-\unicode[STIX]{x2202}_{z}v$ . It should be noted that, throughout this paper, the aspect ratio of the contour plots is chosen so that the units of length in the physical, i.e. unscaled, coordinates $y$ and $z$ are equal on each axis. Thus, the aspect ratios of the rolls as seen in the plots are the physical aspect ratios. Therefore, we see in figure 5 that in the TA, the rolls have their major axes is aligned parallel to the surface. These vortices line up in two vertical stacks, one on either side of the equator. The stacks are centred at the locations predicted by (2.30), i.e. for $A=5$ at $\pm {\tilde{y}}_{m}\approx \pm 1.12$ . The bounds on the unstable region are at ${\tilde{y}}_{\pm }=\pm \sqrt{2\ln A}\approx \pm 1.79$ . By comparing the rolls in the two panels, we see that as $Re$ increases, both the width and the height of the rolls decrease. This shrinking of the rolls in both directions as $Re$ increases is not just for the equatorial problem, but is generally the case in inertial instability.

The major axes of the rolls are parallel to the horizontal instability direction $\boldsymbol{e}_{1}$ (2.29) associated with maximal growth of displacements in Solberg’s argument or the direction of maximal energy release in Rayleigh’s exchange argument (see the end of § 2.1). Moreover, they emerge spontaneously in the instability region with the shape envisioned by Ooyama (Reference Ooyama1966) for growth of meridional kinetic energy: highly elongated along the $\boldsymbol{e}_{1}$ direction.

The evolution in the NTA case is different in several important respects. First of all, it should be noted that the fastest growing mode shown for the wide-jet case ( $L/R_{0}=0.25$ ) in figure 6 for NTA simulations nicely verifies the pattern predicted in figure 2(b). Each roll in the mode is oriented with its long axis perpendicular to the parabolas of constant distance from the planet rotation axis.

Again, this is the $\boldsymbol{e}_{1}$ instability direction (2.22) via Solberg’s/Rayleigh’s argument or Ooyama’s proof of instability but with the full Coriolis force taken into account (see § 2.2). The resulting ‘stack’ of rolls is aligned with the Earth’s rotation axis, as predicted in figure 2. Another important difference from the TA case is that the fastest growing mode in the NTA is confined near the surface (see figure 6). This is due to the flow itself being surface-confined. This mode contracts towards a point located at the surface at the equator, i.e. towards $\{y,z\}=\{0,0\}$ as $Re$ increases. This is the predicted position of maximal growth according to the simple theory in § 2.2.

Figure 7. (a) Amplitude of the $\unicode[STIX]{x1D714}_{x}$ rolls seen in figure 5, as a function of latitude (solid lines). The latitudinal profile is Gaussian: $\exp (-({\tilde{y}}-{\tilde{y}}_{m})^{2}/2h^{2})$ , with ${\tilde{y}}=y_{m}\approx 1.12$ (dashed lines). The best fit for the scale $h$ diminishes with increasing $Re$ : $h=0.255$ , 0.164 and 0.125 respectively. (b) Amplitude of the most intense rolls near the equator like those shown in figure 6, as a function of depth (solid lines). The measured profiles are fitted with $\text{Ai}(-\unicode[STIX]{x1D706}_{0}-\tilde{z}/\unicode[STIX]{x1D701})$ (dashed lines, see text), where Ai is the Airy function, $\unicode[STIX]{x1D706}_{0}\approx 2.338$ and the vertical scale $\unicode[STIX]{x1D701}$ diminishes with increasing $Re$ : $\unicode[STIX]{x1D701}=0.326$ , $\unicode[STIX]{x1D701}=0.180$ and $\unicode[STIX]{x1D701}=0.100$ respectively.

In figure 5 for the TA, in each stack of rolls, the behaviour with depth $z$ is indistinguishable from purely sinusoidal, with scales decreasing with increasing $Re$ . The structure of the rolls with latitude (horizontally) is extremely close to Gaussians of the form $\exp (-({\tilde{y}}-{\tilde{y}}_{m})^{2}/2h^{2})$ , with ${\tilde{y}}_{m}$ the theoretical location of maximal growth and $h$ a horizontal scale that decreases with increasing $Re$ . This is illustrated with figure 7(a). In appendix A, we show that indeed each of the unstable modes in the TA is expected to have a horizontal structure proportional to the product of a Gaussian and a Hermite polynomial of some integral order $n$ , i.e. $H_{n}({\tilde{y}})\exp (-({\tilde{y}}-{\tilde{y}}_{m})^{2}/2h^{2})$ . The fastest growing mode corresponds to $n=0$ .

The cross-stream Gaussian structure of the modes centred about the line of maximal growth shown in figure 7(a) is canonical. Generally, for large $Re$ , in normal-mode analysis, an eigenvalue problem in the form of the quantum harmonic oscillator equation appears (see Dunkerton Reference Dunkerton1981; Bayly Reference Bayly1988; Sipp & Jacquin Reference Sipp and Jacquin2000; Griffiths Reference Griffiths2008a ,Reference Griffiths b and § A.1 in appendix A). It is always predicted for flows $U(y)$ for which  $\unicode[STIX]{x1D706}_{1}=fQ$ has isolated negative minima away from any horizontally confining boundaries. However, we believe that here, for the first time, such theoretical predictions have been confirmed through precise numerical simulations.

In the NTA, the rolls stack along the surface, as seen in figure 6. We find that with depth $z$ , the structure of the modes near the equator in the NTA is well matched by Airy function behaviour (see figure 7 b). They take the form $\text{Ai}(-\unicode[STIX]{x1D706}_{0}-\tilde{z}/\unicode[STIX]{x1D701})$ , with $\unicode[STIX]{x1D701}$ a vertical scale that diminishes with increasing $Re$ . The constant $\unicode[STIX]{x1D706}_{0}\approx 2.338$ is the lowest eigenvalue that sets $\text{Ai}(-\unicode[STIX]{x1D706}_{0})=0$ . This is the analogue of the fact that the fastest growing mode in the TA corresponds to the lowest eigenvalue ( $n=0$ ) of the Hermite functions. This behaviour too follows from the analysis in appendix A when assuming that all variations with latitude $y$ can be ignored. We believe that this is a truly novel result.

The measured amplitudes of the surface-confined patterns in figure 6 also have ever smaller latitudinal scales with increasing $Re$ . Quite remarkable is the observation that below the surface ( $z<0$ ) they cannot be distinguished from $\unicode[STIX]{x1D714}_{x}\propto \cos (l{\tilde{y}})\exp (-{\tilde{y}}^{2}/2\unicode[STIX]{x1D702}^{2})$ behaviour, with $\unicode[STIX]{x1D702}$ a number that decreases with increasing $Re$ , and $l$ a wavenumber that increases with $Re$ . Thus, there is oscillatory behaviour within a Gaussian envelope whose width, measured by $\unicode[STIX]{x1D702}$ , is not determined by the scale of the background jet which varies as $\exp (-y^{2}/2L^{2})=\exp (-{\tilde{y}}^{2}/2)$ .

For example, for $Re=10^{6}$ , we find $\unicode[STIX]{x1D702}=0.0374$ , so then the scale of the envelope is only approximately 4 % of that of the background flow. However, if we ignore this latitudinal variation of the envelope of the surface-confined modes for the NTA case, we can predict how rates of growth scale with the Reynolds number $Re$ and the horizontal/vertical scales. All details are given in appendix A, §§ A.3 and A.4.

4.2 Growth rates and horizontal/vertical scales as functions of $Re$ and $L$

In Kloosterziel et al. (Reference Kloosterziel, Orlandi and Carnevale2015), we studied the nonlinear development of the instability on the traditional $\unicode[STIX]{x1D6FD}$ -plane. Among other flows, the Gaussian barotropic jet was also considered. However, the effect of the variation of $Re$ on the width/height and the growth rate of the rolls has not been considered in any detail by us or anyone else, to the best of our knowledge. This is clearly of interest, and we shall now compare the results with the previously entirely unexplored NTA dynamics.

A series of linearized simulations of the type used to produce figure 5 (TA) and figure 6 (NTA) was performed over a wide range of values of $Re$ . The value of the growth rate of the fastest growing mode was calculated as half of the growth rate of the perturbation energy after all transients had died down and the growth appeared to be due to a single exponentially growing mode, with amplitude growing as $\exp (\unicode[STIX]{x1D6FE}t)$ . The resulting graph of $\unicode[STIX]{x1D6FE}$ versus $Re$ for fixed $A=5$ and $L/R_{0}=0.25$ is shown in figure 8(a). The TA data have been fitted with $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{\infty }-bRe^{-1/3}$ , with $\unicode[STIX]{x1D6FE}_{\infty }$ given by (2.31), which for $A=5$ and $L/R_{0}=0.25$ is $\unicode[STIX]{x1D6FE}_{\infty }=0.362$ . With the NTA, $\unicode[STIX]{x1D6FE}_{\infty }$ is determined by (2.27), which for $A=5$ means $\unicode[STIX]{x1D6FE}_{\infty }=\sqrt{5-1}=2$ . This is well fitted with $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{\infty }-bRe^{-1/4}$ in figure 8(a). Thus, in both the TA and the NTA, the growth rates asymptotically approach the appropriate inviscid upper bound as $Re\rightarrow \infty$ , but in the TA with $Re^{-1/3}$ while in the NTA with $Re^{-1/4}$ . This behaviour is predicted by the asymptotic analyses presented in appendix A.

Figure 8. Statistics for the fastest growing normal mode with $A=5$ . (a) Effects of variation of $Re$ on the growth rate in both the TA and the NTA. In both cases, $\unicode[STIX]{x1D6FE}_{\infty }$ is the theoretical inviscid upper limit (see text): in the TA, $\unicode[STIX]{x1D6FE}_{\infty }=0.362$ , and in the NTA, $\unicode[STIX]{x1D6FE}_{\infty }=\sqrt{1-A}=2$ . The parameter $b$ was estimated as the best least-squares fit to the data over the range $Re\in [10^{5}:10^{6}]$ : 1.04 (TA); 8.94 (NTA). These results are from simulations with $L/R_{0}=0.25$ . (b) Effects of variation of $Re$ on the spatial horizontal and vertical scales of the rolls in the TA. Here, $L_{v}$ and $L_{h}$ are given in units of $L^{2}/R_{0}$ and $L$ respectively, and a power-law fit to the $L_{v}$ and $L_{h}$ data is made with least-squares best-fit proportionality coefficients of $a=31.3$ and $a=1.45$ respectively. These results are from simulations with $L/R_{0}=0.25$ . (c) Effects of variation of $Re$ on the spatial horizontal and vertical scales of the rolls in the NTA. A power-law fit to the simulation data for $L_{v}$ and $L_{h}$ is made with least-squares best-fit proportionality coefficients of $a=4.74$ and $a=0.404$ respectively. These results are from simulations with $L/R_{0}=0.05$ . (d) Effects of variation in the NTA of $L/R_{0}$ on the number of rolls $n=2L/L_{h}$ within the range ${\tilde{y}}=[-1,+1]$ of the current. A least-squares fit over the range $L/R_{0}\in [0.01,0.25]$ to $n=a(L/R_{0})^{-b}$ yields $a=20.3$ and $b=0.584$ . The Reynolds number is $10^{5}$ .

In figure 8(b), we show how the height and width of the rolls vary with $Re$ for the TA dynamics. The measured vertical length scales $L_{v}$ of the typical rolls seen in figure 5 were determined via vertical Fourier transform of $\unicode[STIX]{x1D714}_{x}$ along the vertical centreline at ${\tilde{y}}_{m}$ . From the resulting power spectrum, the mean wavelength was determined and the width $L_{v}$ was computed as half of that wavelength. The measured horizontal scale $L_{h}$ of the rolls (with the Gaussian behaviour plotted in figure 7 a) was obtained as the horizontal half-width (distance between points where $\unicode[STIX]{x1D714}_{x}$ has its maximum value and half of that value) of the cell with the highest peak $|\unicode[STIX]{x1D714}_{x}|$ in the domain. The data in figure 8(b) have been fitted with $L_{v}=aRe^{-1/3}$ and $L_{h}=aRe^{-1/6}$ . These scaling laws were previously uncovered by Kloosterziel et al. (Reference Kloosterziel, Carnevale and Orlandi2007a ) for unstable circular vortices on the $f$ -plane. These exponents in the TA dynamics are expected from the general asymptotic analysis presented in appendix A, with the relevant material in § A.1, § A.2 adapted from Griffiths (Reference Griffiths2008b ) in Kloosterziel & Carnevale (Reference Kloosterziel and Carnevale2008).

For the NTA dynamics, the vertical scale $L_{v}$ was measured as the vertical half-width (with respect to $\unicode[STIX]{x1D714}_{x}$ ), measured vertically downward from the point of maximum $|\unicode[STIX]{x1D714}_{x}|$ in the domain. These scales determine the Airy function behaviour with depth plotted in figure 7(b). The result along with a power-law fit is shown in figure 8(c). To measure the horizontal scales $L_{h}$ in the NTA, we used a strategy similar to that used in the TA case, but here adapted to the fact that the rolls are primarily oriented vertically. The depth $z_{max}$ of the point at which the magnitude of the perturbation in $\unicode[STIX]{x1D714}_{x}$ was greatest was determined. The Fourier transform was taken on the horizontal line passing through this point, i.e. the transform of $\unicode[STIX]{x1D714}_{x}(y,z_{max})$ was performed in $y$ . From the resulting power spectrum, the mean wavelength was determined and the width $L_{h}$ was computed as half of that wavelength. The resulting data are shown as solid circles in figure 8(c) and have been fitted with $L_{v}=aRe^{-1/4}$ and $L_{h}=aRe^{-3/8}$ . It should be noted that these powers are different from what we showed in figure 8(b) for the traditional dynamics. These different exponents are again expected from the asymptotic analysis of appendix A in §§ A.3 and A.4, which, however, ignores all variation with latitude  $y$ .

Finally, we examine how the width of the rolls $L_{h}$ varies with the width of the current $L$ in the NTA. This is in preparation for the next section where we study the nonlinear evolution of the jet with the full Coriolis force. In order to emphasize the difficulty that we will encounter in resolving the fully nonlinear flow as $L/R_{0}$ decreases, in figure 8(d) we plot the number of rolls $n=2L/L_{h}$ that would fit within a distance $L$ from the equator (i.e. within ${\tilde{y}}=[-1,+1]$ ). A fit of the measured data to $n=a(L/R_{0})^{-b}$ gives roughly $a=20.3$ and $b=0.58$ . This suggest that $n$ tends to infinity as $L$ tends to zero. The increase in the number of rolls $n$ as $L/R_{0}$ decreases, or, in other words, the decrease of the latitudinal width of the rolls, makes it increasingly difficult to resolve the full current and still well resolve the inertial instability as $L/R_{0}$ decreases. This increase of the ‘density of rolls’ with decreasing curvature $L/R_{0}$ is as of yet unexplained.

5.1 Non-traditional approximation: vortex dynamics

In what follows, the code is run in the fully nonlinear mode. The initial condition is again as in the linear simulations: random initial perturbation to the basic state with $A=5$ . For the wide jet with $L/R_{0}=0.25$ , we showed the typical rolls near the surface in figure 6 for $Re=10^{5}$ and $Re=10^{6}$ that emerge in the linearized dynamics. In the fully nonlinear dynamics, early in the evolution, the rolls form and grow just as in the linear simulations.

As the rolls continue to amplify, nonlinear effects become evident. The primary change that occurs is that the vortices (the rolls seen in figure 6) start to ‘feel’ the shear created by their neighbours and begin to roll up in a manner already familiar from evolution on the $f$ -plane and the traditional $\unicode[STIX]{x1D6FD}$ -plane (see Kloosterziel et al. Reference Kloosterziel, Carnevale and Orlandi2007a ,Reference Kloosterziel, Orlandi and Carnevale b ; Plougonven & Zeitlin Reference Plougonven and Zeitlin2009; Carnevale et al. Reference Carnevale, Kloosterziel, Orlandi and van Sommeren2011, Reference Carnevale, Kloosterziel and Orlandi2013; Ribstein et al. Reference Ribstein, Plougonven and Zeitlin2014; Kloosterziel et al. Reference Kloosterziel, Orlandi and Carnevale2015). This action involves pairs of rolls of opposite signed vorticity. The tips of the pair roll up into a dipole. For each roll, one dipole tends to form at its tip near the surface and one at its tip deeper in the fluid. Dipoles are self-advecting structures. The dipoles that form near the surface are oriented so as to propagate towards the surface, which blocks their progress. The dipoles that form on the lower tips of the rolls are oriented so that they propagate downward towards greater depth. Through the agency of these propagating dipoles, the instability, which begins near the surface, fills out the region of linear instability and then goes beyond it.

This action can be seen in figure 9 for a case with $L/R_{0}=0.25$ and $Re=10^{5}$ . The red/blue contours correspond to positive/negative values of $\unicode[STIX]{x1D714}_{x}$ . The thick dashed parabola marks the limit of the linear instability region as determined by $\unicode[STIX]{x1D706}_{1}=0$ , which for the NTA via (2.26) gives the parabola $\tilde{z}-{\tilde{y}}^{2}/2=-\ln A$ (with ${\tilde{y}}=y/L$ , $\tilde{z}=zR_{0}/L^{2}$ ). The thick solid parabola marks the limit of the mixing range in the inviscid theory as determined by (B 8) in appendix B.

Figure 9. Contour plots of zonal vorticity $\unicode[STIX]{x1D714}_{x}$ from a fully nonlinear unforced simulation of the wide jet ( $L/R_{0}=0.25$ ) with $A=5$ and $Re=10^{5}$ in the NTA. The $|\unicode[STIX]{x1D714}_{x}|$ maximum in each plot is (a) 19, (b) 19, (c) 8.8, (d) 5.7. The corresponding contour increments $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{x}$ are (a) 4, (b) 2, (c) 1.5, (d) 0.5, and red/blue $=$ positive/negative. The thick black dashed parabola is the limit of the linear instability region: $\tilde{z}-{\tilde{y}}^{2}/2=-\ln A=-1.61$ for $A=5$ . The thick grey parabola is the limit of the predicted equilibrium mixing region (at ${\tilde{y}}=0$ at depth $\tilde{z}=\tilde{z}_{l}=-2.47$ , see text). The aspect ratio of the plots is such that the units of length in the unscaled physical variables $y$ and $z$ are the same in each direction, ensuring that the shape of the vortices as shown is the same as in the physical $y{-}z$ coordinates. Only a portion of the computational domain is shown. In hours/days (see text), the time sequence is (a) 9.5 h ( $t=5$ ), (b) 19 h ( $t=10$ ), (c) 2 days ( $t=25$ ) and (d) 7 days ( $t=90$ ).

In figure 9(a), we see that by non-dimensional time $t=5$ , which corresponds to $5/4\unicode[STIX]{x03C0}\approx 0.4$ days (9.5 h), the linear instability has grown sufficiently for the nonlinear effect of dipolar formation to have begun. In figure 9(b), we see that as the dipoles reach the limit of the theoretical mixing range, they begin to be diverted and some are distorted and destroyed. In figure 9(c), we see that some dipoles propagate somewhat beyond the limit of the mixing range before being stopped or turned around. The size of the dipoles decreases with increasing $Re$ , and so, as the results from our previous studies confirm, the higher the value of $Re$ is, the less these dipoles penetrate beyond the theoretical boundary of the mixing range (Kloosterziel et al. Reference Kloosterziel, Carnevale and Orlandi2007a ,Reference Kloosterziel, Orlandi and Carnevale b ; Carnevale et al. Reference Carnevale, Kloosterziel, Orlandi and van Sommeren2011, Reference Carnevale, Kloosterziel and Orlandi2013; Kloosterziel et al. Reference Kloosterziel, Orlandi and Carnevale2015). Being smaller, the dipoles at higher $Re$ can detect the mixing boundary more precisely, and they can have a smaller turning radius so that their motions can more effectively ‘describe’ where that boundary lies.

As the evolution proceeds, the absolute momentum in the mixing region becomes increasingly uniform (see below), the instability saturates and the perturbations begin to die away. By $t=90$ , the amplitudes of the vortices are down by a factor of approximately 4 from their peak values. We can note at this point the amazing ability of the small-scale structures, formed near the surface and around the equator, to effect the thorough mixing of momentum throughout the mixing region. Both the instability region and the mixing region are defined inviscidly. They do not change as $Re$ increases. On the other hand, the spatial scales of the most unstable mode shrink and the amplitude of the rolls becomes increasingly concentrated near $(y,z)=(0,0)$ as $Re$ increases. However, with increasing $Re$ , the dipoles, although of ever smaller scale, can persist longer and travel farther in the ever less viscous environment. Furthermore, their nonlinear interactions can create increasingly many more generations of vortices the higher the value of $Re$ is. Thus, even with the shrinking of the region of generation, and the shrinking of the scale of the initially generated dipoles, the entire mixing region can still be engulfed with vortices, come to be thoroughly mixed, and, as it turns out, mixed more uniformly and precisely, as $Re$ increases.

Another interesting feature in this simulation is that waves, with their crests roughly parallel to the parabolas drawn representing constant distance from the Earth’s rotation axis, are evident in the region outside the mixing region (see figure 9 d). These waves are created by disturbances in the mixing region and then propagate towards the bottom of the domain. They already exist at least by $t=10$ , but early on their amplitude is small compared with the maximum of $\unicode[STIX]{x1D714}_{x}$ in the vortices, so they are not evident in the early contour plots. However, as the peak amplitude of the vortices diminishes, the waves become evident in the contour plot if we decrease the contour level increment sufficiently, as we have done in figure 9(d). In an infinite domain, the waves could propagate away and completely escape, but in these simulations, in a finite box with free-slip conditions, they simply reflect from the bottom and from the walls and are only slowly dissipated by viscosity. Such waves have also been seen in studies of inertial instability on the traditional $f$ -plane (Kloosterziel et al. Reference Kloosterziel, Carnevale and Orlandi2007a ; Plougonven & Zeitlin Reference Plougonven and Zeitlin2009; Ribstein et al. Reference Ribstein, Plougonven and Zeitlin2014).

In most of the simulations that we have run, the waves are much less evident. This example was chosen to point out their presence. That these waves are simply inertial waves is implied by the facts that they propagate in the range where the main flow is very weak, their amplitude is small, suggesting linearity, there is no stratification and their frequencies of oscillation are in the expected range for internal gyroscopic waves. Since their nature is transitory with regard to our current examination of mixing of momentum  $m$ , we have not examined them further.

5.2 Absolute momentum mixing

Our theory for the nonlinear saturation of inertial instability is based on mixing of the absolute momentum, as discussed in § 2.5. In terms of the decomposition of the velocity that we are using here, the full absolute momentum is written as $m=(U+u)-(1/2)\unicode[STIX]{x1D6FD}y^{2}+2\unicode[STIX]{x1D6FA}z$ . Mixing is driven by the vortex dynamics illustrated in figure 9. By time $t=200$ in this case with $L/R_{0}=0.25$ , most residual vortex and wave activity is at a very low amplitude and the effect on $m$ can be clearly seen in figure 10. Figure 10(a) shows the initial distribution $m=M=U-(1/2)\unicode[STIX]{x1D6FD}y^{2}+2\unicode[STIX]{x1D6FA}z$ . It should be noted that the thick black dashed contour marks the line of maximal $M$ , that is where $\unicode[STIX]{x2202}M/\unicode[STIX]{x2202}z=0$ . This dashed line forms the boundary of the instability region. It should be noted that $\unicode[STIX]{x2202}M/\unicode[STIX]{x2202}z<0$ in the whole region extending from this dashed curve to the surface and that this negative gradient becomes very steep near the surface. In (b), we see that this gradient has been replaced by a nearly constant field. To be accurate, the region from the mixing boundary, the thick grey contour, up to the surface now exhibits only a weak positive vertical gradient of  $m$ . The contours outside the mixing region are essentially unchanged.

Figure 10. Contour plots of initial and late-time absolute momentum $\tilde{m}=m/\unicode[STIX]{x1D6FD}L^{2}$ from a simulation with $Re=10^{5}$ ,  $A=5$ and $L/R_{0}=0.25$ . All contours represent negative values. The contour increment is 0.16. The thick grey parabola is the predicted boundary of the uniformly mixed region expected in the inviscid limit, which at ${\tilde{y}}=0$ is at the depth $\tilde{z}=\tilde{z}_{l}=-2.47$ . The dashed black parabola is the line of maximum  $M$ , which forms the boundary of the instability region. At its lowest point (at ${\tilde{y}}=0$ ), this is at $\tilde{z}=\tilde{z}_{ins}=-1.61$ . The time $t=200$ corresponds to approximately 16 days.

Figure 11. Plots of the vertical profiles of (a) $\tilde{m}=m/\unicode[STIX]{x1D6FD}L^{2}$ and (b)  $\tilde{u} =u/\unicode[STIX]{x1D6FD}L^{2}$ at the equator at $t=200$ for $A=5$ , $L/R_{0}=0.25$ and three values of $Re$ , $10^{4}$ , $3\times 10^{4}$ and $10^{5}$ . Also shown are the initial conditions (thin solid black curve) and the predicted inviscid equilibrium (thick solid grey) for which $\tilde{z}_{l}\approx -2.47$ , $\tilde{M}_{c}\approx -2.89$ (see text). The results for $Re=10^{4}$ , $3\times 10^{4}$ and $10^{5}$ were averaged over 10, 100 and 200 realizations respectively. The time $t=200$ corresponds to 16 days.

Our inviscid prediction has two regions, that below the mixing boundary satisfying $m=M$ , the initial background distribution, and that above satisfying $m=M_{c}$ , a constant. However, for any finite $Re$ , the smoothing effect of viscosity tends to erode the separation between these two regions and to create a positive gradient throughout the mixing region. To get a better idea of this viscous effect, it is useful to examine the vertical profiles of $m$ at the equator for increasing $Re$ . In figure 11(a), we see that as $Re$ increases, the profile of $m$ at $t=200$ systematically approaches that of the prediction shown as the thick grey curve. The flat part of the curve is at the predicted level $\tilde{M}_{c}=M_{c}/\unicode[STIX]{x1D6FD}L^{2}\approx -2.89$ (see § 2.5 and appendix B). In order to eliminate effects of residual vortices or waves from the profiles, we have averaged the profiles over a number of realizations with different initial random seeds, as indicated in the figure caption.

Figure 11(b) shows the ensemble averaged $u$ profiles, indicating how absolute momentum mixing changes the mean flow. It should be noted that the inviscid theory predicts that mixing will cause a drop in speed in the centre of the jet (at $y=0$ , $z=0$ ) from the initial $\tilde{U} _{0}\equiv |U_{0}|/\unicode[STIX]{x1D6FD}L^{2}=-A=-5$ to $\tilde{u} =\tilde{M}_{c}\approx -2.89$ , i.e. a drop of approximately 40 % in magnitude.

We wish now to illustrate how the inertial instability proceeds for a much narrower jet with $L/R_{0}=0.05$ . Choosing this value for our background current, we find that we cannot simulate the full width of the jet and fully resolve the turbulent flow that results at $Re=10^{5}$ as we did in the case of $L/R_{0}=0.25$ treated above. The difficulty stems from the rapid increase in the number of rolls stacked in the latitudinal direction in the fastest growing mode as $L/R_{0}$ decreases, as illustrated in figure 8(d). To get around this, we use the following strategy: we keep the Reynolds number at the same value $Re=10^{5}$ and set $L/R_{0}$ to the desired value 0.05. However, we truncate the latitudinal span of the domain that we will simulate. We do this by placing free-slip vertical walls within the background current at ${\tilde{y}}=y/L=\pm 0.4$ . We resolve the latitudinal span with 2049 grid points. The vertical span of the computational domain is from the surface $\tilde{z}=0$ down to $\tilde{z}=zR_{0}/L^{2}=-8$ , which is covered by 1024 grid points.

Figure 12. Contour plots of $\unicode[STIX]{x1D714}_{x}$ from a fully nonlinear unforced simulation with a narrow jet $A=5$ , $L/R_{0}=0.05$ and $Re=10^{5}$ . The thick black dashed line is the limit of the linear instability region, as in the caption of figure 9. The thick grey solid line is the limit of the predicted equilibrium mixing region, which at ${\tilde{y}}=0$ is at the depth $\tilde{z}=\tilde{z}_{l}=-2.76$ (see the text). That prediction assumes that the unperturbed $M$ varies little in the narrow latitudinal band afforded to the flow by the confining free-slip walls at latitudes ${\tilde{y}}=\pm 0.4$ . The contour increment $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{x}=2$ for all panels (red/blue $=$ positive/negative). Only a portion of the computational domain is shown. The times correspond to (a) 9.5 h ( $t=5$ ), (b) 13 h ( $t=7$ ), (c) 19 h ( $t=10$ ) and (d) 34 h ( $t=18$ ).

Figure 13. Plots of vertical profiles of (a)  $\tilde{m}$ and (b)  $\tilde{u}$ at the equator at $t=25$ for three values of $Re$ , $3\times 10^{4}$ , $5\times 10^{4}$ and $10^{5}$ . Also shown are the initial conditions (thin solid black curve) and the predicted inviscid equilibrium (thick solid grey). The prediction is made assuming that the unperturbed $M$ varies little in the narrow latitudinal band afforded to the flow by the confining free-slip walls at latitudes ${\tilde{y}}=\pm 0.4$ . The predicted values are $\tilde{z}_{l}\approx -2.76$ , $\tilde{M}_{c}\approx -3.08$ (see text). The results for $3\times 10^{4}$ and $5\times 10^{4}$ are from ensembles of 20 realizations. The $Re=10^{5}$ result is averaged over 100 realizations ( $A=5$ , $L/R_{0}=0.05$ ). The time $t=25$ corresponds to 48 h.

The evolution of the vorticity field from a simulation with these parameters is shown in figure 12. The stages are very similar to those already seen in the $L/R_{0}=0.25$ case. There are no major differences save perhaps for the absence of evidence for waves in the latest stage of the flow, but even in the $L/R_{0}=0.25$ case, the strength of those waves varied greatly depending on the particular choice of seed for the initial random perturbation, the time of the contour plot and the contour level chosen.

In the late stages of the evolution, $m$ and $u$ approach the inviscid predictions as $Re$ increases in the equilibration of the narrow jet ( $L/R_{0}=0.05$ ). However, with the free-slip walls restricting the flow to a narrow latitudinal strip around the equator, the mixing prediction has to be altered. With no such latitudinal boundaries, the turbulent phase of the flow is free to mix momentum horizontally over a wide latitudinal range, as well as vertically. In the restricted flow, mixing must be almost entirely vertical. In this strip, the parabolas are nearly flat and there is almost no latitudinal variation in the initial absolute momentum $M(y,z)$ . Approximating $M(y,z)$ as $M(0,z)$ in this strip leads to a simplified relation predicting the new mixing depth and the new mixing value $M_{c}$ (see (B 9) in appendix B). For $A=5$ , this gives $\tilde{z}_{l}=-2.76$ . This is little deeper than for the unrestricted mixing which predicts $\tilde{z}_{l}=-2.47$ . The adjusted relation implies $\tilde{M}_{c}\approx -3.08$ , which is somewhat more negative than that for the unrestricted mixing where $\tilde{M}_{c}\approx -2.89$ .

The variation of $\tilde{m}$ with $Re$ is shown in figure 13(a) for $t=25$ . As in figure 11, these are averages over a number of realizations, as indicated in the figure caption. One can note how the profiles approach the inviscid prediction (thick grey curve) as $Re$ increases. For later times, the effect of viscosity in smoothing the velocity field will result in a positive gradient of $\tilde{m}$ in the mixing region for any finite $Re$ . Figure 13(b) is the corresponding figure for the dependence of $\tilde{u}$ on $Re$ in ‘equilibrium’. It should be noted that the inviscid theory predicts that $\tilde{u}$ will change from the initial value at the centre of the jet $\tilde{u} =-5$ to ${\approx}-3$ , a drop of 40 % of its initial magnitude, just as in the $L/R_{0}=0.25$ case.

Comparing the new equilibrium of the narrow jet shown in figure 13 at $t=25$ (after 2 days) to that of the wide jet in figure 11 at $t=200$ (after 16 days), it is seen that in the former, the ensemble averaged profiles at high $Re$ are closer to the prediction than in the latter. This must not be taken as proof that generally narrow jets reach the equilibrium more quickly than wide jets. In the simulations, small initial perturbations may project to a lesser extent on fast growing modes than on slower ones. This may determine how long it takes to reach the new equilibrium after nonlinear effects have taken over. There is also some arbitrariness in deciding when the flow has equilibrated (there were always some small-amplitude residual vortical and wave motions). Either way, the evolution may be considered to be quite rapid.

The narrow-jet example ( $L/R_{0}=0.05$ ) reveals that at high $Re$ , horizontal scales associated with the instability become quite small. In figure 12(a), we observe of the order of 50 rolls near the surface (this is subjective). The horizontal range corresponds to $0.8\times L$ , and in terms of ‘density of rolls’ $n$ contained within a span $2L$ (see the end of § 4.2), this is approximately $n\approx 2.5\times 50=125$ . In the linearized dynamics, we measured $n=a/(L/R_{0})^{b}$ , with $a=20.3$ and $b=0.584$ for $Re=10^{5}$ (§ 4.2; see figure 8 d). With $L/R_{0}=0.05$ , this empirical formula gives $n\approx 117$ . The ratio $L/R_{0}=0.05$ corresponds to approximately $3^{\circ }$ .

However, at $y=\pm L$ , the Gaussian jet still has approximately 60 % of its peak amplitude at $y=0$ . Only at $y=\pm 2L$ does it drop to 13 %. Defining the width of the jet as $W=4L$ , with $L/R_{0}=0.05$ , $W\approx 12^{\circ }$ . This is fairly wide, and one may consider narrower jets, with widths of say $W=2^{\circ }$ , which is closer to the lateral extent of some of the subsurface equatorial jets in the Pacific Ocean (see Luyten & Swallow Reference Luyten and Swallow1976; Firing Reference Firing1987; Hua, Moore & Le Gentil Reference Hua, Moore and Le Gentil1997). This width corresponds to $L/R_{0}\approx 0.009$ . With this value, the density of rolls is expected to be of the order of $n\approx 320$ . In other words, the horizontal scale of the instability rolls in such a narrow jet is near the surface of the order of $(1/300)^{\circ }$ . On the Earth, with $R_{0}=6380$  km, this means horizontal scales of 300–400 m if $Re=10^{5}$ . These scales get even smaller with larger  $Re$ .

In as much as the Reynolds number is not known or hard to define for oceanic flows, this may appear meaningless. However, it shows that if, in an ocean numerical model, $Re$ , based on $U_{0}$ , $L$ and $\unicode[STIX]{x1D708}$ , is of order $10^{5}$ or larger, the instability in jets with widths of just a few degrees will most likely not be resolved unless the computational domain is constrained to a very narrow latitudinal strip and very high horizontal resolutions are used. Apart from the speculative nature concerning the Reynolds number, there is another caveat: what is the influence of stratification on the instability? We shall briefly discuss this in our final section below.