2.1 Formulation

The expansion is based on a distinguished limit in which the Burger number, denoted in this section by  $\unicode[STIX]{x1D716}^{2}$ , is small. Following Young & Ben Jelloul (Reference Young and Ben Jelloul1997) one uses the complex representation

(2.1a,b ) $$\begin{eqnarray}{\mathcal{U}}\stackrel{\mathit{def}}{=}u+\text{i}v,\quad \text{and}\quad s\stackrel{\mathit{def}}{=}x+\text{i}y.\end{eqnarray}$$

It follows that

(2.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{s}={\textstyle \frac{1}{2}}(\unicode[STIX]{x2202}_{x}-\text{i}\unicode[STIX]{x2202}_{y}),\quad \unicode[STIX]{x2202}_{s^{\ast }}={\textstyle \frac{1}{2}}(\unicode[STIX]{x2202}_{x}+\text{i}\unicode[STIX]{x2202}_{y}),\end{eqnarray}$$

and the horizontal Laplacian is

(2.3) $$\begin{eqnarray}\triangle =\unicode[STIX]{x2202}_{x}^{2}+\unicode[STIX]{x2202}_{y}^{2}=4\unicode[STIX]{x2202}_{s}\unicode[STIX]{x2202}_{s^{\ast }}.\end{eqnarray}$$

Other useful results are $u_{x}+v_{y}={\mathcal{U}}_{s}+{\mathcal{U}}_{s^{\ast }}^{\ast }$ and $v_{x}-u_{y}=\text{i}{\mathcal{U}}_{s^{\ast }}^{\ast }-\text{i}{\mathcal{U}}_{s}$ .

We employ with a multiple time scale expansion

(2.4) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}=\unicode[STIX]{x2202}_{t_{0}}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{2}}+\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x2202}_{t_{4}}+\cdots\end{eqnarray}$$

Then, with the complex notation outlined in (2.1) through (2.3), the scaled horizontal momentum equation is

(2.5) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{t_{0}}+\text{i}f){\mathcal{U}}=-\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{2}}{\mathcal{U}}-\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x2202}_{t_{4}}{\mathcal{U}}-\unicode[STIX]{x1D716}^{2}2p_{s^{\ast }}+\text{ord}(\unicode[STIX]{x1D716}^{6}).\end{eqnarray}$$

The other equations of motion are

(2.6) $$\begin{eqnarray}\displaystyle & {\mathcal{U}}_{s}+{\mathcal{U}}_{s^{\ast }}^{\ast }+w_{z}=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t_{0}}b+wN^{2}=-\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{2}}b+\text{ord}(\unicode[STIX]{x1D716}^{4}), & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & p_{z}-b=0. & \displaystyle\end{eqnarray}$$

All fields are expanded in powers of $\unicode[STIX]{x1D716}^{2}$ : ${\mathcal{U}}={\mathcal{U}}_{0}+\unicode[STIX]{x1D716}^{2}{\mathcal{U}}_{2}+\unicode[STIX]{x1D716}^{4}{\mathcal{U}}_{4}+\cdots$ and so on.

2.2 The expansion

At leading order, $\unicode[STIX]{x1D716}^{0}$ , the pressure is negligible in (2.5) and the general solution is

(2.9) $$\begin{eqnarray}{\mathcal{U}}_{0}=\mathsf{L}\,A\text{e}^{-\text{i}ft_{0}},\end{eqnarray}$$

where $\mathsf{L}$ is the operator in (1.3). Using $\mathsf{L}\,A$ in (2.9) facilitates expression of the other leading-order fields in terms of  $A$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle w_{0}=-\frac{f^{2}}{N^{2}}A_{zs}\text{e}^{-\text{i}ft_{0}}+\text{cc}, & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & b_{0}=\text{i}fA_{zs}\text{e}^{-\text{i}ft_{0}}+\text{cc}, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & p_{0}=\text{i}fA_{s}\text{e}^{-\text{i}ft_{0}}+\text{cc}, & \displaystyle\end{eqnarray}$$

where $\text{cc}$ denotes complex conjugate.

At order $\unicode[STIX]{x1D716}^{2}$ , the horizontal momentum equation is

(2.13) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{t_{0}}+\text{i}f){\mathcal{U}}_{2}=-\unicode[STIX]{x2202}_{t_{2}}\mathsf{L}\,A\text{e}^{-\text{i}ft_{0}}-2\text{i}fA_{ss^{\ast }}\text{e}^{-\text{i}ft_{0}}+2\text{i}fA_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}.\end{eqnarray}$$

Removing the resonant terms gives the leading-order evolution equation

(2.14) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t_{2}}\mathsf{L}\,A+2\text{i}fA_{ss^{\ast }}=0.\end{eqnarray}$$

Solving the remnants of (2.13) one obtains

(2.15) $$\begin{eqnarray}{\mathcal{U}}_{2}=A_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}.\end{eqnarray}$$

The ‘reconstituted’ velocity, ${\mathcal{U}}_{0}+\unicode[STIX]{x1D716}^{2}{\mathcal{U}}_{2}$ , has the slightly elliptic hodograph characteristic of NIWs.

The second-order vertical velocity follows from (2.6)

(2.16) $$\begin{eqnarray}w_{2z}=-A_{sss^{\ast }}\text{e}^{-\text{i}ft_{0}}+\text{cc}.\end{eqnarray}$$

Now we use the leading-order evolution equation in (2.14) to rewrite (2.16) as

(2.17) $$\begin{eqnarray}2\text{i}fw_{2z}=\unicode[STIX]{x2202}_{t_{2}}\mathsf{L}\,A_{s}\text{e}^{-\text{i}ft_{0}}+\text{cc}.\end{eqnarray}$$

This reduces the number of $x$ and $y$ derivatives from three in (2.16) to one in (2.17). Thus YBJ⁺ avoids the high-order differential operator, $\triangle ^{2}$ , in (1.12). Another advantage of (2.17) is that one can now integrate with respect to $z$ to obtain the second-order vertical velocity

