2 Preliminary simulations

Consider a baroclinically unstable zonal flow represented by the two-layer quasi-geostrophic model (e.g. Pedlosky Reference Pedlosky1987)

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}{\displaystyle \frac{\unicode[STIX]{x2202}Q_{1}}{\unicode[STIX]{x2202}t}}+J(\unicode[STIX]{x1D6F9}_{1},Q_{1})=\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D6F9}_{1},\quad Q_{1}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6F9}_{1}+{\displaystyle \frac{f^{2}}{g^{\prime }H_{1}}}(\unicode[STIX]{x1D6F9}_{2}-\unicode[STIX]{x1D6F9}_{1})+\unicode[STIX]{x1D6FD}y,\\ {\displaystyle \frac{\unicode[STIX]{x2202}Q_{2}}{\unicode[STIX]{x2202}t}}+J(\unicode[STIX]{x1D6F9}_{2},Q_{2})=\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D6F9}_{2}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6F9}_{2},\quad Q_{2}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6F9}_{2}+{\displaystyle \frac{f^{2}}{g^{\prime }H_{2}}}(\unicode[STIX]{x1D6F9}_{1}-\unicode[STIX]{x1D6F9}_{2})+\unicode[STIX]{x1D6FD}y,\end{array}\right\}\end{eqnarray}$$

where $(\unicode[STIX]{x1D6F9}_{1},\unicode[STIX]{x1D6F9}_{2})$ are the streamfunctions, $(Q_{1},Q_{2})$ are the potential vorticities and $(H_{1},H_{2})$ are the depths of the upper and lower layers respectively; $g^{\prime }=(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C})g$ is the reduced gravity, $\unicode[STIX]{x1D70C}$ is the density, $f$ is the Coriolis parameter, $\unicode[STIX]{x1D6FD}=(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}y)$ , $\unicode[STIX]{x1D708}$ is the lateral viscosity and $\unicode[STIX]{x1D6FE}$ is the bottom drag coefficient. The flow field $(\unicode[STIX]{x1D6F9}_{1},\unicode[STIX]{x1D6F9}_{2})$ is separated into the Phillips basic state $(\bar{\bar{\unicode[STIX]{x1D713}}}_{1},\bar{\bar{\unicode[STIX]{x1D713}}}_{2})$ , representing the laterally homogeneous zonal current, and the perturbation $(\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D713}_{2})$ . Without loss of generality, we consider a basic state in which the lower layer is motionless ( $\bar{\bar{\unicode[STIX]{x1D713}}}_{2}=0$ ) and express the governing equations (2.1) in terms of $(\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D713}_{2})$ .

To reduce the number of governing parameters, the system is non-dimensionalized using $R_{d\,1}$ , $|U|$ and $R_{d\,1}/|U|$ as the scales of length, velocity and time respectively; $U$ is the basic velocity of the upper layer and $R_{d\,1}=\sqrt{g^{\prime }H_{1}}/f$ is the radius of deformation of the upper layer. The resulting non-dimensional system takes the form

(2.2) $$\begin{eqnarray}\hspace{-10.00002pt}\left.\begin{array}{@{}l@{}}{\displaystyle \frac{\unicode[STIX]{x2202}q_{1}}{\unicode[STIX]{x2202}t}}+J(\unicode[STIX]{x1D713}_{1},q_{1})+(\unicode[STIX]{x1D6FD}_{nd}+s){\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x}}+s{\displaystyle \frac{\unicode[STIX]{x2202}q_{1}}{\unicode[STIX]{x2202}x}}=\unicode[STIX]{x1D708}_{nd}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{1},\quad q_{1}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{1}+(\unicode[STIX]{x1D713}_{2}-\unicode[STIX]{x1D713}_{1}),\\ {\displaystyle \frac{\unicode[STIX]{x2202}q_{2}}{\unicode[STIX]{x2202}t}}+J(\unicode[STIX]{x1D713}_{2},q_{2})+(\unicode[STIX]{x1D6FD}_{nd}-sr){\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}x}}=\unicode[STIX]{x1D708}_{nd}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{2}-\unicode[STIX]{x1D6FE}_{nd}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{2},\quad q_{2}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{2}+r(\unicode[STIX]{x1D713}_{1}-\unicode[STIX]{x1D713}_{2}),\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{nd}=(\unicode[STIX]{x1D6FD}R_{d\,1}^{2})/|U|$ , $\unicode[STIX]{x1D708}_{nd}=\unicode[STIX]{x1D708}/(R_{d\,1}|U|)$ , $\unicode[STIX]{x1D6FE}_{nd}=(\unicode[STIX]{x1D6FE}R_{d\,1})/|U|$ , $r=(H_{1}/H_{2})$ are the key non-dimensional parameters, $s=U/|U|$ is the sign of the background velocity and $(q_{1},q_{2})$ are the perturbation PV fields. Thus, in non-dimensional units, the Phillips basic state is $\bar{\bar{u}}_{1,2}=-(\unicode[STIX]{x2202}\bar{\bar{\unicode[STIX]{x1D713}}}_{1,2}/\unicode[STIX]{x2202}y)=(s,0)=(\pm 1,0)$ .

Figure 1. Typical instantaneous patterns of the upper-layer streamfunction at various evolutionary stages. (a,c,e) The flow evolution in the high-drag experiment ( $\unicode[STIX]{x1D6FE}_{nd}=0.5$ ). The corresponding low-drag simulation ( $\unicode[STIX]{x1D6FE}_{nd}=0.05$ ) is shown in (b,d,f). The other governing parameters are identical: $(\unicode[STIX]{x1D6FD}_{nd},r,\unicode[STIX]{x1D708}_{nd},s)=(0.1,1/3,0.05,1)$ and $(L_{x},L_{y})=(400,200)$ . The reduction of the bottom drag coefficient results in the spontaneous generation of large-scale patterns.

In the following numerical experiments, we shall assume doubly periodic boundary conditions for $(\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D713}_{2})$ in $x$ and $y$ , and integrate the governing system (2.2) using the dealiased pseudospectral model employed in Radko et al. (Reference Radko, Peixoto de Carvalho and Flanagan2014). Figure 1(a,c,e) presents a typical simulation in the baroclinically unstable regime. For this experiment, we have used

(2.3) $$\begin{eqnarray}(\unicode[STIX]{x1D6FD}_{nd},r,\unicode[STIX]{x1D6FE}_{nd},\unicode[STIX]{x1D708}_{nd},s)=(0.1,13,0.5,0.05,1).\end{eqnarray}$$

