2. Coupled model of NIW–mean flow interaction on a $\beta$-plane

For simplicity and numerical efficiency and to capture the key feature of wave–mean flow interaction, we study the interaction between a single-vertical-mode NIW and a barotropic QG mean flow using the 2-D model derived by Xie & Vanneste (Reference Xie and Vanneste2015) on a $\beta$-plane:

(2.1a)$$\begin{gather} \partial_t q+{J}(\psi,q)=0 , \end{gather}$$

(2.1b)$$\begin{gather}\partial_t\phi+{J}(\psi,\phi)+\mathrm{i} g(y)\phi-\frac{\mathrm{i} f_0}{2}{L}^2\phi+\frac{\mathrm{i}}{2}\nabla^2\psi\phi=0, \end{gather}$$

where

(2.2)\begin{equation} q=\nabla^2\psi+ g(y) +\frac{\mathrm{i} f_0}{2}{J}(\phi^*,\phi)+\frac{f_0}{4}\nabla^2|\phi|^2\end{equation}

is the NIW-modified quasi-geostrophic potential vorticity (QGPV). Here, $\psi$ is the Lagrangian mean streamfunction and the horizontal velocity of QG mean flows is $(u,v)=(-\partial _y \psi, \partial _x \psi )$; $f_0=\text {const.}$ is the local Coriolis frequency; and $\phi$ is the complex wave amplitude and links to the two horizontal velocities $(u_0,v_0)$ through $u_0+\mathrm {i} v_0=-\mathrm {i} f_0 \phi \, \mathrm {e}^{-\mathrm {i} f_0 t+\mathrm {i} m z}$ with $m$ the vertical wavenumber of the single-vertical-mode NIW. The asterisk denotes the complex conjugate and ${J}(a,b)=a_xb_y-a_yb_x$ is the Jacobian. Also, $g(y)$ is the variation of the Coriolis parameter $f$ in the meridional direction such that $f=f_0+g(y)$; therefore, the $\beta$-effect is captured by $\beta =\partial _y g$. This model consists of the equation describing the material-invariant QGPV (2.1a) and the linear equation capturing the slow modulation of NIW amplitude (Young & Ben Jelloul Reference Young and Ben Jelloul1997). In (2.1b) we consider a modified dispersion (Asselin & Young Reference Asselin and Young2019) with the operator

(2.3)\begin{equation} {L}^2=\frac{\nabla^2}{\dfrac{f_0^2}{N^2}m^2-\dfrac{1}{4}\nabla^2}, \end{equation}

where $N$ is the buoyancy frequency. Here, the Laplacian operator in the denominator should be understood in the spectral space. This modified dispersion term captures the high-order corrections of the NIW dispersion relation and improves the numerical efficiency by bounding the frequency as the wavenumber increases.

For simplicity, we introduce a modified meridional-dependent $\beta$-effect, $\beta (y)$, which is periodic in $y$ instead of the traditional choice of a constant, to perform numerical simulations in a periodic domain. Practically, to extend the region with a constant $\beta$, we prescribe

(2.4) \begin{equation} \beta(y)= \begin{cases} \beta_0 \dfrac{\mathrm{e}^{l_b y}-\mathrm{e}^{{-}l_b y}}{\mathrm{e}^{l_b y}+\mathrm{e}^{{-}l_b y}}, & 0\leq y < L_y/4 ,\\ \beta_0 \dfrac{\mathrm{e}^{l_b (L_y/2-y)}-\mathrm{e}^{{-}l_b (L_y/2-y)}}{\mathrm{e}^{l_b (L_y/2-y)}+\mathrm{e}^{{-}l_b (L_y/2-y)}}, & L_y/4\leq y < 3L_y/4 ,\\ \beta_0 \dfrac{\mathrm{e}^{l_b ({-}L_y+y)}-\mathrm{e}^{{-}l_b ({-}L_y+y)}}{\mathrm{e}^{l_b ({-}L_y+y)}+\mathrm{e}^{{-}l_b ({-}L_y+y)}}, & 3L_y/4\leq y < L_y , \end{cases} \end{equation}

where $\beta _0=\text {const.}$, $l_b$ is a constant that controls the range of the variation in $\beta$ and $L_y$ is the domain size in the meridional direction. Our analysis uses the data obtained from the region $\approx [1/l_b, L_y/2-1/l_b]$, where $\beta$ is almost a constant. The other half of the domain is only introduced for numerical simulation, and it is unphysical because $\beta$ should be positive in the northern hemisphere. We use the current periodic set-up mainly because of its simplicity and computational efficiency. The expressions of $\beta (y)$ and $g(y)$ are illustrated in figure 1. An alternative set-up is to consider a constant-$\beta$ system in a channel, which is widely used. For example, Abernathey, Marshall & Ferreira (Reference Abernathey, Marshall and Ferreira2011) numerically simulated a $\beta$-plane channel with a sponge layer at the northern boundary to study the Antarctic Circumpolar Current.

Figure 1. Illustration of (a) $\beta (y)$ and (b) $g(y)$ with $L_y=2{\rm \pi}$ and $l_b=10$ in (2.4).

The coupled system (2.1) conserves energy, potential enstrophy and wave action, which reads

(2.5a) $$\begin{gather} \mathcal{E}=\int \left\{\underset{\mathcal{E}_{QG}}{\underbrace{\frac{1}{2}|\boldsymbol{\nabla} \psi|^2}}+ \underset{\mathcal{E}_{\beta}}{\underbrace{\frac{1}{2}f_0g(y) |\phi|^2}}+ \underset{\mathcal{E}_{NIW}}{\underbrace{\frac{f_0^2}{4} |L\phi|^2}} \right\}\mathrm{d}\boldsymbol{x}, \end{gather}$$

(2.5b)$$\begin{gather}\mathcal{P}=\int q^2\,\mathrm{d}\boldsymbol{x}, \end{gather}$$

(2.5c)$$\begin{gather}\mathcal{A}=\int |\phi|^2\,\mathrm{d}\boldsymbol{x}. \end{gather}$$

Here we decompose the total energy into three components corresponding to the mean, $\beta$-effect and wave potential energy for convenience of later discussion.

In this paper, we focus on the statistically steady states of the NIW–mean flow coupled system with external forcing and dissipation. Following Xie (Reference Xie2020), for a straightforward quantification of energy injection rate to the mean and wave components, forcing and dissipation are added to the mean vorticity equation instead of the QGPV equation and the wave amplitude equation. Therefore we obtain the system