(2.18) $$\begin{eqnarray}2\text{i}fw_{2}=\frac{f^{2}}{N^{2}}\unicode[STIX]{x2202}_{t_{2}}A_{zs}\text{e}^{-\text{i}ft_{0}}+\text{cc}.\end{eqnarray}$$

One rejoices that the boundary condition $A_{z}=0$ ensures that both $w_{0}$ and $w_{2}$ vanish. This integration in $z$ avoids the inverse operator $\mathsf{L}^{-1}$ that complicates (1.12). For all these reasons, the transition from (2.16) to (2.17) is a crucial technical step in the derivation of YBJ⁺.

With $w_{2}$ in hand, the other second-order fields follow

(2.19a,b ) $$\begin{eqnarray}b_{2}={\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{t_{2}}A_{zs}\text{e}^{-\text{i}ft_{0}}+\text{cc},\quad \text{and}\quad p_{2}={\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{t_{2}}A_{s}\text{e}^{-\text{i}ft_{0}}+\text{cc}.\end{eqnarray}$$

At fourth order, the horizontal momentum equation is

(2.20) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x2202}_{t_{0}}+\text{i}f){\mathcal{U}}_{4} & = & \displaystyle -\unicode[STIX]{x2202}_{t_{4}}{\mathcal{U}}_{0}-\unicode[STIX]{x2202}_{t_{2}}{\mathcal{U}}_{2}-2p_{2s^{\ast }},\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & = & \displaystyle -\unicode[STIX]{x2202}_{t_{4}}\mathsf{L}\,A\text{e}^{-\text{i}ft_{0}}-\unicode[STIX]{x2202}_{t_{2}}A_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}-\unicode[STIX]{x2202}_{t_{2}}A_{ss^{\ast }}\text{e}^{-\text{i}ft_{0}}-\unicode[STIX]{x2202}_{t_{2}}A_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}.\end{eqnarray}$$

To remove the resonant terms from (2.21), we require

(2.22) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t_{4}}\mathsf{L}\,A+\unicode[STIX]{x2202}_{t_{2}}A_{ss^{\ast }}=0.\end{eqnarray}$$

2.3 Reconstitution and remodelling

Following the philosophy of Roberts (Reference Roberts1985), we now ‘reconstitute’ the perturbation series by forming the linear combination (2.14)  $+$   $\unicode[STIX]{x1D716}^{2}$ (2.22):

(2.23) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{t_{2}}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{4}})\mathsf{L}\,A+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{2}}A_{ss^{\ast }}+2\text{i}fA_{ss^{\ast }}=\text{ord}(\unicode[STIX]{x1D716}^{4}).\end{eqnarray}$$

On the right-hand side above we indicate that there are higher-order terms that have not been calculated explicitly. One such term is $\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x2202}_{t_{4}}A_{ss^{\ast }}$ , which can be rescued from the garbage on the right and moved to the left-hand side of (2.23) so that

(2.24) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{t_{2}}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{4}})(\mathsf{L}\,A+\unicode[STIX]{x1D716}^{2}A_{ss^{\ast }})+2\text{i}fA_{ss^{\ast }}=\text{ord}(\unicode[STIX]{x1D716}^{4}).\end{eqnarray}$$

Recalling $\unicode[STIX]{x2202}_{s}\unicode[STIX]{x2202}_{s^{\ast }}=\frac{1}{4}\triangle$ , setting $\unicode[STIX]{x1D716}=1$ , throwing out the $\unicode[STIX]{x1D716}^{4}$ garbage and consolidating the time derivatives with

(2.25) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t_{2}}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{t_{4}}\mapsto \unicode[STIX]{x2202}_{t},\end{eqnarray}$$

we thus obtain the $\unicode[STIX]{x1D6F9}=0$ version of (1.4).

2.4 Non-zero $\unicode[STIX]{x1D6F9}$

Now we consider the higher-order terms associated with non-zero $\unicode[STIX]{x1D6F9}$ . TSB show that there are many such terms and most have no clear physical interpretation. Moreover, some of the higher-order terms that do have a clear physical interpretation either do very little, or even degrade the performance of the model e.g. see the discussion of the ‘very weak dispersion regime’ in TSB.

In formulating (1.4) we have taken a minimalist approach by including just two higher-order $\unicode[STIX]{x1D6F9}\neq 0$ terms with clear physical effects

(2.26a,b ) $$\begin{eqnarray}J\left(\unicode[STIX]{x1D6F9},\frac{1}{4}\triangle A\right),\quad \text{and}\quad \frac{\text{i}}{2}\triangle \unicode[STIX]{x1D6F9}\frac{1}{4}\triangle A.\end{eqnarray}$$

To preserve Galilean invariance, and to obtain the correct Doppler shifted dispersion relation in the special case $\unicode[STIX]{x1D6F9}=-Uy$ , with constant  $U$ , one must have the first term in (2.26). The second term in (2.26) also ensures that the YBJ⁺ has the well-known $\unicode[STIX]{x1D701}/2$ -frequency shift (Mooers Reference Mooers1975; Kunze Reference Kunze1985).

We also demand that the YBJ⁺ equation should have an action conservation law. The second term in (2.26) ensures action conservation and other desirable conservation laws: see § 3.

2.5 The reconstituted velocity and phase-averaged kinetic energy

In the discussion of conservation laws in § 3 we will need the backrotated velocity correct to order  $\unicode[STIX]{x1D716}^{4}$ . The field ${\mathcal{U}}_{4}$ is obtained by solving (2.21) after the right-hand side has been simplified by removal of resonant terms. Combining the resulting expression for ${\mathcal{U}}_{4}$ with earlier expressions for ${\mathcal{U}}_{0}$ and  ${\mathcal{U}}_{2}$ , one obtains