In dimensional units, these parameters describe, for instance, an eastward basic current with $U=0.05~\text{m}~\text{s}^{-1}$ , $\unicode[STIX]{x1D6FD}=10^{-11}~\text{m}^{-1}~\text{s}^{-1}$ , $H_{1}=1~\text{km}$ , $H_{2}=3~\text{km}$ , $\unicode[STIX]{x1D6FE}=10^{-6}~\text{s}^{-1}$ , $\unicode[STIX]{x1D708}=60~\text{m}^{2}~\text{s}^{-1}$ and $R_{d\,1}=25~\text{km}$ – scales that are representative of typical mid-ocean flows. The size of the computational domain in figure 1(a,c,e) is $(L_{x},L_{y})=(400,200)$ , which is equivalent to $10\,000~\text{km}\times 5000~\text{km}$ , and it is resolved by a uniform mesh with $N_{x}\times N_{y}=768\times 384$ elements. The simulation was initialized from rest by a small-amplitude random computer-generated $(\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D713}_{2})$ distribution. The active statistically steady eddying motion driven by baroclinic instability was established after a few growth rate periods. Figure 1(a) shows a typical instantaneous ( $t=147$ ) streamfunction field ( $\unicode[STIX]{x1D713}_{1}$ ) in the initial stage of linear growth. As expected from linear stability theory, the most rapidly growing perturbations take the form of meridional harmonics. Figure 1(c) ( $t=227$ ) presents the second evolutionary stage – development of secondary instabilities, which act to distort the primary modes, adversely affecting their growth. Finally, by $t=308$ , figure 1(e), the system enters the quasi-equilibrium stage, characterized by irregular transient mesoscale patterns. The simulation was extended further in time (to $t=3000$ ), but no additional qualitative changes were observed.

A substantially different pattern emerged in the simulation shown in figure 1(b,d,f). This experiment is identical to that in figure 1(a,c,e) in all respects except for the bottom drag coefficient, which was reduced by an order of magnitude to $\unicode[STIX]{x1D6FE}_{nd}=0.05$ . While the initial stages of both simulations are similar, the long-term evolution of the low-drag experiment (figure 1 f) is characterized by the emergence of large-scale zonally elongated patterns. They gradually intensify and ultimately evolve into quasi-steady and largely barotropic structures dominating the final streamfunction field. The dramatically different outcomes of the experiments in figures 1(e) and 1(f) suggest that the spontaneous generation of LEDPs is controlled by the competition between the destabilizing action of mesoscale eddies and the damping tendency of bottom drag.

Since this study is focused on cross-scale interactions, it becomes necessary at this point to introduce a specific and convenient operational definition of a ‘large-scale’ pattern. In the following analysis, we assume that LEDPs consist of Fourier harmonics with wavelengths exceeding the cutoff scale $L_{cr}$ . This critical scale is set to $L_{cr}=20$ and this convention will be used consistently throughout our study. For instance, to further isolate the role of multiscale interactions in LEDP generation, we present a ‘filtered’ simulation (figure 2) in which all Fourier components with meridional wavelengths exceeding $L_{cr}=20$ are reset to zero at each time step, except for the fundamental meridional harmonic

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{1,2}=\text{Re}\left[A_{1,2}\exp \left(\text{i}{\displaystyle \frac{2\unicode[STIX]{x03C0}}{L_{y}}}y\right)\right]. & & \displaystyle\end{eqnarray}$$

Figure 2. The filtered simulation in which all Fourier components with meridional wavelengths exceeding $L_{cr}=20$ are reset to zero at each time step, except for the fundamental meridional harmonic (2.4). The upper- and lower-layer streamfunctions are shown in (a,c,e) and (b,d,f) respectively. The fundamental harmonic grows in time, eventually dominating the flow field, which reflects the destabilizing nature of mesoscale variability.

Such a wide spectral gap precludes the possibility of LEDP generation through the sequential coalescence of eddies into larger and larger structures, associated with the upscale energy cascade in predominantly two-dimensional turbulence. Additionally, the bottom drag in this experiment is applied to mesoscale components but not to the large-scale mode (2.4). Neglecting the large-scale bottom drag ( $\unicode[STIX]{x1D6FE}_{ls}=0$ ) but retaining the substantial drag on mesoscale components ( $\unicode[STIX]{x1D6FE}_{ms}=0.5$ ) allows us to isolate and quantify the impact of eddies on the large-scale current. The filtered simulation (figure 2) is characterized by the monotonic growth of (2.4) – the only large-scale mode present. The mode (2.4) represents the zonal flow and therefore it is baroclinically stable (recall that the basic Phillips flow is also zonal). Thus, its amplification is entirely due to the action of eddies. Since the cascade is precluded, the conclusion that one might draw from this experiment is that spectrally non-local interactions are sufficiently vigorous and efficient to induce the spontaneous generation of LEDPs. The evolutionary time scales and the resulting flow patterns in figure 2 are similar to those realized in its unfiltered counterpart (figure 1 b,d,f), in which bottom drag is uniformly low at all scales ( $\unicode[STIX]{x1D6FE}_{nd}=0.05$ ).

Figure 3(a) presents the time record of the upper-layer amplitude ( $|A_{1}|$ ) of the large-scale harmonic (2.4) as a function of time for the experiment in figure 2. The amplitude is plotted in logarithmic coordinates, and therefore the straight segment of the curve in figure 2 at moderate values of $|A_{1}|$ is indicative of the exponential growth of the perturbation. This, in turn, suggests that the emergence of the LEDP in this simulation can be interpreted as a form of instability driven by interactions between the large-scale harmonic and the mesoscale eddy field. To be more specific, one might argue that the acceleration of the large-scale current (figure 2 a,c,e) is a manifestation of the negative eddy viscosity phenomenon (Starr Reference Starr1968). This suggestion is also supported by the observation that the growth rate of the large-scale perturbation is proportional to its wavenumber $(m=2\unicode[STIX]{x03C0}/L_{y})$ squared, as expected for the linear diffusion equation with negative diffusivity.

Figure 3. (a) The amplitude of the fundamental harmonic (2.4) is plotted on a logarithmic scale as a function of time (solid curve) for the filtered calculation in figure 2. The relatively straight segment between $t=1000$ and $t=4000$ corresponds to exponential amplification, which, in turn, suggests linear instability of the large-scale mode (2.4). The straight dashed line represents the amplification corresponding to the uniform negative viscosity of $K=-1.2$ . (b) The growth rate of the large-scale perturbation is plotted as a function of its wavenumber ( $m$ ) in logarithmic coordinates. The numerical results are indicated by the plus signs and the straight line corresponds to the anti-diffusive scaling  $\unicode[STIX]{x1D706}\propto m^{2}$ .