(2.6a)$$\begin{gather} \partial_t\nabla^2\psi+\mathcal{N}(\psi,\phi) = {D}_1\nabla^2\psi+M_{F}k_{f}^{1/2}\mathsf{F}_1, \end{gather}$$

(2.6b)$$\begin{gather}\partial_t\phi+{J}(\psi,\phi)+\mathrm{i} g(y)\phi-\frac{\mathrm{i} f_0}{2}L^2\phi+\frac{\mathrm{i}}{2}\nabla^2\psi\phi = {D}_2\phi+RM_{F}\frac{m}{2N k_{f}^{3/2}}\mathsf{F}_2, \end{gather}$$

where the artificial dissipation operators take the form of ${D}_1=-\alpha _1+\nu _1\nabla ^6$ and ${D}_2=\alpha _2\nabla ^{-2}+\nu _2\nabla ^6$ with constants $\alpha _i$ and $\nu _i$ ($i=1,2$). The external forcings $\mathsf {F}_i$ ($i=1,2$), which are white-noise in time, are isotropic in space and centre around wavenumber $|\boldsymbol {k}|=k_{f}$, where $\boldsymbol {k}=(k,l)$ is the 2-D wavenumber vector in spectral space. Therefore,

(2.7)\begin{equation} \langle\mathsf{F}_i(\boldsymbol{x}_1,t_1) \mathsf{F}_i(\boldsymbol{x}_2,t_2)\rangle= C_0{J}_0(k_{f}|\boldsymbol{x_1}-\boldsymbol{x_2}|)\delta(t_1-t_2), \end{equation}

where $\langle \cdot \rangle$ is an ensemble average, $C_0$ is the covariance and ${J}_0$ is the zeroth-order Bessel function. Here $M_{F}$ and $R$ control the external forcing magnitude and $k_{f}^{1/2}$ and $m/(Nk^{3/2}_{f})$ are the normalized coefficients. Combining (2.5a) and (2.6), we estimate the energy injection rates of QG mean flow and NIW as $\epsilon _{QG}\sim M_{F}^2/2$ and $\epsilon _{NIW}\sim R^2M_{F}^2f_0^2m^2/(16f_0^2m^2+4N^2k_{f}^2)$. Therefore, we define the ratio between energy injection rates of NIW and mean flow as

(2.8)\begin{equation} \lambda =\frac{R^2m^2f_0^2}{8m^2f_0^2+2N^2k_{f}^2}.\end{equation}

Here, the denominator is consistent with the modification by Asselin & Young (Reference Asselin and Young2019), and when this modification is absent, $\lambda = R^2/8$.

We decompose the nonlinear term $\mathcal {N}(\psi,\phi )$ in (2.6a) into different components for convenience in later discussions as

(2.9)\begin{equation} \mathcal{N}(\psi,\phi)=\mathcal{N}_{\beta}(\psi)+ \mathcal{N}_m(\psi)+\mathcal{N}_{W}(\psi,\phi),\end{equation}

where

(2.10a)$$\begin{gather} \mathcal{N}_{\beta}(\psi)=g_y\psi_x+\frac{f_0}{2}g_y(\phi^*_x\phi +\phi^*\phi_x), \end{gather}$$

(2.10b)$$\begin{gather}\mathcal{N}_m(\psi)={J}(\psi,\nabla^2\psi), \end{gather}$$

(2.10c)$$\begin{gather}\mathcal{N}_{W}(\psi,\phi)=\frac{- f_0^2}{4}({J}(\phi^*,L^2\phi) - {J}(L^2\phi^*,\phi)) -\frac{\mathrm{i} f_0^2}{8}\nabla^2(\phi L^2\phi^*-\phi^*L^2\phi)\nonumber\\ -\frac{f_0}{2}\boldsymbol{\nabla} \boldsymbol{\cdot} {J}(\boldsymbol{\nabla} \psi,|\phi|^2). \end{gather}$$

3. Basic properties of the NIW–QG coupled system on a $\beta$-plane

We run the numerical simulations using a Fourier pseudospectral method with $2/3$ dealiasing in space and a resolution of $512\times 512$ in a domain size of $2{\rm \pi} \times 2{\rm \pi}$. We apply a fourth-order explicit Runge–Kutta temporal scheme (cf. Cox & Matthews Reference Cox and Matthews2002) where the linear terms are solved by an integrating factor method, and the nonlinear terms are explicitly approximated. We choose $\beta _0=80$ and $f_0 m/N=32$ in all simulations. The QG mean flows and NIWs are driven by external forcing at a forcing wavenumber $k_{f}=32$. We vary two parameters $M_{F}$ and $R$ to create a matrix of results. Thus, the energy injection rates are $\epsilon _{QG}=8\times 10^{-5},3.2\times 10^{-4},7.2\times 10^{-4}$ with $M_{F}=1000,2000,3000$. Parameter $R$ ranges from $0$ to $3.6$, and the energy injection ratio of NIWs to mean flows $\lambda$ varies from $0$ to $1.296$ by (2.8).

We first run numerical simulations without NIWs as references to study the wave effect. Figure 2 shows a snapshot of the vorticity field $\nabla ^2\psi$ of a NIW-free simulation with $\epsilon _{QG}=8\times 10^{-5}$ at a statistically steady state when jets form. Then we use this vorticity field as the initial condition for other NIW-involved simulations with non-zero $\lambda$.

Figure 2. A snapshot of QG mean flow vorticity fields in the statistically steady state of the NIW-free simulation with $\epsilon _{QG}=8\times 10^{-5}$.

In this section, we focus on qualitative descriptions of characteristic features and leave quantitative studies to §§ 4–6.

3.1. Energy evolution and statistically steady states

In this subsection, we run the coupled system into statistically steady states. We display NIWs’ impact on the total, wave and mean energy evolution. Figure 3(a) shows the evolution of mean energy with different energy injection ratios $\lambda$. When NIWs are injected, the mean energy $\mathcal {E}_{QG}$ decreases and gradually reaches a statistically steady state. With a stronger NIW energy injection rate, the steady value of $\mathcal {E}_{QG}$ becomes weaker, and the relaxation time to reach the statistically steady states increases and then decreases. This non-monotonic dependence implies the existence of a phase transition associated with the disappearance of the jet (cf. figure 4d). In figure 3(b), the $\beta$-related NIW energy, $\mathcal {E}_{\beta }$, oscillates around zero and the amplitude of oscillation increases as $\lambda$ increases. Compared with $\mathcal {E}_{QG}$, $\mathcal {E}_{\beta }$ is relatively small. Figure 3(c) shows that the evolution of wave potential energy $\mathcal {E}_{NIW}$ has three stages. First, there is a sudden wave energy adjustment at a short time interval, then the wave energy increases in a comparable time scale to the mean energy changes and, finally, a statistically steady state is reached. Since $\mathcal {E}_{\beta }$ and $\mathcal {E}_{NIW}$ are small compared with $\mathcal {E}_{QG}$, figure 3(d) shows that the total energy $\mathcal {E}$ is dominated by $\mathcal {E}_{QG}$ but with stronger fluctuations.