(2.27) $$\begin{eqnarray}{\mathcal{U}}=\mathsf{L}\,A\text{e}^{-\text{i}ft_{0}}+\unicode[STIX]{x1D716}^{2}A_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}+\unicode[STIX]{x1D716}^{4}\frac{\text{i}}{f}\unicode[STIX]{x2202}_{t_{2}}A_{s^{\ast }s^{\ast }}^{\ast }\text{e}^{+\text{i}ft_{0}}+\text{ord}(\unicode[STIX]{x1D716}^{6}).\end{eqnarray}$$

The NIW kinetic energy is therefore

(2.28) $$\begin{eqnarray}{\textstyle \frac{1}{2}}|{\mathcal{U}}|^{2}={\textstyle \frac{1}{2}}|\mathsf{L}\,A|^{2}+\unicode[STIX]{x1D716}^{4}{\textstyle \frac{1}{2}}|A_{ss}|^{2}+\cdots\end{eqnarray}$$

where $\cdots$ indicates terms of order $\unicode[STIX]{x1D716}^{6}$ , but also terms of order $\unicode[STIX]{x1D716}^{2}$ and $\unicode[STIX]{x1D716}^{4}$ containing the second harmonics $\text{e}^{\pm 2\text{i}ft}$ . The second harmonics are destroyed by a phase average and

(2.29) $$\begin{eqnarray}|A_{ss}|^{2}=A_{ss}A_{s^{\ast }s^{\ast }}^{\ast }=A_{ss^{\ast }}A_{s^{\ast }s}^{\ast }+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{F}={\textstyle \frac{1}{16}}|\triangle A|^{2}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{F},\end{eqnarray}$$

where $\boldsymbol{F}$ is a horizontal flux whose detailed form is unimportant. Thus phase averaging (2.28) yields

(2.30) $$\begin{eqnarray}\text{phase average}~\left\{{\textstyle \frac{1}{2}}|{\mathcal{U}}|^{2}\right\}={\textstyle \frac{1}{2}}|\mathsf{L}\,A|^{2}+\unicode[STIX]{x1D716}^{4}{\textstyle \frac{1}{2}}{\textstyle \frac{1}{16}}|\triangle A|^{2}+\unicode[STIX]{x1D716}^{4}{\textstyle \frac{1}{2}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{F}+\text{ord}(\unicode[STIX]{x1D716}^{6}).\end{eqnarray}$$

3.1 Conservation laws of YBJ⁺

In this section we discuss the conservation laws of YBJ⁺. To obtain conservation of action, we form $\mathsf{L}^{\!+}A^{\ast }$ (1.4) $+$   $\mathsf{L}^{\!+}A$ (1.4)^* and find

(3.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\frac{1}{2}\left\langle |\mathsf{L}^{\!+}A|^{2}\right\rangle =0,\end{eqnarray}$$

where $\langle \,\rangle$ denotes a spatial average. The action itself can be written as

(3.2) $$\begin{eqnarray}\frac{1}{2}\left\langle |\mathsf{L}^{\!+}A|^{2}\right\rangle =\frac{1}{2}\left\langle |\mathsf{L}\,A|^{2}+\frac{1}{2}\frac{f^{2}}{N^{2}}|\unicode[STIX]{x1D735}A_{z}|^{2}+\frac{1}{16}|\triangle A|^{2}\right\rangle .\end{eqnarray}$$

The first and last terms on the right of (3.2) correspond to the phase-averaged kinetic energy in (2.30) and the middle term on the right of (3.2) is the potential energy, $b^{2}/2N^{2}$ . We conclude that in YBJ⁺, the action is an approximation to the total phase-averaged NIW energy, correct to order $\mathit{Bu}^{2}$ .

Danioux, Vanneste & Bühler (Reference Danioux, Vanneste and Bühler2015) noted that if $\unicode[STIX]{x1D6F9}_{t}=0$ (steady eddies) then the YBJ equation has a second conservation law that can be used to make interesting deductions concerning the concentration of NIW energy into anticyclones. The YBJ⁺ analogue is obtained by forming $\langle \!\text{i}\mathsf{L}^{\!+}A_{t}^{\ast }$ (1.4)  $-$   $\text{i}\mathsf{L}^{\!+}A_{t}$ (1.4)^* $\rangle$ , leading to

(3.3) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\left\langle \frac{1}{4}\frac{f^{2}}{N^{2}}|\unicode[STIX]{x1D735}A_{z}|^{2}+\frac{1}{16}|\triangle A|^{2}\right\rangle =-\left\langle \frac{\text{i}\unicode[STIX]{x1D6F9}}{2f}\unicode[STIX]{x2202}_{t}J(\mathsf{L}^{\!+}A^{\ast },\mathsf{L}^{\!+}A)+\frac{\triangle \unicode[STIX]{x1D6F9}}{4f}\unicode[STIX]{x2202}_{t}|\mathsf{L}^{\!+}A|^{2}\right\rangle .\end{eqnarray}$$

If $\unicode[STIX]{x1D6F9}$ does not depend on $t$ then the right of (3.3) is a time derivative and the linear conservation law of Danioux et al. (Reference Danioux, Vanneste and Bühler2015) follows.

3.2 Xie and Vanneste’s coupled energy

Xie & Vanneste (Reference Xie and Vanneste2015, XV hereafter) derived a generalized-Lagrangian-mean model which couples the YBJ equation to the quasi-geostrophic equation; see also Salmon (Reference Salmon2016) and Wagner & Young (Reference Wagner and Young2016). As an experimental modification of XV, we propose here that eddy evolution is governed by the quasi-geostrophic potential vorticity (PV) equation

(3.4) $$\begin{eqnarray}q_{t}+J(\unicode[STIX]{x1D6F9},q)=0,\end{eqnarray}$$

where the PV of the ‘XV⁺ system’ is

(3.5) $$\begin{eqnarray}q=\triangle \unicode[STIX]{x1D6F9}+\mathsf{L}\unicode[STIX]{x1D6F9}+\frac{\text{i}}{2f}J(\mathsf{L}^{\!+}A^{\ast },\mathsf{L}^{\!+}A)+\frac{1}{4f}\triangle |\mathsf{L}^{\!+}A|^{2}.\end{eqnarray}$$

The XV PV is recovered by $\mathsf{L}^{\!+}\mapsto \mathsf{L}$ in (3.5).

Justification for extension of the $\mathsf{L}\mapsto \mathsf{L}^{\!+}$ rule to XV via an asymptotic expansion is beyond the scope of this paper. Nonetheless, this seemingly ad hoc modification of XV is impressive because the resulting system consisting of (1.4), (3.4) and (3.5) has the expected nonlinear conservation law for ‘coupled energy’: forming  one finds

(3.6) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\left\langle \frac{1}{2}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}|^{2}+\frac{1}{2}\frac{f^{2}}{N^{2}}\unicode[STIX]{x1D6F9}_{z}^{2}\right\rangle =\left\langle \frac{\text{i}\unicode[STIX]{x1D6F9}}{2f}\unicode[STIX]{x2202}_{t}J(\mathsf{L}^{\!+}A^{\ast },\mathsf{L}^{\!+}A)+\frac{\triangle \unicode[STIX]{x1D6F9}}{4f}\unicode[STIX]{x2202}_{t}|\mathsf{L}^{\!+}A|^{2}\right\rangle .\end{eqnarray}$$

Summing (3.3) and (3.6), we obtain the ‘coupled energy’

(3.7) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\left\langle \frac{1}{2}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}|^{2}+\frac{1}{2}\frac{f^{2}}{N^{2}}\unicode[STIX]{x1D6F9}_{z}^{2}+\frac{1}{4}\frac{f^{2}}{N^{2}}|\unicode[STIX]{x1D735}A_{z}|^{2}+\frac{1}{16}|\triangle A|^{2}\right\rangle =0.\end{eqnarray}$$