Figure 3(b) presents a series of simulations which are identical to that in figure 2 in all respects except for the meridional extent of the domain ( $L_{y}$ ) – the latter is systematically varied from $L_{y}=100$ to $L_{y}=1600$ . The growth rate ( $\unicode[STIX]{x1D706}$ ) of the fundamental mode (2.4) is then diagnosed from numerics and plotted as a function of the large-scale wavenumber $m$ in logarithmic coordinates. The data points (figure 3 b) align along the straight line with slope corresponding to the power law $\unicode[STIX]{x1D706}\propto m^{2}$ . The effective eddy viscosity associated with this amplification is $K=-1.2$ . It is interesting that if the spectral gap is not enforced in the simulation, while the direct influence of bottom drag on the large-scale mode (2.4) is still excluded ( $\unicode[STIX]{x1D6FE}_{ls}=0$ ), the effective viscosity reduces to $K=-0.9$ . The limited sensitivity ( ${\sim}25\,\%$ ) of eddy viscosity to the presence/absence of intermediate scales $(L_{cr}<L<L_{y})$ is suggestive. It offers yet another indication that, in the parameter regime explored here, the evolution of large-scale patterns is controlled by spectrally non-local processes.

While the influence of eddies on large-scale perturbations is destabilizing, leading to the amplification of large-scale patterns, bottom drag and beta-effect act in the opposite sense. Figure 4 combines a series of simulations performed with various values of $(s,\unicode[STIX]{x1D6FD}_{nd},\unicode[STIX]{x1D6FE}_{nd})$ . The bottom drag in these experiments is spatially uniform at all scales and spectral filtering is not applied. For each simulation, we quantify the relative magnitude of mesoscale eddies and LEDPs using the discriminator variable

(2.5) $$\begin{eqnarray}R=\ln \left({\displaystyle \frac{E_{LEDP}}{E_{MS}}}\right),\end{eqnarray}$$

where $E_{MS}$ and $E_{LEDP}$ represent the perturbation energy contained in the mesoscale and large-scale eddy fields respectively. The mean perturbation energy of the two-layer system is

(2.6) $$\begin{eqnarray}E=\left\langle {\displaystyle \frac{r}{2(1+r)}}\left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{1}\right|^{2}+{\displaystyle \frac{1}{2(1+r)}}\left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{2}\right|^{2}+{\displaystyle \frac{r}{2(1+r)}}(\unicode[STIX]{x1D713}_{2}-\unicode[STIX]{x1D713}_{1})^{2}\right\rangle ,\end{eqnarray}$$

where the angled brackets represent averaging in $(x,y)$ . The mesoscale (large-scale) components are represented by a superposition of Fourier harmonics with the wavelengths $L<L_{cr}=20$ ( $L>L_{cr}$ ). The discriminator $R$ , diagnosed at the end of each simulation ( $t=2000$ ), is plotted as a function of $s\unicode[STIX]{x1D6FD}_{nd}$ and $\unicode[STIX]{x1D6FE}_{nd}$ in figure 4. For each value of $(s,\unicode[STIX]{x1D6FD}_{nd})$ , the fraction of energy contained in large scales monotonically increases with decreasing $\unicode[STIX]{x1D6FE}_{nd}$ , and the system eventually transitions to the LEDP-dominated regime. For given $\unicode[STIX]{x1D6FE}_{nd}$ , LEDPs are more likely to emerge at low values of  $\unicode[STIX]{x1D6FD}_{nd}$ .

Figure 4. The discriminator variable ( $R$ ) defined in (2.5) is plotted as a function of $s\unicode[STIX]{x1D6FD}_{nd}$ and $\unicode[STIX]{x1D6FE}_{nd}$ . The positive values of $R$ (shown in red) represent the LEDP-dominated simulations and the negative $R$ (shown in blue) reflect the predominantly mesoscale dynamics. The logarithmic (linear) axis is used for $\unicode[STIX]{x1D6FE}_{nd}$ ( $s\unicode[STIX]{x1D6FD}_{nd}$ ). The heavy solid curve represents the theoretical estimate (§ 6) of the boundary between the LEDP-dominated and mesoscale-dominated regimes.

The results presented in this section are generally consistent with the conventional views on the dynamics of mesoscale turbulence (Rhines Reference Rhines and Goldberg1977; Larichev & Held Reference Larichev and Held1995; Thompson & Young Reference Thompson and Young2006, Reference Thompson and Young2007). For instance, the link between weak bottom drag and spontaneous generation of large-scale structures has been noted in several studies (Nadiga Reference Nadiga2006; Nadiga & Straub Reference Nadiga and Straub2010; Radko et al. Reference Radko, Peixoto de Carvalho and Flanagan2014). The upscale barotropic cascade is known to be suppressed by large bottom drag as well as by large beta (e.g. Held & Larichev Reference Held and Larichev1996; Thompson & Young Reference Thompson and Young2006, Reference Thompson and Young2007; Berloff et al. Reference Berloff, Karabasov, Farrar and Kamenkovich2011). These effects may be invoked to broadly rationalize salient features of the presented experiments: (i) the limited amplitude of large-scale components and (ii) the dominance of spectrally non-local interactions when the bottom drag is substantial. The following model attempts to explain the eddy/mean-flow interaction on a more quantitative and mechanistic level (§§ 3–5). We shall focus on the moderately high-drag regime, which is readily amenable to multiscale treatment, and attempt to explain the transition to LEDP-favourable conditions in certain regions of the $(s\unicode[STIX]{x1D6FD}_{nd},\unicode[STIX]{x1D6FE}_{nd})$ parameter space (§ 6).

3 Generation of LEDPs as a multiscale problem

In this section, the problem of spontaneous generation of LEDPs is treated using techniques of asymptotic multiscale analysis (Kevorkian & Cole Reference Kevorkian and Cole1996; Mei & Vernescu Reference Mei and Vernescu2010). Details are relegated to the Appendix, and here we summarize the key steps. The interaction between vastly dissimilar flow components is described using new temporal and spatial variables ( $T$ , $Y$ ) over which the background eddy pattern is modulated:

(3.1a,b ) $$\begin{eqnarray}Y=\unicode[STIX]{x1D700}y,\quad T=\unicode[STIX]{x1D700}^{2}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D700}\ll 1$ is a measure of the relative spatial scales of the background eddies and large-scale flow. We are concerned by the ability of small-scale eddies, represented by the streamfunction fields $\bar{\unicode[STIX]{x1D713}}_{1,2}(x,y,t)$ , to affect the slow evolution of a broad zonal current. The primary eddy field is assumed to be fully developed and statistically homogeneous on large scales. In terms of magnitude, it represents the dominant temporally variable component, and $\bar{\unicode[STIX]{x1D713}}_{1,2}\sim O(1)$ in the expansion. Thus, the strength of primary eddies is comparable to the Phillips basic flow $\bar{\bar{u}}_{1,2}=(\pm 1,0)$ that generates them, and therefore the eddy dynamics are fundamentally nonlinear.

We now examine how this large-amplitude small-scale eddy field interacts with the large-scale zonal perturbation. The latter is taken to be barotropic at the leading order, which is consistent with our preliminary simulations (figure 2) and can be rationalized theoretically (see the Appendix):

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{1,2}^{(0)}=\unicode[STIX]{x1D713}^{(0)}(Y,T).\end{eqnarray}$$

Since the zonal jet (3.2) varies on large scales, its velocity is weak, $u^{(0)}=-\unicode[STIX]{x1D700}\unicode[STIX]{x1D713}_{Y}^{(0)}$ , and the shear induced by the jet is even weaker, $sh^{(0)}=\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D713}_{YY}^{(0)}$ . Both displacement and shear induced by the large-scale jet perturb and modulate the background eddy field $\bar{\unicode[STIX]{x1D713}}_{1,2}$ . These dynamics are reflected by expressing the total streamfunction field ( $\unicode[STIX]{x1D713}_{1,2}$ ) as follows:

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{1,2} & = & \displaystyle \unicode[STIX]{x1D713}^{(0)}(Y,T)+\bar{\unicode[STIX]{x1D713}}_{1,2}(x,y,t)+\unicode[STIX]{x1D700}\unicode[STIX]{x1D713}_{Y}^{(0)}(Y,T)\cdot \tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)}(x,y,t)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D713}_{YY}^{(0)}(Y,T)\cdot \tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)}(x,y,t)+O(\unicode[STIX]{x1D700}^{3}).\end{eqnarray}$$

Equation (3.3) can be readily obtained (see the Appendix) using a formal expansion of $\unicode[STIX]{x1D713}_{1,2}$ in powers of $\unicode[STIX]{x1D700}$ . As is common to all modulational stability analyses, the perturbation components $(\unicode[STIX]{x1D713}_{1,2}^{(1)},\unicode[STIX]{x1D713}_{1,2}^{(2)})$ in (3.3) are expressed as products of modulating functions $(\unicode[STIX]{x1D713}_{Y}^{(0)},\unicode[STIX]{x1D713}_{YY}^{(0)})$ , which vary on large spatial/temporal scales, and the corresponding auxiliary functions $(\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)},\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)})$ , varying on small scales:

(3.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{1,2}^{(1)}=\unicode[STIX]{x1D713}_{Y}^{(0)}(Y,T)\cdot \tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)}(x,y,t),\quad \unicode[STIX]{x1D713}_{1,2}^{(2)}=\unicode[STIX]{x1D713}_{YY}^{(0)}(Y,T)\cdot \tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)}(x,y,t).\end{eqnarray}$$

The auxiliary functions $(\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)},\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)})$ are related to the background pattern $\bar{\unicode[STIX]{x1D713}}_{1,2}$ by a set of partial differential equations (A 6), (A 9) and (A 10) in small-scale variables $(x,y,t)$ . Determining $\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)}$ and $\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)}$ for given $\bar{\unicode[STIX]{x1D713}}_{1,2}$ – the background small-scale eddy field – is referred to as solving the first and second auxiliary problems respectively.

The large-scale modulation of the eddy field in (3.4) is essential for triggering feedbacks of eddies on large-scale jets. The interaction of the jet (represented by the first term on the right-hand side of (3.3)) with the background eddy field (the second term in (3.3)) generates the modulated perturbations (terms three and four). These perturbations $(\unicode[STIX]{x1D713}_{1,2}^{(1)},\unicode[STIX]{x1D713}_{1,2}^{(2)})$ , in turn, interact with the basic field $\bar{\unicode[STIX]{x1D713}}_{1,2}$ , which modifies the large-scale pattern $\unicode[STIX]{x1D713}^{(0)}$ . The system evolution resulting from the interplay between the large-scale current and modulated eddy perturbations is represented by the diffusion equation

(3.5) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}T}}u^{(0)}=K{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}Y^{2}}}u^{(0)},\end{eqnarray}$$

where $K$ is the effective eddy viscosity (see the Appendix). For the sake of simplicity, in the present formulation we have linearized all asymptotic equations with respect to $\unicode[STIX]{x1D713}^{(0)}$ terms. However, it should be emphasized that multiscale models can incorporate the nonlinearities in the evolutionary large-scale equations in a rather straightforward manner (e.g. Gama et al. Reference Gama, Vergassola and Frisch1994; Novikov & Papanicolau Reference Novikov and Papanicolau2001).

In view of (3.1), it is apparent that (3.5) retains its form even when it is written using the original (small-scale) temporal and spatial variables $(t,y)$ . The multiscale model makes it possible to express $K$ in terms of the primary eddy field $\bar{\unicode[STIX]{x1D713}}_{1,2}$ and the auxiliary functions $(\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)},\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)})$ :

(3.6) $$\begin{eqnarray}K=-{\displaystyle \frac{1}{1+r}}\left[\left\langle r{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{1}}{\unicode[STIX]{x2202}x}}\left(2{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{1}^{(2)}}{\unicode[STIX]{x2202}y}}+\tilde{\unicode[STIX]{x1D713}}_{1}^{(1)}\right)+{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{2}}{\unicode[STIX]{x2202}x}}\left(2{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{2}^{(2)}}{\unicode[STIX]{x2202}y}}+\tilde{\unicode[STIX]{x1D713}}_{2}^{(1)}\right)\right\rangle \right],\end{eqnarray}$$

where angled brackets represent averaging in $(x,y)$ and square brackets represent averaging in $t$ . The instantaneous counterpart of (3.6), in which averages are taken in space only, will be denoted by $K_{inst}(t)$ :