Figure 3. Evolutions of (a) QG mean, (b) $\beta$-related, (c) $\mathrm {NIW}$ and (d) total energy (cf. (2.5a)) in simulations with $\lambda =0.016, 0.144, 0.4, 0.784, 1.296$. In (b), for clarity we vertically shift the curves following $\mathcal {E}_{\beta }(\lambda =0.144)+0.02$, $\mathcal {E}_{\beta }(\lambda =0.4)+0.05$, $\mathcal {E}_{\beta }(\lambda =0.784)+0.1$, $\mathcal {E}_{\beta }(\lambda = 1.296)+0.15$. The legend in (a) applies to all four panels.

Figure 4. Snapshots of QG mean vorticity fields, $\nabla ^2\psi$, at statistically steady states with (a) $\lambda =0.144$, (b) $\lambda =0.4$, (c) $\lambda =0.784$ and (d) $\lambda =1.296$. The mean energy injection rate is $\epsilon _{QG}=8\times 10^{-5}$.

Snapshots of statistically steady QG mean flow vorticity field $\nabla ^2\psi$ are presented in figure 4. When the energy injection rate of NIWs is weak (figure 4a with $\lambda =0.144$), zonal jets remain similar in magnitude and shape compared with those in the NIW-free case. As the energy injection ratio of NIW to mean flows increases ($\lambda =0.4$ and $0.784$ shown in figures 4b and 4c), jets weaken and the jet spacing increases. When the energy injection ratio is large enough (figure 4d with $\lambda =1.296$), jets disappear. Compared with the vorticity field without NIWs in figure 2, zonal jets are weakened by NIWs and their geometries are modified, which we study quantitatively in §§ 4–6.

3.2. Energy spectra

To present the impact of NIWs on the spectrum of QG mean flow and, in particular, the jets, based on the 2-D energy spectrum

(3.1)\begin{equation} E(\boldsymbol{k}) =\frac{1}{2}|\boldsymbol{k}|^2 \hat{\psi}(\boldsymbol{k}) \hat{\psi}^*(\boldsymbol{k}), \end{equation}

where the hat $\hat {\cdot }$ denotes the Fourier transform, we define the one-dimensional spectra of the total mean energy, the jet energy and the perturbation energy as

(3.2a)$$\begin{gather} E_T(K) = \sum_{|\boldsymbol{k}|=K-1/2}^{K+1/2} E(\boldsymbol{k}), \end{gather}$$

(3.2b)$$\begin{gather}E_Z(K) = \sum_{|\boldsymbol{k}|=K-1/2}^{K+1/2} E(0,l), \end{gather}$$

(3.2c)$$\begin{gather}E_R(K) = \sum_{|\boldsymbol{k}|=K-1/2}^{K+1/2} E(k\neq0,l). \end{gather}$$

Here, $E_T(K)=E_Z(K)+E_R(K)$. In the classic 2-D $\beta$-plane turbulence, the jet and perturbation energy spectra follow the scalings $E_Z(K)=C_Z\beta ^2K^{-5}$ and $E_R(K)=C_K\epsilon ^{2/3}K^{-5/3}$, where $\epsilon$ is the viscous dissipation rate (Galperin, Sukoriansky & Dikovskaya Reference Galperin, Sukoriansky and Dikovskaya2010). The 2-D potential energy spectrum of NIWs is defined as

(3.3)\begin{equation} E_W(K) =\frac{f_0^2}{4}\sum_{|\boldsymbol{k}|=K-1/2}^{K+1/2} \widehat{{L}\phi}(\boldsymbol{k}) \widehat{{L}\phi}^*(\boldsymbol{k}). \end{equation}

First, we show the spectra of the statistically steady state in the NIW-free simulation in figure 5. At scales between the spectral peak and the forcing wavenumber, the mean and perturbation energy spectra follow $K^{-5/3}$ and $K^{-4}$ scalings, respectively. This $K^{-4}$ scaling was derived by Danilov & Gurarie (Reference Danilov and Gurarie2004) with consideration of a sawtooth profile of zonal mean vorticity, making it different from the classic $K^{-5}$ scaling proposed by Rhines (Reference Rhines1975). To show the wave effect, we present energy spectra with different energy injection ratios $\lambda$ in figure 6. From figure 6(a), we find that as $\lambda$ increases from $0$ to $1.024$, the peaks of $E_Z(K)$ move to smaller wavenumber; in the intermediate scale between spectral peak and the energy injection scale, the jet energy drops but keeps the $K^{-4}$ scaling. However, when $\lambda =1.296$, jets disappear and $E_Z(K)$ deviates from the classic $K^{-5}$ scaling. Figure 6(b) shows that as the NIW energy injection rate gradually becomes comparable with the mean energy injection rate, the perturbation energy drops at the large scale and increases at the small scale; the peak of $E_R(K)$ also moves to larger wavenumbers; and the range of the perturbation energy spectrum with the $K^{-5/3}$ scaling shortens and finally disappears. When jets completely disappear with $\lambda = 1.296$, for wavenumbers smaller than $k_{f}$, $E_R(K)$ follows a $K^{}$ scaling. The combination of jet and perturbation energy spectra, $E_Z(K)$ and $E_R(K)$, respectively, explains the variation of total energy spectrum scaling from $K^{-4}$ to $K^{-5/3}$ and finally to $K^{}$ under the influence of increasing energy injection ratio $\lambda$, which is shown in figure 6(c). However, when the NIW energy injection rate is large enough, such as $\lambda = 1.296$, figures 6(a)–6(c) show that at scales between the spectral peak and the forcing wavenumber, the jet and perturbation spectra differ from the scalings in the theories without wave impact. As to the potential energy spectrum of NIWs in figure 6(d), as the energy injection ratio $\lambda$ increases, for wavenumbers smaller than the forcing wavenumber, it reaches its largest value when $\lambda =1.024$, corresponding to jet disappearance. If we treat the jet strength as an order parameter and the NIW as a perturbation, the increase in perturbation strength implies the existence of a phase transition, which we study in § 4.