The first two terms in $\langle \,\rangle$ above constitute the standard quasi-geostrophic energy, the third term is the NIW potential energy and the final term, $|\triangle A|^{2}/16$ , is the result of the $\mathsf{L}\mapsto \mathsf{L}^{\!+}$ modification of XV.

4.1 Formulation

To test the new proposed model, we compare numerical solutions of YBJ and YBJ⁺ to solutions of the non-hydrostatic Boussinesq (BQ) equations, the latter of which we take as ‘truth’. The three models are initialized with a single vertical mode inertial wave,

(4.1a,b ) $$\begin{eqnarray}u=u_{0}\cos mz,\quad \text{and}\quad v=w=b=p=0.\end{eqnarray}$$

Following TSB, this wave is superposed on a barotropic geostrophic flow with streamfunction

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}=\ell U_{0}\sin (x/\ell )\sin (y/\ell );\end{eqnarray}$$

we use the convention $\boldsymbol{U}=(-\unicode[STIX]{x1D6F9}_{y},\unicode[STIX]{x1D6F9}_{x})$ . The pressure is defined so that this flow is a steady solution of the Boussinesq equations (but note that the fully nonlinear equation set is nevertheless integrated). We verified that the barotropic component departs from its steady solution by much less than one part in $10^{8}$ over the time of interest. Moreover, this very weak evolution of the geostrophic component was found to be comparable over the entire range of parameters considered, thus eliminating it as a source of error. The wave component solution of the Boussinesq model is obtained by subtracting the steady solution, equation (4.2), from the total solution. Because we use a very small wave-to-flow amplitude factor, $u_{0}/U_{0}=10^{-4}$ , this is a good approximation.

For simplicity, $N$ and $f$ are assumed to be constant, the domain is horizontally periodic $2\unicode[STIX]{x03C0}\ell \times 2\unicode[STIX]{x03C0}\ell$ , and is bounded vertically by rigid lids located at $z=0$ and $-\unicode[STIX]{x03C0}h$ . The dimensional wavenumbers are therefore quantized with $k=k^{\prime }/\ell$ and $m=m^{\prime }/h$ where $k^{\prime }$ and $m^{\prime }$ are integers. The wave in (4.1) has $k^{\prime }=0$ and variable  $m^{\prime }$ . Notice that the total horizontal wavenumber of the geostrophic flow in (4.2) is $\sqrt{2}/\ell$ . Because the geostrophic flow is barotropic, the vertical wavenumber $m$ of the NIWs is set by the initial condition in (4.1) and does not change throughout the integration.

Important dimensionless numbers emerging from this formulation are the domain-scale Burger number of TSB

(4.3) $$\begin{eqnarray}\mathit{Bu}^{1,1}\stackrel{\mathit{def}}{=}\left(\frac{Nh}{f\ell }\right)^{2},\end{eqnarray}$$

and the Rossby number

(4.4) $$\begin{eqnarray}\mathit{Ro}\stackrel{\mathit{def}}{=}\frac{U_{0}}{f\ell }.\end{eqnarray}$$

(The superscript $1,1$ in (4.3) indicates that $k^{\prime }=1$ and $m^{\prime }=1$ .) There is a third non-dimensional number, $u_{0}/U_{0}$ , which we set to $10^{-4}$ ; this is so small that wave–wave interactions and wave feedback on the geostrophic flow can be ignored over the duration of the integration (typically 50 inertial periods).

The main interaction is distortion of the wave field in (4.1) by the barotropic geostrophic flow in (4.2). The initial Burger number of the wave field in (4.1) is zero because the flow has infinite horizontal scale, $k^{\prime }=0$ . The geostrophic flow quickly impresses its horizontal scale $\ell$ , onto the waves. Subsequent evolution results in the wave field developing a spectrum of horizontal wavenumbers, and depending on  $m$ , the high $k$ components of this spectrum might have large values of the effective wave Burger number, $\mathit{Bu}$ in (1.1). We comment further on this important point in § 4.5.