(3.7) $$\begin{eqnarray}K_{inst}=-{\displaystyle \frac{1}{1+r}}\left\langle r{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{1}}{\unicode[STIX]{x2202}x}}\left(2{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{1}^{(2)}}{\unicode[STIX]{x2202}y}}+\tilde{\unicode[STIX]{x1D713}}_{1}^{(1)}\right)+{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{2}}{\unicode[STIX]{x2202}x}}\left(2{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{2}^{(2)}}{\unicode[STIX]{x2202}y}}+\tilde{\unicode[STIX]{x1D713}}_{2}^{(1)}\right)\right\rangle .\end{eqnarray}$$

It is interesting and perhaps somewhat counterintuitive that while the displacement-induced and shear-induced effects appear at different orders in the streamfunction series (3.3), both processes are represented at the leading order in the expression for eddy viscosity (3.6). The sign of $K$ determines whether eddies exert positive or negative feedback on large-scale patterns. If $K<0$ , we can expect eddy-induced amplification of large-scale patterns, but $K>0$ would imply that eddies act to suppress them.

One of the complications that commonly arises in multiscale modelling is that the interaction of primary eddies with large-scale flows may not be the only mechanism producing the perturbation fields $(\unicode[STIX]{x1D713}_{1,2}^{(1)},\unicode[STIX]{x1D713}_{1,2}^{(2)})$ . An alternative source of these perturbations is small-scale instability. This possibility is particularly relevant for the baroclinically unstable flows considered in our study. The irregular and time-dependent mesoscale variability at high Reynolds numbers is inherently susceptible to various secondary instabilities. These small-scale instabilities, by themselves, may not systematically affect large-scale patterns. Nevertheless, their presence complicates the analysis of the solutions obtained, since the unstable modes amplify in time and can eventually mask the processes of interest – the generation of perturbations through the interaction of the primary eddy field $\bar{\unicode[STIX]{x1D713}}_{1,2}$ with the large-scale flow $\unicode[STIX]{x1D713}^{(0)}$ and the feedback of these perturbations on the large-scale current. The majority of classical multiscale models deal with this complication either by restricting analysis to sufficiently low Reynolds numbers required to suppress the secondary small-scale instabilities (Meshalkin & Sinai Reference Meshalkin and Sinai1961; Nepomnyashchy Reference Nepomnyashchy1976; Sivashinsky Reference Sivashinsky1985; Frisch, She & Sulem Reference Frisch, She and Sulem1987; Gama et al. Reference Gama, Vergassola and Frisch1994; Wirth, Gama & Frisch Reference Wirth, Gama and Frisch1995; among many others) or by a priori ignoring the fast time variation of small-scale components (e.g. Novikov & Papanicolau Reference Novikov and Papanicolau2001; Radko Reference Radko2014). Both approaches are legitimate, as a starting point, and could be easily applied to our problem. However, this study provides us with the opportunity to develop an alternative non-invasive multiscale method, resulting in a more realistic description of eddy/LEDP interactions, which we describe next.

Figure 5. Typical upper/lower streamfunction patterns $\bar{\unicode[STIX]{x1D713}}_{1,2}$ (a,b) along with the corresponding solutions of the first (c,d) and second (e,f) auxiliary problems.

4 The AE model

Modulational stability analysis constitutes an effective analytical technique, commonly used to elucidate the dynamics of cross-scale interactions in most physical sciences (e.g. Mei & Vernescu Reference Mei and Vernescu2010). Perhaps the biggest obstacle for even wider application of multiscale models is the high sensitivity of model predictions to the assumed eddy structure, which severely limits the practical value of an otherwise very promising technique. The approach advocated in this study may offer a simple and effective remedy for the problem. Instead of focusing on idealized backgrounds (e.g. harmonic, hexagonal or dipolar) we propose to consider more realistic and dynamically consistent patterns that could be readily obtained from simulations. In the following examples, the primary eddy pattern $\bar{\unicode[STIX]{x1D713}}_{1,2}(x,y,t)$ is obtained by integrating the governing equations (2.2) subject to doubly periodic boundary conditions in $(x,y)$ . These calculations are accompanied by the integrations of auxiliary problems (A 9) and (A 10), making it possible to obtain $(\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1)},\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(2)})$ and, ultimately, the eddy viscosity (3.6). Experimentation with various computational domains indicated that the results are not particularly sensitive to the size as long as it significantly exceeds the typical eddy scale. Multiscale models tend to be slightly more accurate when the meridional domain size is increased, and therefore the solutions presented here were obtained using relatively large areas $(L_{x},L_{y})=(80,320)$ .

Figure 5 shows typical fully developed streamfunction patterns of the primary eddy field and the corresponding auxiliary functions, obtained as follows. First, the governing equations (2.2) were integrated in time starting from a small-amplitude random perturbation for $\bar{\unicode[STIX]{x1D713}}_{1,2}$ . The model parameters for the experiment in figure 5 are the same as in figure 1(a,c,e), and are listed in (2.3). It should be recalled that these parameters correspond to stable large-scale perturbations (figure 1 a,c,e) if bottom drag is assumed to be uniform at all scales ( $\unicode[STIX]{x1D6FE}_{ms}=\unicode[STIX]{x1D6FE}_{ls}=\unicode[STIX]{x1D6FE}_{nd}$ ) and to the amplification of large-scale modes (figures 2 and 3) if direct effects of drag on large scales are ignored ( $\unicode[STIX]{x1D6FE}_{ms}=\unicode[STIX]{x1D6FE}_{nd},\unicode[STIX]{x1D6FE}_{ls}=0$ ). By $t_{00}=300$ , the system evolved to the statistically steady quasi-equilibrium state. The experiment was then extended by $\unicode[STIX]{x0394}t=5$ and, during this interval, it was accompanied by concurrent integration of the auxiliary problems (A 9) and (A 10) starting from $\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1,2)}(t_{00})=0$ . Figure 5(a,b) presents the primary eddy fields $\bar{\unicode[STIX]{x1D713}}_{1,2}$ at $t=t_{00}+\unicode[STIX]{x0394}t$ , and the corresponding solutions of the first and second auxiliary problems are shown in figures 5(c,d) and 5(e,f) respectively. During the period $t_{00}<t<t_{00}+\unicode[STIX]{x0394}t$ , the eddy viscosity (3.7) gradually decreased from zero to $K_{inst}=-0.69$ , indicating that the eddy transfer of momentum occurred in the up-gradient sense.