Figure 5. The QG mean energy spectrum without NIWs at a statistically steady state with $\epsilon _{QG}=8\times 10^{-5}$.

Figure 6. The NIW-modified QG mean and NIW energy spectra at statistically steady states with $\lambda$ ranging from $0$ to $1.296$: (a) $E_Z(K)$, (b) $E_R(K)$, (c) $E_T(K)$ and (d) $E_W(K)$. The legend in (c) applies to all four panels. The black dash-dotted lines are power functions for reference. The black arrows show the direction of increasing $\lambda$.

3.3. Energy flux across scales

Since the energy flux of the QG mean flows across scales relates to jet formation (Rhines Reference Rhines1975; Vallis & Maltrud Reference Vallis and Maltrud1993; Chekhlov et al. Reference Chekhlov, Orszag, Sukoriansky, Galperin and Staroselsky1996; Sukoriansky, Dikovskaya & Galperin Reference Sukoriansky, Dikovskaya and Galperin2007), this subsection focuses on the QG mean energy flux influenced by NIWs to explain jet weakening and annihilation observed in vorticity fields and mean energy spectra. The spectral energy flux of QG mean flow is defined as

(3.4)\begin{equation} F_{QG}(K)=\sum_{|\boldsymbol{k}|=0}^{K} \frac{-|\boldsymbol{k}|^2}{2}(\hat{\psi}\hat{\mathcal{N}}^*+{\rm c.c.}), \end{equation}

and the positive and negative values of $F_{QG}$ correspond to downscale and upscale fluxes, respectively.

To show the influence of mean flow and NIWs separately, we decompose the energy flux of mean flow as follows:

(3.5a)$$\begin{gather} F_{QG} = F_{QG}^{m} + F_{QG}^{NIW}, \quad \mathrm{with} \end{gather}$$

(3.5b)$$\begin{gather}F_{QG}^{m}=\sum_{|\boldsymbol{k}|=0}^{K}\frac{1}{2}(\hat{\psi} \hat{\mathcal{N}}^*_{m}+ {\rm c.c.}), \end{gather}$$

(3.5c)$$\begin{gather}F_{QG}^{NIW}=\sum_{|\boldsymbol{k}|=0}^{K}\frac{1}{2}(\hat{\psi} \hat{\mathcal{N}}^*_{NIW}+ {\rm c.c.}), \end{gather}$$

where the superscripts ‘$m$’ and ‘${NIW}$’ denote the effects of mean flow and NIWs, respectively.

When NIWs are absent, the coupled system (2.1) is identical to the 2-D turbulence system on a $\beta$-plane; therefore, the mean energy should transfer upscale with a constant value, which is shown in figure 7. Here, the finite hyperviscosity brings about a weak downscale flux.

Figure 7. The QG mean energy flux without NIWs at a statistically steady state with $k_{f}=32$ and $\epsilon _{QG}=8\times 10^{-5}$. The negative and positive values denote the upscale and downscale energy fluxes, respectively.

In figure 8 we present the energy flux of mean flow with different energy injection rates of NIWs. The NIWs induce bidirectional flux of mean-flow energy, and as $\lambda$ increases, the portion of downscale flux increases, as shown in figures 8(a)–8(c). Therein, both the energy flux induced by mean flow, $F_{QG}^{m}$, and the energy flux induced by NIWs, $F_{QG}^{NIW}$, increase and the cancellation of these two parts leaves a bidirectional mean-flow energy flux. When $\lambda$ is large enough (cf. figure 8d), all the mean energy transfers downscale. The presence of mean-flow inertial ranges with constant energy fluxes above and below the forcing wavenumber implies negligible conversion between QG mean flow and NIWs, which is consistent with the catalytic effect of NIWs on the mean energy flux (Xie Reference Xie2020; Barkan et al. Reference Barkan, Srinivasan, Yang, McWilliams, Gula and Vic2021).

Figure 8. Energy flux of QG mean flow in NIW-containing cases of (a) $\lambda =0.144$, (b) $\lambda =0.4$, (c) $\lambda =0.784$ and (c) $\lambda =1.296$ with $k_{f}=32$ and $\epsilon _{QG}=8\times 10^{-5}$.

4. Dependence of jet magnitude on the strength of NIWs

In this section, we study the dependence of jet strength on the strength of NIWs, which $\lambda$ parametrizes. We define jet energy as $\mathcal {E}_{jet}=\int \bar {u}^2/2\,\mathrm {d} \boldsymbol {x}$, where $\bar {u}$ is the zonal-averaged zonal QG mean velocity. Figure 9(a) shows the dependence of $\mathcal {E}_{jet}/\mathcal {E}_{QG}$ on $\lambda$ at statistically steady states with different mean energy injection rates. When the mean energy injection rate, $\epsilon _{QG}$, is fixed, there exists a critical value $\lambda _c$: when $0<\lambda <\lambda _c$, $\mathcal {E}_{jet}/\mathcal {E}_{QG}$ is almost a constant, and when $\lambda$ goes across the critical points $\lambda _c$, it abruptly drops to a value close to $0$. Therefore, as the NIW energy injection rate increases from $0$, the QG–NIW coupled system experiences a first-order phase transition, across which jets suddenly disappear.

Figure 9. Ratio between jet energy and the QG mean-flow energy with different mean energy injection rates. In (b), the horizontal coordinate is normalized as $\tilde {\lambda }=\lambda (\epsilon _{QG}/\epsilon _{QG}^0)^{1/2}$.

Next, we explore the dependence of QG mean energy injection rate $\epsilon _{QG}$ on critical values $\lambda _c$. Since near the critical point the jet magnitude is small, in the mean vorticity equation (2.6a), by ignoring the mean-flow impact, we obtain the quasilinear balance between nonlinear wave term in (2.10c) and the external forcing, i.e. $\mathcal {N}_{W}\sim M_{F}k_{f}^{1/2}\mathsf {F}_1$, leading to $|\phi |^2\sim M_{F}k_{f}^{-1/2}f_0^{-2}$. The wave magnitude can be estimated by $|\phi |\sim RM_{F} m N^{-1} f_0^{-1}k_{f}^{-3/2}$ from the wave equation (2.6b). Combining them we obtain $R^2 \sim k_{f}^{3/2}N^2M_{F}^{-1}m^{-2}$. Since $\lambda \sim R^2$ and $\epsilon _{QG}\sim M_{F}^2$, the critical value scales as $\lambda _c \sim \epsilon _{QG}^{-1/2}N^2k_{f}f_0^{-1/2}m^{-2}$. Therefore in figure 9(b), we introduce a normalized parameter $\tilde {\lambda }=\lambda (\epsilon _{QG}/\epsilon _{QG}^0)^{1/2}$ with $\epsilon _{QG}^0=8\times 10^{-5}$ to collapse data with different mean energy injection rates. Hence, we find that the quasilinear dynamics controls the phase transition of the NIW–QG coupled system.