Table 1 outlines the parameter set explored. In a nutshell, simulations span non-dimensional vertical wavenumbers $m^{\prime }=1$ to 12, and three Rossby numbers. To fix ideas, the bottom entries of the table give one possible set of realistic dimensional parameters corresponding to the above dimensionless parameters. Of course, any combination of parameters yielding a given dimensionless numbers set leads to the same YBJ and YBJ⁺ solution. This is made clear in the non-dimensional form of the YBJ⁺ equation. Using $\ell$ to scale $x$ and  $y$ , $h$ to scale  $z$ and $(\mathit{Bu}^{1,1}f)^{-1}$ to scale  $t$ , the non-dimensional YBJ⁺ equation is

(4.5) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\mathsf{L}^{\!+}A+\frac{\mathit{Ro}}{\mathit{Bu}^{1,1}}\left[J(\unicode[STIX]{x1D6F9},\mathsf{L}^{\!+}A)+\frac{\text{i}}{2}\triangle \unicode[STIX]{x1D6F9}\mathsf{L}^{\!+}A\right]+\frac{\text{i}}{2}\triangle A=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D6F9}=\sin x\sin y$ and, since $N$ is constant,

(4.6) $$\begin{eqnarray}\mathsf{L}^{\!+}A=A_{zz}+{\textstyle \frac{1}{4}}\mathit{Bu}^{1,1}\triangle A.\end{eqnarray}$$

Once again, the plain YBJ version is obtained by removing the term ${\textstyle \frac{1}{4}}\mathit{Bu}^{1,1}\triangle A$ from the definition of $\mathsf{L}^{\!+}A$ . The above equation set highlights an important distinction between the two approximate models: YBJ only depends on the ratio $\mathit{Bu}^{1,1}/\mathit{Ro}$ , whereas YBJ⁺ depends on $\mathit{Bu}^{1,1}$ and $\mathit{Ro}$ independently. This signals a superior parameter sensitivity for YBJ⁺. We comment on this in § 4.6.

Table 1.

Key parameters of the numerical simulations.

The testing procedure goes as follows: the YBJ and YBJ⁺ versions of (4.5) and (4.6) are integrated over 50 inertial periods using the parameter set in table 1. The accuracy of the two models is then assessed by comparing their solutions to that of the non-hydrostatic Boussinesq model used by Asselin, Bartello & Straub (Reference Asselin, Bartello and Straub2018). To put emphasis on the models rather than numerics, all three equation sets are integrated using the same numerical methods and parameters (cf. table 1). That is, the YBJ, YBJ⁺ and BQ models are pseudo-spectral, allowing horizontal derivatives to be computed with spectral accuracy. The $2/3$ rule is used to remove aliased modes. Vertical derivatives are approximated with second-order centred finite differences. Time integration is accomplished with the leap-frog scheme with a diffusion parameter $\unicode[STIX]{x1D6FE}=10^{-2}$ (Asselin Reference Asselin1972). To facilitate comparison between the models, no spatial diffusion is applied (except in a supporting simulation in § 4.5). Both temporal and vertical resolution were increased to convergence to eliminate them as sources of error. Horizontal resolution is much lower than vertical resolution to allow explicit treatment of the dispersive term for both YBJ and YBJ⁺. Indeed, with a horizontal resolution of $64\times 64$ , explicit treatment of the dispersive term in plain YBJ requires a time step of only 20 s. We verified that results below are unchanged as horizontal resolution is increased.

4.2 Overview comparison of YBJ and YBJ⁺

In both YBJ and YBJ⁺ the wave kinetic energy, estimated using the leading-order term in the expansion (2.28), is

(4.7) $$\begin{eqnarray}\text{WKE}\stackrel{\mathit{def}}{=}{\textstyle \frac{1}{2}}|\mathsf{L}\,A|^{2}.\end{eqnarray}$$

To assess the error relative to the Boussinesq solution we introduce

(4.8) $$\begin{eqnarray}E(t)\stackrel{\mathit{def}}{=}\sqrt{\frac{\langle (\text{WKE}-\text{WKE}_{BQ})^{2}\rangle }{\langle \text{WKE}_{BQ}^{2}\rangle }}.\end{eqnarray}$$

Above, the brackets $\langle \,\rangle$ indicate spatial averaging over the horizontal coordinates and subscript BQ indicates the solution of the Boussinesq solutions. The error is independent of the vertical level because WKE and $\text{WKE}_{BQ}$ share the same vertical structure.

(We emphasize that to define WKE for YBJ⁺ we use $\mathsf{L}$ , not $\mathsf{L}^{\!+}$ , in (4.7). We have verified by examination of the numerical solutions that the difference between $\mathsf{L}\,A$ and $\mathsf{L}^{\!+}A$ is small and thus there are only small differences between the two definitions of WKE. Moreover, the discussion following (3.2) provides a physical rationale in support of (4.7): $|\mathsf{L}^{\!+}A|^{2}$ contains a term corresponding to NIW potential energy.)

Figure 2 shows $\bar{E}$ , the $E(t)$ time averaged over the first 50 inertial periods of integration. The error is plotted against $m^{\prime }$ for the three Rossby numbers (one for each panel). In virtually all cases YBJ⁺ is more accurate that YBJ, as one would expect from the superior accuracy of YBJ⁺. Furthermore, for fixed $m^{\prime }$ the error decreases monotonically with $\mathit{Ro}$ for both models. In the lowest $\mathit{Ro}$ case (bottom panel), both models are accurate with just a few per cent error. This is consistent with the improved accuracy of their formal expansion as $\mathit{Ro}\rightarrow 0$ .

Figure 2.

Error in WKE, defined in (4.8), and averaged over the first 50 inertial periods. Each marker represents a different simulation, with blue dots and green pluses for the error in the YBJ and YBJ⁺ models, respectively. Each panel has a fixed Rossby number.