Extensions of the auxiliary integrations for much longer periods of time, however, revealed a serious complication: the unbounded growth of the auxiliary functions $(\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1,2)})$ , caused by the instability of the linear operator $A$ in (A 6). As discussed earlier (§ 3), the appearance of new small-scale instabilities in the irregular and time-dependent eddying flows is natural and expected. The auxiliary problems (A 9) and (A 10) are structurally similar to the system (2.2) that governs the evolution of background eddies, and therefore they are subject to similar (baroclinic) instabilities. The interaction of the auxiliary functions with the background eddy field – the interaction represented by the Jacobian terms in $A$ – tends to intensify rather than suppress these instabilities. Unless this issue is addressed, the indefinite growth of the auxiliary functions can limit the utility of the multiscale model.

The crux of the proposed approach is a simple observation that even when operator $A$ is unstable, the asymptotic solution (3.3) still remains valid and accurate for relatively long periods of time. If $\unicode[STIX]{x1D706}_{A}$ represents the largest growth rate of $A$ , then the expansion leading to (A 9) and (A 10) fails (in the worst-case scenario) only after a period of $t_{f}\sim (1/\unicode[STIX]{x1D706}_{A})\ln (\unicode[STIX]{x1D700}^{-1})\gg 1$ . The time of adjustment of the eddy field to the large-scale forcing, on the other hand, is an order-one quantity: $t_{adj}=O(1)$ – the time scale set by a typical eddy turnaround period. The asymptotic disparity of the scales of adjustment and of the model failure ( $t_{adj}\ll t_{f}$ ) justifies the use of the multiscale framework for calculating – and physically explaining the selection of – the effective eddy viscosity.

The technical problem that arises at this point is that the amplification of auxiliary functions in time makes it difficult to accurately evaluate the time-mean values of $K$ in (3.6) from a single experiment. This complication, however, is readily resolved by ensemble averaging of the results over a large number ( $N$ ) of realizations. For a given set of parameters, the auxiliary problems were repeatedly integrated over time intervals $[t_{0},t_{0}+\unicode[STIX]{x0394}T]$ greatly exceeding the eddy adjustment time scale ( $\unicode[STIX]{x0394}T\gg t_{adj}$ ), where $t_{0}=t_{00}+n\unicode[STIX]{x0394}T,n=1,\ldots ,N$ . The initial conditions for each integration of auxiliary problems were taken to be $\tilde{\unicode[STIX]{x1D713}}_{1,2}^{(1,2)}(t_{0})=0$ . Then, the $N$ time records of $K_{inst}(\unicode[STIX]{x1D70F})$ , where $\unicode[STIX]{x1D70F}=t-t_{0}$ , obtained at each integration were averaged over all realizations, which effectively isolated the dynamically significant signal.

Figure 6 presents the ensemble-averaged time record of eddy viscosity $K_{inst}(\unicode[STIX]{x1D70F})$ for $N=445\,000$ and $\unicode[STIX]{x0394}T=25$ . As expected, the evolution of $K_{inst}(\unicode[STIX]{x1D70F})$ is a two-stage process. First, it gradually decreases from zero to its preferred equilibrium value ( $1<\unicode[STIX]{x1D70F}<5$ ) and subsequently remains largely uniform ( $\unicode[STIX]{x1D70F}>5$ ). The mean value evaluated for the second stage is interpreted as the eddy viscosity suggested by the multiscale model. While the ensemble-averaged $K_{inst}(\unicode[STIX]{x1D70F})$ converges with increasing $\unicode[STIX]{x1D70F}$ to a well-defined equilibrium level, the individual realizations are characterized by irregular amplifying oscillations about this mean value.

Figure 6. The ensemble-averaged eddy viscosity as a function of time. After the initial period of adjustment ( $\unicode[STIX]{x1D70F}<5$ ), the instantaneous eddy viscosity $K_{inst}$ obtained with the multiscale model (indicated by the solid curve) becomes largely uniform and closely matches the corresponding numerical estimate, represented by the dashed horizontal line.

In order to assess the accuracy of the multiscale model, the effective eddy viscosity was also diagnosed numerically. For this, we integrated the governing equations (2.2) using the version of the spectral code (cf. § 2) in which the bottom drag operator was applied to all spectral modes except for the fundamental harmonic (2.4). As a result, this mode amplified exponentially due to the action of eddies. The growth rate $(\unicode[STIX]{x1D706})$ realized in this simulation was then used to infer the eddy viscosity $K_{num}=-\unicode[STIX]{x1D706}/m^{2}$ , where $m=2\unicode[STIX]{x03C0}/L_{y}$ is the fundamental wavenumber. The numerical estimate of eddy viscosity is indicated in figure 6 by the horizontal dashed line. The close agreement between the numerical and theoretical estimates leaves no doubt that the formalism of multiscale modelling can be successfully applied to realistic and dynamically consistent eddy patterns.

To explore the parameter space, we performed a series of simulations analogous to that in figure 6 for a range of $s\unicode[STIX]{x1D6FD}_{nd}$ and compared (figure 7) the multiscale results with the corresponding numerical estimates of eddy viscosity. Figure 7 indicates that the multiscale model not only correctly produces the negative sign and typical magnitude of eddy viscosity but also accurately captures the variation of $K$ with $s\unicode[STIX]{x1D6FD}_{nd}$ . The pattern of $K(s\unicode[STIX]{x1D6FD}_{nd})$ is of interest in its own right. While the intensity of eddies, as measured for instance by the root mean square perturbation velocity, tends to be significantly higher for the westward propagating basic current ( $s=-1$ ), the eddy viscosity pattern is reversed. The largest values of eddy viscosity $(|K|\sim 0.9)$ are found for $s=1$ , which suggests that baroclinic instability of eastward flows is much more effective in terms of the spontaneous generation of LEDPs.

Figure 7. The effective eddy viscosity ( $K$ ) diagnosed from numerical simulations (solid curve) is plotted as a function of $s\unicode[STIX]{x1D6FD}_{nd}$ along with the corresponding theoretical estimate based on the ensemble-averaged multiscale model (plus signs).

5 The anatomy of counter-gradient momentum transport

The ability of the multiscale model to reproduce eddy viscosity (figures 6 and 7) opens an attractive opportunity to (i) examine various processes involved in the selection of $K$ , (ii) identify the dominant balances realized at each order in the asymptotic expansion and (iii) physically interpret the chain of events leading to counter-gradient momentum transport.

The first question arising after inspection of (3.7) is the relative contributions of displacement-induced and shear-induced perturbations in the top/bottom layers to the net eddy viscosity $K_{inst}=K_{1}^{(1)}+K_{1}^{(2)}+K_{2}^{(1)}+K_{2}^{(2)}$ , where