Figure 10 shows the $\lambda$ dependence of the ratio between the upscale mean energy flux and mean energy injection rate, $\epsilon _{up}/\epsilon _{QG}$. Here, $\epsilon _{up} = \max _{K< k_{f}} |F_{QG}|$ is the upscale mean energy flux. In figure 10(a), we present the results with three mean energy injection rates. The upscale energy flux ratio gradually decreases as $\lambda$ increases and finally shows a second-order phase transition across $\lambda _c$, which contrasts with the behaviour of $\mathcal {E}_{jet}/\mathcal {E}_{QG}$ shown in figure 9, but these two quantities share the same critical values. This behaviour is similar to that on an $f$-plane (Xie Reference Xie2020). Figure 10(b) shows that after normalizing $\lambda$ by $\epsilon _{QG}^{-1/2}$, curves with different mean energy injection rates collapse. This implies that the above-mentioned quasilinear argument controls the phase transition of the NIW–QG coupled system. The concurrence of the disappearance of jet and upscale energy flux implies that we can understand jet annihilation from the perspective of NIW-modified mean energy flux. The injection of NIWs induced downscale flux of mean energy (Xie Reference Xie2020); therefore, the upscale energy flux decreases and the jets weaken. Thus, from the same perspective, other properties, such as jet spacing, should also be controlled by upscale energy flux, which we study in detail in § 5.

Figure 10. Ratio between the upscale energy flux and the total mean energy injection rate with different mean energy injection rates. In (b), the horizontal coordinate is normalized as $\tilde {\lambda }=\lambda (\epsilon _{QG}/\epsilon _{QG}^0)^{1/2}$.

5. Jet spacing

Another important character of jets is their spacing, $L_{jet}$. We may define $L_{jet}$ from the jet energy spectrum:

(5.1)\begin{equation} L_{jet}=\frac{\displaystyle \int E_Z(k)\, {\rm d} k}{\displaystyle \int k E_Z(k)\, {\rm d} k}. \end{equation}

Proposed by Rhines (Reference Rhines1975), jet spacing can be captured by a characteristic scale $L_{Rh}=\sqrt {{U}/{\beta }}$, which is named the Rhines scale. By choosing the r.m.s. of jet velocity, $U_{jet} =\sqrt {2\mathcal {E}_{jet}}$, as the characteristic velocity, the Rhines scale can be expressed as

(5.2)\begin{equation} L_{Rh}=\frac{(2\mathcal{E}_{jet})^{1/4}}{\beta^{1/2}}. \end{equation}

In 2-D $\beta$-plane turbulence, energy transfers upscale, and with further frictional damping, the kinetic energy can be obtained from the balance between energy injection and dissipation as $\mathcal {E} = \epsilon /(2\alpha )$ in a statistically steady state. In our system, the QG mean energy transfers bidirectionally. Since the jet is a large-scale structure, one may propose that the upscale energy flux $\epsilon _{up}$ is important for jets (cf. Scott & Dritschel Reference Scott and Dritschel2012). Therefore, following the above-mentioned classic 2-D turbulence argument, we obtain a characteristic scale:

(5.3)\begin{equation} L_{\alpha}=\frac{\epsilon_{up}^{1/4}}{\alpha^{1/4}\beta^{1/2}}. \end{equation}

We further consider that in the NIW–QG coupled system, QG mean energy flux $F_{QG}$ can be decomposed into the mean part $F_{QG}^m$ and the NIW-related flux $F_{QG}^{NIW}$ (cf. (3.5)). Since a jet is a zonal-mean quantity, we conjecture that only $F_{QG}^m$ is relevant to the jet spacing. Therefore, we propose a characteristic scale estimated by the upscale energy flux induced by mean flow only, $\epsilon _{adv} = \max _{K< k_{f}} |F_{QG}^m|$:

(5.4)\begin{equation} L_{adv}=\frac{\epsilon_{adv}^{1/4}}{\alpha^{1/4}\beta^{1/2}}. \end{equation}

When no NIWs are injected into the coupled system (i.e. $\lambda =0$), ${L}_{adv}$ equals $L_{\alpha }$ since $\epsilon _{up} = \epsilon _{adv}$.

We show the normalized jet spacing and three normalized characteristic scales in figure 11(a) with $\lambda$ ranging from $0$ to $\lambda _c$. As $\lambda$ increases, jet spacing increases. However, both $L_{Rh}$ and $L_{\alpha }$ decrease, so they cannot describe jet spacing. Interestingly, even though the presence of NIWs reduces the upscale flux of the mean energy, the jet spacing increases with $\lambda$ increasing. In contrast, ${L}_{adv}$ increases as $\lambda$ increases, which can be seen from the increase of $F_{QG}^m$ (cf. figure 6). Figure 11(b) shows that the jet spacing is well described by ${L}_{adv}$, which implies that the upscale energy flux induced by mean advection controls jet formation. Thus, by measuring the mean quantities only, we can capture the jet spacing using the classic Rhines scale as if the NIW is absent, which again points out the catalytic role of NIWs.

Figure 11. (a) Jet spacings and different definitions of Rhines scale with $k_{f}=32$, $\epsilon _{QG}=8\times 10^{-5}$ and varying $\lambda$. The tilde denotes the scale normalized by NIW-free jet spacing ($\lambda =0$). (b) Comparison between the numerically obtained jet spacing $\tilde {L}_{jet}$ and the mean-advection-related Rhines scale $\tilde {L}_{adv}$ with different $\lambda$. In these simulations $\epsilon _{QG}=8\times 10^{-5}$.

6. Near-inertial-wave-induced symmetry breaking

The NIWs break the temporal symmetry by rotating in the direction opposite to that of the background rotation. This temporal symmetry breaking leads to spatial asymmetry, which we study in this section from the perspectives of mean-velocity profile and jet meridional migration.