The error dependence on $m^{\prime }$ is more puzzling. With increasing $m^{\prime }$ in the initial condition (4.1), the wave Burger number $\mathit{Bu}$ in (1.1) decreases as $m^{\prime -2}$ . The formal expansion indicates that both YBJ and YBJ⁺ become more accurate as $\mathit{Bu}$ is reduced. One might thus expect that the error in figure 2 should decrease monotonically and rapidly with increasing  $m^{\prime }$ . In fact the dependence of error on $m^{\prime }$ in figure 2 is non-monotonic. Moreover, while the error does eventually decrease monotonically for large  $m^{\prime }$ , this decrease is quite slow. To understand, or at least rationalize, the non-intuitive dependence of error on  $m^{\prime }$ , we examine two illustrative cases: one with low ( $m^{\prime }=1$ ) and one with high ( $m^{\prime }=8$ ) vertical wavenumber.

4.3 A detailed look at the case $m^{\prime }=1$ , $\mathit{Ro}=0.05$ and $\mathit{Bu}^{1,1}=1/4$

We first examine the case that one would naively assume to be most challenging: the solution with $\mathit{Ro}=0.05$ , $\mathit{Bu}^{1,1}=1/4$ and $m^{\prime }=1$ . In this strongly dispersive case, it is useful to introduce the parameter

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D702}\stackrel{\mathit{def}}{=}m^{\prime 2}\frac{\mathit{Ro}}{\mathit{Bu}^{1,1}}.\end{eqnarray}$$

In the limit $\unicode[STIX]{x1D702}\rightarrow 0$ , the YBJ⁺ equation (4.5) can be solved perturbatively with the initial condition in (4.1) (Young & Ben Jelloul Reference Young and Ben Jelloul1997; Thomas et al. Reference Thomas, Smith and Bühler2017). For YBJ⁺ this solution is

(4.10) $$\begin{eqnarray}A=\text{e}^{\text{i}\unicode[STIX]{x1D714}t}[1+\unicode[STIX]{x1D702}(1-\text{e}^{-\text{i}\unicode[STIX]{x1D6FA}t})\sin x\sin y+\text{ord}(\unicode[STIX]{x1D702}^{2})]\cos (m^{\prime }z),\end{eqnarray}$$

where

(4.11a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}\stackrel{\mathit{def}}{=}\frac{2}{2m^{\prime 2}+\mathit{Bu}^{1,1}},\quad \text{and}\quad \unicode[STIX]{x1D714}\stackrel{\mathit{def}}{=}\frac{2m^{\prime 2}+\mathit{Bu}^{1,1}}{8}\left(\frac{\mathit{Ro}}{\mathit{Bu}^{1,1}}\right)^{2}.\end{eqnarray}$$

(In their ‘strong dispersion regime’, TSB considered an asymptotic limit of the rotating shallow-water equations in which the Burger number is fixed and of order unity, while the Rossby number limits to zero. The more exact TSB version of (4.10) has the same structure: the leading-order term is a spatially uniform near-inertial oscillation and the higher-order terms contain the spatial dependence of the solution.)

The wave kinetic energy of the $\unicode[STIX]{x1D702}\ll 1$ solution in (4.10) is

(4.12) $$\begin{eqnarray}\text{WKE}(t)=\text{WKE}(0)[1+2\unicode[STIX]{x1D702}(1-\cos \unicode[STIX]{x1D6FA}t)\sin x\sin y+O(\unicode[STIX]{x1D702}^{2})].\end{eqnarray}$$

The YBJ version of (4.12) is obtained by neglecting $\mathit{Bu}^{1,1}$ relative to $2m^{\prime 2}$ in the definition of $\unicode[STIX]{x1D6FA}$ above. In the present case $\unicode[STIX]{x1D702}=0.2$ and (4.12) predicts $\text{WKE}(t)$ has modulations with amplitude reaching 80 % of the initial level $\text{WKE}(0)$ . This modulation has the same horizontal structure as the streamfunction (or vorticity, since in this case $\triangle \unicode[STIX]{x1D6F9}=-2\unicode[STIX]{x1D6F9}$ ). Indeed, the effect of the refractive term $\triangle \unicode[STIX]{x1D6F9}\mathsf{L}^{\!+}A$ in (4.5) is to imprint the length scale of $\triangle \unicode[STIX]{x1D6F9}$ onto the wave field. Now the squared horizontal wavenumber of $\unicode[STIX]{x1D6F9}$ is 2 and therefore, although the domain-scale Burger number is $\mathit{Bu}^{1,1}={\textstyle \frac{1}{4}}$ , the effective Burger number of the NIW immediately after this initial imprinting is

(4.13) $$\begin{eqnarray}\mathit{Bu}_{init}={\textstyle \frac{1}{2}}.\end{eqnarray}$$

The formal expansion used to derive both YBJ and YBJ⁺ assumes that $\mathit{Bu}\ll 1$ and thus this test with an initial Burger number $1/2$ should be challenging.

Figure 3.

Case: $m^{\prime }=1$ , $\mathit{Ro}=0.05$ and $\mathit{Bu}^{1,1}=1/4$ . Distribution of normalized wave kinetic energy as predicted by the YBJ, YBJ⁺ and Boussinesq (BQ) models at intervals of one inertial period (IP).

Figure 3 displays the horizontal distribution of wave kinetic energy in the first ten inertial periods of integration. Each row is associated with one of the three models, with the middle row BQ regarded as truth. These numerical solutions are broadly consistent with the asymptotic solution, (4.12). First, wave energy concentrates in regions of negative vorticity, i.e. following the predicted spatial structure $\sin x\sin y$ in (4.12). YBJ predicts that the period of the oscillation, $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FA}$ , is 4 inertial periods (IP); YBJ⁺predicts 4.5 IP. Figure 3 shows that the YBJ⁺ prediction is closer to the period of the Boussinesq solution. Finally, (4.12) predicts a modulation reaching $4\unicode[STIX]{x1D702}=80\,\%$ of $\text{WKE}(0)$ , which is comparable to that of the BQ solution. We verified that both YBJ and YBJ⁺ become more accurate as $\mathit{Ro}\rightarrow 0$ (and thus $\unicode[STIX]{x1D702}\rightarrow 0$ ).