(5.1) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}K_{1}^{(1)}=-{\displaystyle \frac{r}{1+r}}\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{1}}{\unicode[STIX]{x2202}x}}\tilde{\unicode[STIX]{x1D713}}_{1}^{(1)}\right\rangle ,\quad K_{1}^{(2)}=-{\displaystyle \frac{2r}{1+r}}\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{1}}{\unicode[STIX]{x2202}x}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{1}^{(2)}}{\unicode[STIX]{x2202}y}}\right\rangle ,\\ K_{2}^{(1)}=-{\displaystyle \frac{1}{1+r}}\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{2}}{\unicode[STIX]{x2202}x}}\tilde{\unicode[STIX]{x1D713}}_{2}^{(1)}\right\rangle ,\quad K_{2}^{(2)}=-{\displaystyle \frac{2}{1+r}}\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{2}}{\unicode[STIX]{x2202}x}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D713}}_{2}^{(2)}}{\unicode[STIX]{x2202}y}}\right\rangle .\end{array}\right\}\end{eqnarray}$$

These components were evaluated from the ensemble-averaged simulation in figure 6. The values of the eddy viscosity components (5.1) recorded during the period of eddy adjustment to the large-scale forcing ( $\unicode[STIX]{x1D70F}=2$ ) and during the quasi-equilibrium stage ( $\unicode[STIX]{x1D70F}=10$ ) are listed in table 1. In both regimes, the dominant component of (3.6), which controls the magnitude and sign of eddy viscosity, is the upper-layer displacement term $K_{1}^{(1)}<0$ . The lower-layer displacement term $K_{2}^{(1)}$ is also negative, but it is much weaker and therefore plays a secondary role in LEDP generation. The upper- and lower-layer shearing terms ( $K_{1,2}^{(2)}$ ) are both positive and act to reduce the counter-gradient flux driven by the upper-layer displacement mode ( $K_{1}^{(1)}$ ). It should be noted, however, that the magnitude of the amplifying terms ( $K_{1,2}^{(1)}$ ) is close to that of the stabilizing components ( $K_{1,2}^{(2)}$ ). Thus, even small variations in each of the four terms in (5.1) could change the sign of the net viscosity.

Table 1. Components of the ensemble-averaged eddy viscosity from the multiscale model at the adjustment ( $\unicode[STIX]{x1D70F}=2$ ) and the equilibrium ( $\unicode[STIX]{x1D70F}=10$ ) stages. The dominant contributor to eddy viscosity is the displacement-driven upper-layer component $K_{1}^{(1)}$ .

To explain the origin of negative viscosity, we focus on the first auxiliary problem (A 9), which governs the destabilizing displacement mode. First, we note that the exact analytical solution of (A 9) is given by

(5.2) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D713}}_{i}^{(1)}={\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{i}}{\unicode[STIX]{x2202}x}}(t-t_{0}),\quad i=1,2.\end{eqnarray}$$

The validity of (5.2) can be readily ascertained by differentiating the governing equations for the primary eddies in $x$ , multiplying the result by $(t-t_{0})$ and subtracting from the corresponding first auxiliary equations. The physical interpretation of (5.2) is also straightforward: at the leading order, the large-scale perturbation results in the kinematic shift of the eddy pattern. Despite linear instability of the first auxiliary problem, the analytical solution (5.2) was found to be close to the corresponding numerical realizations. The relative differences between the analytical and numerical solutions accumulating over the finite integration interval $\unicode[STIX]{x0394}T=$ 25 were of the order of 0.2 %.

The structure of the first auxiliary solution (5.2) yields insight into the mechanisms controlling the meridional potential vorticity flux associated with the displacement mode. The perturbation PV flux associated with the displacement mode in each layer contains two distinct components – advection of the basic PV by the perturbation velocity $\langle \bar{q}_{i}v_{i}^{(1)}\rangle$ and advection of the perturbation PV by the basic eddy velocity $\langle q_{i}^{(1)}\bar{v}_{i}\rangle$ , where $i=1,2$ . However, it is argued in the Appendix that $\langle \bar{q}_{i}v_{i}^{(1)}\rangle$ does not contribute to the net PV flux and the relevant component is $F_{qi}^{(1)}=\langle q_{i}^{(1)}\bar{v}_{i}\rangle$ . Furthermore, at the leading order, this flux reduces to

(5.3) $$\begin{eqnarray}F_{qi}^{(1)}=\unicode[STIX]{x1D700}^{3}\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D713}}_{i}}{\unicode[STIX]{x2202}x}}\tilde{\unicode[STIX]{x1D713}}_{i}^{(1)}\right\rangle {\displaystyle \frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}^{(0)}}{\unicode[STIX]{x2202}Y^{3}}}.\end{eqnarray}$$

In view of (5.2), PV fluxes in both layers are positively correlated with $(\unicode[STIX]{x2202}q^{(0)}/\unicode[STIX]{x2202}Y)$ , where $q^{(0)}=(\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}^{(0)}/\unicode[STIX]{x2202}Y^{2})$ is the vorticity of the large-scale flow. Therefore, we conclude that the sign of the depth-weighted vertical average of the PV fluxes in both layers, $F_{q}=(rF_{q1}^{(1)}+F_{q2}^{(1)})/(1+r)$ , is also controlled by the direction of the large-scale vorticity gradient:

(5.4) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}F_{q}>0\quad \text{for }{\displaystyle \frac{\unicode[STIX]{x2202}q^{(0)}}{\unicode[STIX]{x2202}y}}>0,\\ F_{q}<0\quad \text{for }{\displaystyle \frac{\unicode[STIX]{x2202}q^{(0)}}{\unicode[STIX]{x2202}y}}<0.\end{array}\right\}\end{eqnarray}$$

Thus, the barotropic vorticity flux is up-gradient, which results in the amplification of LEDPs. The physical mechanism of amplification is illustrated in figure 8, which depicts the large-scale zonal current ( $\bar{\bar{u}}+u^{(0)}$ ) with $y$ -dependent velocity and the corresponding vorticity pattern $q^{(0)}=-(\unicode[STIX]{x2202}u^{(0)}/\unicode[STIX]{x2202}y)$ . According to (3.1a,b ), the meridional eddy-induced PV flux is positive for $y_{1}<y<y_{2},y_{3}<y<y_{4}$ and negative for $y_{2}<y<y_{3}$ , where the $y_{i}$ denote the extrema of the large-scale vorticity. Thus, the divergence of the potential vorticity flux $((\unicode[STIX]{x2202}F_{q}/\unicode[STIX]{x2202}y)>0)$ results in the reduction of potential vorticity in the regions where the large-scale potential vorticity is already negative $(q^{(0)}<0)$ . In regions where the large-scale potential vorticity is positive, it will be further increased. This positive feedback mechanism explains the eddy-induced amplification of large-scale perturbations and the associated negative eddy viscosity.