6.1. Mean-velocity profile

In this subsection, we focus on the mean-velocity profile of zonal jets. We show the jet velocity profile $\bar {u}$ obtained by taking a zonal average of zonal velocity $u$ in figure 12(a). Similar to the classic $\beta$-plane jet profile, the zonal velocity has an east–west asymmetry: the east jets have cusps while the west jets are smooth. In figure 12(a), the presence of NIWs introduces an additional asymmetry of the jet with steeper left flanks (red rectangles) than the right flanks (blue rectangles). This asymmetry manifests itself in the asymmetric profile of zonal-mean vorticity $\bar {Q}=-\bar {u}_y$ with a greater magnitude of positive vorticity, which is shown in figure 12(b). Figure 12(c) shows that the zonal-mean wave action $|\bar {\phi }|^2$ concentrates around regions with negative vorticity, $\bar {Q} < 0$, which is consistent with NIW concentration on anticyclones (Kunze Reference Kunze1985; Danioux et al. Reference Danioux, Vanneste and Bühler2015). Since NIWs concentrate on the right flank with negative vorticity, the intensity of NIW–QG interaction there should be more potent than that at the left flank with positive vorticity. The NIW–QG interaction induces a downscale mean-energy flux and increases dissipation at small scales (cf. § 3). Therefore, the mean shear on the left flank is stronger than that on the right flank.

Figure 12. Zonal-mean profiles of (a) zonal velocity, (b) vorticity and (c) wave action with $\epsilon _{QG}=8\times 10^{-5}, \lambda =0.4$ at $t=2000$. The blue and red rectangles denote the regions with positive and negative vorticity, respectively.

To show the impact of NIW strength on mean-velocity profiles, in figure 13 we compare jet profiles with different $\lambda$ by moving the profiles to the same peak location. It shows that the asymmetry of jet flanks enhances as $\lambda$ increases. Since jets are annihilated when $\lambda >\lambda _c$, we only study the jet profile in the parameter regime with $\lambda <\lambda _c$.

Figure 13. Jet velocity profiles with $\epsilon _{QG}=8\times 10^{-5}$ and changing $\lambda$. Here, $\tilde {u}$ denotes the moved profiles sharing the same peak location for comparison; $\tilde {y}$ denotes the moved meridional coordinate.

To quantify the jet flank asymmetry, we study the skewness of the zonal-mean vorticity $\bar {Q}$. We define $Q_{3}=\langle \overline {Q}^3\rangle$ and $Q_{2}=\langle \overline {Q}^2\rangle$, and use the skewness $S_{k}=Q_{3}/Q_{2}^{3/2}$ to quantify the asymmetry. We present $S_{k}$ with different $\epsilon _{QG}$ and $\lambda$ in figure 14(a). The skewness of the mean vorticity first increases with $\lambda$ due to NIWs concentrating on the negative vorticity, and then it suddenly drops as $\lambda$ reaches $\lambda _c$ because jets are annihilated. Similar to figures 9(b) and 10(b), figure 14(b) shows that the same normalization $\tilde {\lambda }=\lambda (\epsilon _{QG}/\epsilon _{QG}^0)^{1/2}$ also works for $S_{k}$, which again reveals the significance of quasilinear dynamics (cf. § 4).

Figure 14. The skewness of zonal-mean vorticity with different mean energy injection rates. In (b) the horizontal coordinate is normalized as $\tilde {\lambda }=\lambda (\epsilon _{QG}/\epsilon _{QG}^0)^{1/2}$.

Taken together, these results suggest that the mean flow on the right flank is reduced more than that on the left flank due to NIWs concentrating on the right flanks, where anticyclones reside and NIW–mean interaction grows. Therefore, NIWs catalyse the asymmetry of jet flanks. Besides, in the current parameter regime, we find that the Stokes drift (Bühler Reference Bühler2009) is negligible compared with the Lagrangian mean velocity. Details are presented in Appendix A.

6.2. Jet migration

In this subsection, we focus on the impact of NIWs on jet migration. We show the Hovmöller diagram of the zonal-mean zonal velocity, $\bar {u}$, with $\lambda =0.4$ and $\epsilon _{QG}=8\times 10^{-5}$ in figure 15, where jets migrate with a constant speed in the meridional direction. From this Hovmöller diagram, we can define a jet migration speed $\mathcal {V}=\delta y/\delta t$, where $\delta y$ is the latitudinal shift and $\delta t$ is the elapsed time.

Figure 15. Hovmöller diagram of the zonal-mean zonal velocity with $\epsilon _{QG}=8\times 10^{-5}$ and $\lambda =0.4$.

Figure 16 shows the dependence of the measured jet migration speed on the energy injection ratio. The zonal jets without NIWs stay static, and as the NIW energy injection rate increases, jet meridional migration speed gradually increases. To find the mechanism behind zonal jet migration, we start from the zonal-mean momentum equation derived from (2.6a):

(6.1)\begin{equation} \partial_t \bar{u} +\partial_y \overline{u'v'} -\bar{\varLambda} ={-}\alpha_1 \bar{u}+\nu_1 \frac{\partial^6\bar{u}}{\partial y^{6}},\end{equation}

where $(u',v')$ is the fluctuation with respect to the zonal-mean velocity (cf. Chemke & Kaspi Reference Chemke and Kaspi2015b). Parameter $\varLambda$ is derived from the wave-related term in (2.10c) with $\hat {\varLambda }=-\mathrm {i}\hat {\mathcal {N}}_{NIW}/l$ in Fourier space. Using Taylor's identity

(6.2)\begin{equation} -\partial_y ( \overline{u'v'})=\overline{v'Q'}, \end{equation}

equation (6.1) becomes

(6.3)\begin{equation} \partial_t \bar{u} =\overline{v'Q'} + \bar{\varLambda} -\alpha_1 \bar{u} +\nu_1 \frac{\partial^6 \bar{u}}{\partial y^{6}}.\end{equation}

Here, $Q'= Q-\bar {Q}$. Assuming that jets migrate in the meridional direction with a constant speed, we consider a frame of reference moving with jets by introducing a coordinate $\varphi =y-\mathcal {V}^*t$ with $\mathcal {V}^*$ the jet migration speed. Thus (6.3) becomes

(6.4)\begin{equation} -\mathcal{V}^*\frac{\mathrm{d} \bar{u}}{\mathrm{d} \varphi}= \underset{\mathcal{T}(\varphi)}{\underbrace{ \overline{v'Q'} + \bar{\varLambda} -\alpha_1 \bar{u} +\nu_1 \frac{\partial^6 \bar{u}}{\partial y^{6}}}},\end{equation}