Figure 4.

Case: $m^{\prime }=1$ , $\mathit{Ro}=0.05$ and $\mathit{Bu}^{1,1}=1/4$ . (a) Time series of the error, (4.8), for both YBJ (blue) and YBJ⁺ (green) models. (b) Time series of the wave kinetic energy density at the core of the bottom left quarter for YBJ (blue), YBJ⁺ (green) and Boussinesq (BQ, red) solutions.

Although the overall horizontal structure of the BQ solution is well captured by the YBJ and YBJ⁺ models, the amplitude of the YBJ/YBJ⁺ oscillation is roughly 50 % larger than that of the BQ solution. This can be appreciated in the bottom panel of figure 4, which displays the time series of WKE at one of its maxima, the centre of the bottom left quarter, $x^{\prime }=y^{\prime }=\unicode[STIX]{x03C0}/2$ . YBJ and YBJ⁺ equally overestimate the amplitude of the BQ oscillation and miss some of its detailed structure. The plain YBJ solution also quickly shifts out of phase with the BQ solution. By comparison, YBJ⁺ predicts the phase of the oscillation significantly better than YBJ, as expected from its more accurate dispersion relation. This is reflected by an overall smaller error for YBJ⁺ during most of the integration (cf. figure 4 a), until the YBJ solution coincidentally comes back into phase after 40 inertial periods.

To recap, in this challenging high- $\mathit{Ro}~m^{\prime }=1$ case, WKE undergoes a weak modulation whose horizontal structure is well captured by both YBJ and YBJ⁺ models. The main source of error comes from an amplitude mismatch. YBJ⁺ captures the phase of oscillation much better than does YBJ, and this accounts for its higher accuracy.

4.4 A detailed look at case $m^{\prime }=8$ , $\mathit{Ro}=0.05$ and $\mathit{Bu}^{1,1}=1/4$

Let us now turn our attention to a case with a larger vertical wavenumber, $m^{\prime }=8$ , keeping the Rossby number to its most challenging value, $\mathit{Ro}=0.05$ . Relative to the $m^{\prime }=1$ solution of figure 4, $m^{\prime }=8$ reduces the Burger number by a factor of 64: the analogue of (4.13) is an initially imprinted Burger number

(4.14) $$\begin{eqnarray}\mathit{Bu}_{init}={\textstyle \frac{1}{128}}.\end{eqnarray}$$

One expects that with this small Burger number, YBJ and YBJ⁺ should be more faithful to the BQ solution. Alas, the error summarized in figure 2 is not much less at $m^{\prime }=8$ than at $m^{\prime }=1$ .

Figure 5.

Case: $m^{\prime }=8$ , $\mathit{Ro}=0.05$ and $\mathit{Bu}^{1,1}=1/4$ . (a) Normalized WKE distribution for the first 50 inertial periods of integration. A logarithmic scale is used to reveal both the fine-scale structure of the solution and the strong focussing of energy in the centre of negative vorticity cores. (b) Corresponding time series of error, equation (4.8).

To understand why, we refer to figure 5, which displays the horizontal WKE distribution (a) along with the associated error time series (b) for $m^{\prime }=8$ . As in the $m^{\prime }=1$ case, WKE initially concentrates in regions of negative vorticity. But there are two striking differences: (i) the anomaly is almost two orders of magnitude stronger: up to 3000 % as opposed to about 60 % in the $m^{\prime }=1$ case (note that a logarithmic scale is used for $m^{\prime }=8$ ); (ii) the $m^{\prime }=8$ solution has a much finer horizontal structure. Because $\mathit{Bu}\propto k^{2}$ , a fine-scale structure implies that the effective $\mathit{Bu}$ is, in fact, much larger than the initially imprinted Burger number in (4.14). This crucial result deserves further discussion.

4.5 Evolution of effective Burger number

As shown above, a larger $m^{\prime }$ does not imply that $\mathit{Bu}$ remains small throughout the integration. As $m^{\prime }$ is increased, $\mathit{Bu}_{init}$ tends to zero and the dispersive term, $\text{i}\triangle A/2$ in (4.5), is initially very small relative to advection and refraction. Weakening dispersion results in the NIWs developing finer scale features, with higher horizontal wavenumber  $k$ , and a larger effective $\mathit{Bu}\propto k^{2}$ . This rationalizes the slow decrease of error as $m^{\prime }$ is increased, but it also begs the question: does weakening dispersion fatally compromise the YBJ and YBJ⁺ models, which are based on the assumption of small $\mathit{Bu}$ ? At least in our test computations, the answer is ‘no’ because dispersive wave propagation is effective in opposing advection and refraction.

Rocha, Wagner & Young (Reference Rocha, Wagner and Young2018) show that lateral straining increases the horizontal group velocity of NIWs so that the waves accelerate and escape from straining regions. This ‘wave escape’ ensures that NIWs do not suffer a forward cascade of action to the small horizontal scales on which dissipation is effective. In a flow that rapidly stirs passive scalar into diffusive oblivion, NIWs concentrate in vortex centres and wave action is robustly conserved, even as passive-scalar variance disappears in a few eddy turnover times.

Figure 6.

Cumulative wave energy as a function of the Burger number. The different colours represent YBJ⁺ solutions for different  $m^{\prime }$ , all evaluated after 50 inertial periods. Horizontal wavenumbers are binned by integer steps, and therefore the spacing between points in $\mathit{Bu}$ thus depends on  $m^{\prime }$ .