Figure 8. Schematic diagram illustrating the acceleration of a large-scale zonal flow by eddies. Extrema of the large-scale vorticity are indicated by $y_{i}$ . The regions where the large-scale vorticity is positive are indicated by shading. The multiscale model – see (3.1a,b ) – suggests that the eddy PV flux ( $F_{q}$ ) increases with the PV gradient $(\unicode[STIX]{x2202}q_{0}/\unicode[STIX]{x2202}y)$ . Therefore, we expect that the divergence/convergence of fluxes will further increase the vorticity in the regions where it is already positive ( $q^{(0)}>0$ ) and reduce it where it is negative ( $q^{(0)}<0$ ).

6 Conditions for the spontaneous generation of LEDPs

The success of the multiscale model in reproducing numerical results raises the question of whether it can also be used to explain the spontaneous emergence of LEDPs in certain regions of the parameter space and the lack of them in others (e.g. figure 4). The formation of LEDPs in simulations is apparently controlled by the competition between the destabilizing action of mesoscale eddies and the stabilizing tendency of the bottom drag acting directly on large scales. The following analysis attempts to identify and rationalize the onset of the LEDP-favourable regime as the bottom drag is gradually decreased from relatively high values, sufficient to stabilize large-scale perturbations, towards marginally unstable conditions.

The multiscale analysis (§§ 4 and 5) suggests that the effects of mesoscale eddies on the barotropic zonal current $u^{(0)}(y)$ can be represented by negative eddy viscosity (3.6). The impact of the bottom drag can be assessed by vertically averaging the governing equations (2.2), which in the absence of eddies amounts to $(\unicode[STIX]{x2202}u^{(0)}/\unicode[STIX]{x2202}t)=-(\unicode[STIX]{x1D6FE}_{nd}/(1+r))u^{(0)}$ . Assuming that the net effect can be estimated by linearly adding the two tendencies, the evolutionary equation for the large-scale flow reduces to

(6.1) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}u^{(0)}}{\unicode[STIX]{x2202}t}}=K{\displaystyle \frac{\unicode[STIX]{x2202}^{2}u^{(0)}}{\unicode[STIX]{x2202}y^{2}}}-{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{nd}}{1+r}}u^{(0)}.\end{eqnarray}$$

The existing evidence suggests that the linearly most unstable LEDPs take the form of zonal jets. This property is clearly revealed by numerical simulations (e.g. Berloff et al. Reference Berloff, Kamenkovich and Pedlosky2009) and can be attributed to the generic tendency of the beta-effect to suppress large-scale meridional displacements. Therefore, it is reasonable to assume that the point of destabilization of large-scale zonal modes represents the transition between LEDP-dominated and mesoscale-dominated regimes.

The stability of large-scale zonal flows is analysed using the linear model (6.1). The substitution of normal modes $u^{(0)}=A\exp (\text{i}my+\unicode[STIX]{x1D706}t)$ in (6.1) reveals that for $K<0$ , the large-scale flow is stable at relatively long wavelengths ( $m<m_{0}$ ) and unstable at short ones ( $m>m_{0}$ ), where

(6.2) $$\begin{eqnarray}m_{0}=\sqrt{-{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{nd}}{K(1+r)}}}.\end{eqnarray}$$

If mesoscale eddies are contained within the range $m>m_{cr}$ , then the condition for spontaneous formation of LEDPs becomes

(6.3) $$\begin{eqnarray}m_{0}<m_{cr}.\end{eqnarray}$$

If (6.3) is satisfied, then there is a finite range of sufficiently long wavelengths $m_{0}<m<m_{cr}$ that are unstable and therefore are likely to amplify.

The criterion (6.3) was used to formulate an algorithm for identifying the regions in the parameter space $(s\unicode[STIX]{x1D6FD}_{nd},\unicode[STIX]{x1D6FE}_{nd})$ that are susceptible to spontaneous LEDP generation as follows. We have assumed (§ 2) that the separation between meso- and large scales occurs at a wavelength $L_{cr}=20$ , and therefore the transitional wavenumber is $m_{cr}=(2\unicode[STIX]{x03C0}/L_{cr})=0.314$ . For any given value of $s\unicode[STIX]{x1D6FD}_{nd}$ , we performed a series of calculations in which bottom drag (assumed here to be uniform at all scales) was systematically decreased. The experiments started from a sufficiently large initial value ( $\unicode[STIX]{x1D6FE}_{nd}=0.5$ ), which ensures that the eddy field is mesoscale-dominated. For each $\unicode[STIX]{x1D6FE}_{nd}$ , we evaluated the eddy viscosity $K$ using the multiscale model and computed $m_{0}$ from (6.2). As $\unicode[STIX]{x1D6FE}_{nd}$ was decreased, $m_{0}$ decreased as well. Eventually, $m_{0}$ reached the mesoscale/large-scale threshold ( $m_{cr}$ ), at which point the calculation was stopped. Further decrease in $\unicode[STIX]{x1D6FE}_{nd}$ would imply the existence of a finite range of long unstable wavelengths ( $m_{0}<m<m_{cr}$ ). Therefore, the critical value of the bottom drag coefficient ( $\unicode[STIX]{x1D6FE}_{cr}$ ) for which $m_{0}=m_{cr}$ was assumed to represent the point of transition between the mesoscale- and LEDP-dominated regimes. This procedure was repeated for a series of $s\unicode[STIX]{x1D6FD}_{nd}$ , and the relation $\unicode[STIX]{x1D6FE}_{cr}(s\unicode[STIX]{x1D6FD}_{nd})$ obtained in this manner is shown by the heavy curve in figure 4. The theoretical curve is in general agreement with the numerics, indicating that the LEDP-dominated region lies at relatively low values of $\unicode[STIX]{x1D6FD}_{nd}$ and $\unicode[STIX]{x1D6FE}_{nd}$ . The ability of the proposed model to reproduce conditions for the spontaneous emergence of LEDPs in simulations is suggestive. It lends credence to all theoretical elements on which the algorithm is based, most notably the idea of competition between stabilizing bottom drag and destabilizing negative eddy viscosity (6.1), multiscale analysis, and the criterion for LEDP formation in (6.3).