where we introduce the momentum flux convergence $\mathcal {T}$ for convenience.

Figure 16. Dependence of jet migration speed $\mathcal {V}$ on $\lambda$ with $\epsilon _{QG}=8\times 10^{-5}$.

In figure 17, we show the NIW-induced asymmetric jet profile (cf. § 6.1) and its corresponding momentum flux convergence $\mathcal {T}$. We define the width of the jet flank as the meridional range with $\bar {u}>0$, i.e. $\varphi _0^+-\varphi _0^-$, and $\varphi _{jet}$ is the latitude of the maximum zonal-mean zonal velocity $\bar {u}_{max}$. Integrating (6.4) from $\varphi _0^-$ to $\varphi _{jet}$ and $\varphi _{jet}$ to $\varphi _0^+$ and then subtracting the resultant equations, the jet migration speed $\mathcal {V}^*$ is expressed as

(6.5)\begin{equation} \mathcal{V}^*=\frac{1}{2\bar{u}_{max}} \left( \int^{\varphi_0^+}_{\varphi_{jet}} \mathcal{T}\,\mathrm{d} \varphi - \int^{\varphi_{jet}}_{\varphi_0^-} \mathcal{T}\,\mathrm{d} \varphi \right).\end{equation}

Figure 17. Profiles of the normalized zonal-mean velocity $\bar {u}/\bar {u}_{max}$ and the normalized momentum flux convergence $\mathcal {T}/\mathcal {T}_{max}$ with $\epsilon _{QG}=8\times 10^{-5}$ and $\lambda =0.4$.

To find the dominant effect on jet migration, we introduce the following decomposition:

(6.6)\begin{equation} \int^{\varphi_0^+}_{\varphi_{jet}}\mathcal{T}\,\mathrm{d} \varphi -\int^{\varphi_{jet}}_{\varphi_0^-}\mathcal{T}\,\mathrm{d} \varphi =\mathcal{I}_{adv}+\mathcal{I}_{NIW} +\mathcal{I}_{\alpha}+\mathcal{I}_{\nu},\end{equation}

where

(6.7a)$$\begin{gather} \mathcal{I}_{adv}=\int^{\varphi_0^+}_{\varphi_{jet}}\overline{v'Q'}\,\mathrm{d} \varphi - \int^{\varphi_{jet}}_{\varphi_0^-}\overline{v'Q'}\,\mathrm{d} \varphi, \end{gather}$$

(6.7b)$$\begin{gather}\mathcal{I}_{NIW}=\int^{\varphi_0^+}_{\varphi_{jet}}\bar{\varLambda}\,\mathrm{d} \varphi - \int^{\varphi_{jet}}_{\varphi_0^-}\bar{\varLambda}\,\mathrm{d} \varphi, \end{gather}$$

(6.7c)$$\begin{gather}\mathcal{I}_{\alpha}=\int^{\varphi_0^+}_{\varphi_{jet}}-\alpha_1 \bar{u}\,\mathrm{d} \varphi + \int^{\varphi_{jet}}_{\varphi_0^-}\alpha_1 \bar{u}\,\mathrm{d} \varphi, \end{gather}$$

(6.7d)$$\begin{gather}\mathcal{I}_{\nu}=\int^{\varphi_0^+}_{\varphi_{jet}}\nu_1 \frac{\partial^6\bar{u}}{\partial \varphi^{6}} \,\mathrm{d} \varphi - \int^{\varphi_{jet}}_{\varphi_0^-}\nu_1 \frac{\partial^6\bar{u}}{\partial \varphi^{6}} \,\mathrm{d} \varphi. \end{gather}$$

Figure 18(a) shows the $\lambda$ dependence of these terms. We observe that the total net momentum flux, $\mathcal {I}_{total}=\mathcal {I}_{adv}+\mathcal {I}_{NIW}+ \mathcal {I}_{\alpha }+ \mathcal {I}_{\nu }$, increases with $\lambda$. Therein, the dissipation-related net momentum fluxes, $\mathcal {I}_{\alpha }$ and $\mathcal {I}_{\nu }$, are negligible. The advection-related net momentum flux, $\mathcal {I}_{adv}$, is dominant and increases with $\lambda$, while the wave-related net momentum flux, $\mathcal {I}_{NIW}$, has an opposite sign to $\mathcal {I}_{total}$ and therefore prohibits jet migration. Figure 18(b) shows that when $\lambda$ ranges from $0$ to $\lambda _c$, the jet migration speeds calculated from (6.5) capture the migration speeds obtained from numerical simulations. The impact of NIWs on jet migration again shows its catalytic effect – the direct effect of NIWs, $\mathcal {I}_{NIW}$, restrains jet migration, but NIW-modified mean flow dominates and drives the migration.

Figure 18. (a) The dependences of the net momentum flux $\mathcal {I}$ on $\lambda$. (b) Numerically obtained jet migration speed $\mathcal {V}$ against the migration speed $\mathcal {V}^*$ calculated from momentum balance (6.5). Here, $\epsilon _{QG}=8\times 10^{-5}$.

The direction of jet migration can be understood as a consequence of NIW concentration in anticyclones (Danioux et al. Reference Danioux, Vanneste and Bühler2015). Since NIWs accumulate in anticyclones, more energy is transferred downscale by negative vorticity and finally dissipates (cf. § 3), so the negative vorticity of the mean flow is weaker than the positive vorticity. Then the vortex pairs propagate poleward. This phenomenon was observed in simulations of both primitive equations (Hernandez-Duenas, Smith & Stechmann Reference Hernandez-Duenas, Smith and Stechmann2014) and reduced model (Xie & Vanneste Reference Xie and Vanneste2015). Similarly, figure 19 shows that both jets and waves migrate poleward and NIWs concentrate in regions with negative vorticity, so jets migrate poleward.

Figure 19. Hovmöller diagram of the zonal-mean (a) vorticity and (b) NIW amplitude with $\epsilon _{QG}=8\times 10^{-5}$ and $\lambda =0.4$.

We need to note that the coupled model was initially derived in the ocean context. But most jets in ocean models are found to propagate equatorward (Chan et al. Reference Chan, Plumb and Cerovecki2007; Ashkenazy & Tziperman Reference Ashkenazy and Tziperman2016). Nevertheless, poleward migration was also observed over a short range and less systematically (Chan et al. Reference Chan, Plumb and Cerovecki2007). Thus, the impact of NIWs cannot fully explain the migration of oceanic jets.