As evidence that wave escape is operating in our solutions, figure 6 shows the cumulative distribution of wave energy, ${\textstyle \frac{1}{2}}|\mathsf{L}^{\!+}A|^{2}$ , as a function of the Burger number. That is, the ordinate at a given $\mathit{Bu}$ corresponds to the sum of wave energy at all Burger numbers equal to or smaller than $\mathit{Bu}$ . Remarkably, the case $m^{\prime }=1$ actually contains less energy at $\mathit{Bu}\sim 1$ than all other cases, including $m^{\prime }=8$ . Indeed, in the former case more of the energy remains in the $k^{\prime }=0$ mode, or equivalently, oscillations of the solution have a much weaker amplitude. In that sense, the effective $\mathit{Bu}$ for $m^{\prime }=8$ is more challenging than for $m^{\prime }=1$ , even though $\mathit{Bu}_{init}$ at $m^{\prime }=8$ is 64 smaller than at $m^{\prime }=1$ .

The numerical size of $\mathit{Bu}$ is important in this discussion. Suppose that the initial $\mathit{Bu}_{init}$ is very small, say $1/128$ . If advection increases the effective $\mathit{Bu}$ to $1/5$ , then that is a large increase relative to $\mathit{Bu}_{init}=1/128$ . But waves with $Bu=1/5$ have $\unicode[STIX]{x1D714}=1.09\,f$ and are therefore still near inertial. Thus it is striking that in figure 6 the range $0\leqslant \mathit{Bu}\leqslant 0.2$ contains at least 90 % of wave energy for all $m^{\prime }$ shown. This also means that 90 % of waves oscillate with a frequency between $f$ and $1.09\,f$ . The amount of energy above $\mathit{Bu}=1$ (equivalently above $1.4\,f$ ) is minute, reaching a maximum of approximately 2 % in the case $m^{\prime }=8$ . Thus although small horizontal scales are created by advection and refraction, and although there is a very large increase in $\mathit{Bu}$ relative to $\mathit{Bu}_{init}$ , the amount of energy reaching Burger numbers large enough to challenge YBJ⁺, such as $\mathit{Bu}=1$ in figure 1, is small.

We emphasize that although dispersion becomes weak at high $m^{\prime }$ , it is by no means negligible. Figure 7 illustrates this point by showing equilibrium wave energy spectra for YBJ⁺ at $m^{\prime }=8$ and $\mathit{Ro}=0.05$ , with and without the dispersive term $\text{i}\triangle A/2$ in (4.5). For these simulations only, the resolution was increased to $256^{3}$ and a small amount of hyperdiffusion was included to smooth small horizontal scales. The cases with and without dispersion are strikingly different. In the absence of dispersion (green), wave energy tends to equipartition between all the available modes, approaching a $k^{1}$ spectrum till an abrupt hyperdiffusive cutoff. In other words, without the small dispersive term, large scales are rapidly purged of their energy, which cascades to smaller scales. But with small dispersive case (purple) energy remains at larger scales with small $\mathit{Bu}$ . After a few hundred inertial periods, more than 90 % of the initial energy remains, whereas about only 10 % remains in the case without dispersion. (We have verified that these results are robust to changes in resolution and the strength of the hyperdiffusion coefficient.)

Figure 7.

Wave energy spectra as a function horizontal wavenumber. YBJ⁺ solutions with $Ro=0.05$ and $m^{\prime }=8$ are shown with (purple) and without (olive) dispersion. Spectra are captured at equilibria after 300 inertial periods, and are averaged over the 50 surrounding inertial periods. The spectra are normalized so that the initial condition sums to unity. These simulations have the same parameters as in the rest of this paper, except that the resolution is $256^{3}$ and hyperdiffusion was added to remove small horizontal scales.

4.6 Dimensional reduction of the parameter space

The non-dimensional equation set (4.5) and (4.6) also deserves further discussion. For a fixed initial condition, the solution of a given equation system, YBJ, YBJ⁺ or BQ, is uniquely defined by its set of dimensionless numbers. For instance, the solution of the inviscid non-hydrostatic Boussinesq equations is determined by three control parameters: $\mathit{Ro}$ , $\mathit{Bu}$ and $Ar$ , where $Ar$ is the aspect ratio. That is, for a fixed initial condition, there is only one solution associated with each point of the cube of dimensions $\mathit{Ro}\times \mathit{Bu}\times Ar$ . Both the YBJ and YBJ⁺ model assume hydrostatic balance, which imply $Ar\rightarrow 0$ : the space of possible solutions shrinks to a $\mathit{Ro}$ – $\mathit{Bu}$ plane. Furthermore, (4.5) and (4.6) highlight that YBJ in fact only depends on Burger-to-Rossby ratio, and thus admits only one solution per straight line in the $\mathit{Bu}$ – $\mathit{Ro}$ plane. This was appreciated by Thomas et al. (Reference Thomas, Smith and Bühler2017), who defined their various dispersion regimes based on the slope of the line crossing the $\mathit{Bu}$ – $\mathit{Ro}$ plane (cf. their figure 3).

Figure 8.

Normalized wave kinetic energy density at the surface, $m^{\prime }=6$ , for the three models (three rows). Each column displays a solution with a different set of $\mathit{Ro}$ and $\mathit{Bu}^{1,1}$ , but with their ratio held fixed. As such, when time is appropriately stretched the YBJ model has the same output (upper panels). On the other hand, the output of the Boussinesq and YBJ⁺ models are influenced by the independent changes in $\mathit{Ro}$ and $\mathit{Bu}^{1,1}$ .

YBJ⁺ has more degrees of freedom than YBJ: YBJ⁺ is sensitive to $\mathit{Bu}$ and $\mathit{Ro}$ independently. This is illustrated in figure 8. As usual, the three rows are associated with the three models. But instead of showing a time series from left to right, figure 8 shows snapshots from two simulations with different Rossby and Burger numbers but fixed $\mathit{Bu}^{1,1}/\mathit{Ro}$ . Here, the right panel has its Burger and Rossby numbers multiplied by 3, and is thus shown at a dimensional time divided by 3. As expected, the YBJ solution is the same in both cases. However, both BQ and YBJ⁺ change rather importantly between the two cases, leaving YBJ behind. Unsurprisingly, YBJ⁺ does better as $\mathit{Bu}$ and $\mathit{Ro}$ are increased while keeping their ratio fixed.