7. Summary and discussion

Using a 2-D single-vertical-mode NIW–QG coupled model on a $\beta$-plane (2.1), we study the impact of NIWs on zonal jets in forced-dissipative statistically steady states. We explore the parameter regime with different mean energy injection rates and the ratio between NIW and QG mean energy injection rate, which $\lambda$ parametrizes.

The presence of NIWs leads to bidirectional mean-flow energy transfer, which shows a second-order phase transition with a critical value $\lambda _c$. When $0<\lambda <\lambda _c$, the energy of balanced flow transfers to both large and small scales. As $\lambda$ increases, the portion of upscale transferring energy decreases, while the proportion of zonal-jet energy in the energy of the balanced flow remains almost unchanged. When $\lambda > \lambda _c$, the QG inverse cascade is eliminated and zonal jets cannot form. Since there is a negligible energy exchange between mean flows and NIWs, the existence of the catalytic wave induction mechanism initially proposed on an $f$-plane (Xie Reference Xie2020) is justified on a $\beta$-plane. This mechanism differs from the spontaneous and stimulated loss of balance (Vanneste Reference Vanneste2013; Xie & Vanneste Reference Xie and Vanneste2015), where the mean energy converts to that of NIWs and finally dissipates.

Based on the picture that the $\beta$-effect arrests the upscale energy flux, previous studies (Rhines Reference Rhines1975; Williams Reference Williams1978) showed that the Rhines scale could capture the zonal jet spacing. Nevertheless, we find that the Rhines scale calculated from the r.m.s. velocity (Rhines Reference Rhines1975) and the upscale energy transfer rate $\epsilon _{up}$ (Vallis Reference Vallis2006) cannot sufficiently describe jet spacing when NIWs are present. Instead, the Rhines scale calculated from the upscale-mean-energy transfer rate induced by the mean advection, $L_{adv}=(2\epsilon _{adv})^{1/4}\alpha ^{-1/4}\beta ^{-1/2}$, captures the jet spacing before jet total annihilation by NIWs.

The concentration of NIWs in anticyclones induces symmetry breaking of zonal jets. The NIWs lead to an asymmetry of the left and right flanks about the jet profile cusp, where the right flank's slope is flatter than that of the left flank. Moreover, this asymmetry enhances as the energy injection ratio $\lambda$ increases. The concentration of NIWs in anticyclones or regions with negative vorticity (Danioux et al. Reference Danioux, Vanneste and Bühler2015) can explain the mechanism behind such asymmetry. Therefore at the right flank of a jet cusp, stronger NIW–mean flow interaction leads to stronger downscale mean energy flux induced by waves, which enhances the mean flow dissipation, so the right flank is flatter than the left flank. The NIW-induced symmetry breaking also manifests in jets’ poleward migration. The interaction between NIWs and mean flows also influences the jet migration speed in the meridional direction. Based on an expression of migration speed calculated from the momentum equation in a frame of reference moving with jets (Chemke & Kaspi Reference Chemke and Kaspi2015b), we find that the distribution of the momentum flux convergence to the cusps in the velocity profile is asymmetric, which leads to jet migration. In this expression, the mean-advection net momentum flux dominates, while the NIW-related term impedes jet migration but is of secondary importance. Furthermore, we numerically find that the migration speed increases as the energy injection of NIWs increases. Even though the asymmetric jet profile and meridional migration coexist, we do not know if this coexistence is inevitable. Such coexistence is also observed in jet migration without wave impact (cf. Cope Reference Cope2020), so there is more to be understood on jet asymmetry.

Considering that the strong vortex regions are close to the peak of the cusp of eastward jets, this poleward migration mimics the left turning of vortex dipoles (Hernandez-Duenas et al. Reference Hernandez-Duenas, Smith and Stechmann2014; Xie & Vanneste Reference Xie and Vanneste2015). Even though the reduced model is derived in an oceanic context (Xie & Vanneste Reference Xie and Vanneste2015), the jet migration direction coincides with the results of observations (Riehl et al. Reference Riehl, Yeh and La Seur1950; Dickey et al. Reference Dickey, Marcus and Hide1992) and numerical simulations (James & Dodd Reference James and Dodd1996; Chemke & Kaspi Reference Chemke and Kaspi2015b) in the atmosphere. However, in the ocean, poleward migration of zonal flow is observed less frequently (Chan et al. Reference Chan, Plumb and Cerovecki2007), implying that the impact of NIWs on zonal jets may not be a dominant effect that determines the direction of oceanic jet migration.

The catalytic wave impact on zonal jets manifests itself from three aspects. (i) In the mean-flow bidirectional energy transfer across scales, NIWs impact the relative strength of the downscale flux without energy conversion between NIW and mean flow. (ii) The jet spacing is described by the Rhines scale obtained from the upscale mean energy flux induced by mean flow alone as if there is no NIW. (iii) The jet meridional migration speed is dominated by the asymmetry of mean momentum flux induced by mean flow, but in our set-up, this asymmetry is induced by NIWs. We summarize the effects of NIW concentration on the right flank, the asymmetry of the jet and the meridional migration in figure 20. However, we do not understand why these wave effects are catalytic, which is an interesting future direction.

Figure 20. Illustration of the jet–NIW interaction.

We close this paper by pointing out several limitations of the present study. To solve the coupled NIW–QG model with the pseudospectral method, we prescribe a periodic $\beta (y)$ (cf. 2.4). Therefore, around one-half of the meridional extension approximates a $\beta$-plane, while other regions are only introduced to construct a periodic domain for the spectral method. To avoid the Gibbs effect, the $\beta$-effect is artificially prescribed with a narrow but finite range of non-constant value, which leads to unphysical effects, but we assume that this effect is only confined in this narrow region, which is the drawback of this set-up. In addition, the NIW–QG model in this paper is asymptotically derived in the regime of small Rossby and small Burger numbers (cf. Xie & Vanneste Reference Xie and Vanneste2015), which limits the validity of the current results; therefore, it would be interesting to check if our results persist in primitive equations in future work. Besides, we applied a simple turbulence phenomenology based on energy flux to explain jet formation, but we did not distinguish the local and non-local effects in spectral space. So in future work, we plan to study the impact of waves on jets by applying a quasilinear approximation to the balanced motion (cf. Marston & Tobias Reference Marston and Tobias2023).