2.1. Non-local form of the LANS- $\alpha$ equations

The LANS- $\alpha$ equations with the Boussinesq approximation (the Boussinesq- $\alpha$ equations) are given by (Holm Reference Holm1999; Holm & Titi Reference Holm and Titi2005)

(2.1a) \begin{align} \frac {\partial \boldsymbol{v}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^T + f\hat {\boldsymbol{z}}\times \boldsymbol{u} + \boldsymbol{\nabla }\phi +\frac {\rho }{\rho _0}g\hat {\boldsymbol{z}} &= \mathscr{D}_{\nu }\boldsymbol{v} + \mathcal{F}, \end{align}

(2.1b) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho -\boldsymbol{u}\boldsymbol{\cdot }b\hat {\boldsymbol{z}} &=\mathscr{D}_\kappa \rho , \end{align}

(2.1c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0, \end{align}

(2.1d) \begin{align} \boldsymbol{v} &= \boldsymbol{u} - \alpha ^2 \Delta \boldsymbol{u}, \end{align}

where $\boldsymbol{v}=(v_1,v_2,v_3)$ is the three-dimensional velocity, $f$ is the Coriolis parameter, $\phi$ is a modified pressure, $\rho$ is the buoyancy variable (density), $\rho _0$ is a reference density, $g$ is the gravitational constant, $b$ is the strength of the underlying density stratification, $\alpha \in \mathbb{R}$ is a parameter with units of length, $\mathcal{F}$ is a forcing term and $\mathscr{D}_\nu$ and $\mathscr{D}_\kappa$ are diffusive operators. For the theoretical developments of this work, we will set these to the usual $\mathscr{D}_\nu \boldsymbol{v} =\nu \Delta \boldsymbol{v}$ and $\mathscr{D}_\kappa \rho = \kappa \Delta \rho$ , where $\nu$ is the viscosity coefficient and $\kappa$ is buoyancy diffusivity, but will use hyperviscosity in the simulations done in § 4, as this is a standard method to increase the inertial range (see e.g. Smith & Waleffe Reference Smith and Waleffe2002).

While the length scale $\alpha$ can be interpreted in many ways – such as the typical deviation of a particle trajectory from its time-averaged path (Holm et al. Reference Holm, Jeffery, Kurien, Livescu, Taylor and Wingate2005) or a filter width (Lunasin et al. Reference Lunasin, Kurien, Taylor and Titi2007) – we adopt the interpretation that it is the length scale below which the dynamics is regularised. This regularisation is accomplished through the introduction of a new velocity field, $\boldsymbol{u}=(u_1, u_2, u_3)$ , as defined in (2.1d ) (note that $\alpha ^2\geq 0$ ). Since the advecting velocity is regularised, the development of sharp gradients at scales smaller than $\alpha$ is mitigated. The field $\boldsymbol{v}$ is referred to as the associated circulation velocity, i.e. the quantity whose line integral appears in the Kelvin circulation theorem along loops transported by $\boldsymbol{u}$ (see e.g. Holm & Titi Reference Holm and Titi2005). The relationship between $\boldsymbol{u}$ and $\boldsymbol{v}$ in (2.1d ) can be otherwise interpreted within the Lagrangian-averaging closure; we summarise this interpretation below.

The difference between the Boussinesq- $\alpha$ equations and the usual Boussinesq equations is essentially due to the use of $\boldsymbol{u}$ as the transport velocity. Denoting the Helmholtz operator as

(2.2) \begin{equation} \mathcal{S}\overset {{\text{def}}}{=} \big(1-\alpha ^2\varDelta \big), \end{equation}

$\boldsymbol{u}$ is constructed via $\mathcal{S}^{-1}\boldsymbol{v}$ (2.1d ). This results in an additional advective term, $\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^T$ (also written as $(\boldsymbol{\nabla }\boldsymbol{u})^T \boldsymbol{v}$ or $v_j\partial _i u_j$ ), in (2.1a ) and necessitates a modification to the pressure. The usual pressure $p$ has been replaced with the modified pressure $\phi$ in (2.1a ), where

(2.3) \begin{equation} \phi \overset {{\text{def}}}{=} \frac {p}{\rho _0} - \frac {1}{2}|\boldsymbol{u}|^2-\frac {\alpha ^2}{2}|\boldsymbol{\nabla }\boldsymbol{u}|^2. \end{equation}

From the definition of $\mathcal{S}$ (2.2), the incompressibility of $\boldsymbol{u}$ implies that of $\boldsymbol{v}$ , and vice versa, since the order of the divergence and the Laplacian can be interchanged. Note also that $\mathcal{S}$ can be applied to scalar functions, with the interpretation that $\varDelta$ is the scalar Laplacian rather than the vector Laplacian. Using the Helmholtz operator to modify the advecting velocity does not introduce dissipative effects, since (2.1a ) and (2.1b ) are time-reversible when $\nu =\kappa =0$ .

In (2.1a – d ), the density has been decomposed as

(2.4) \begin{equation} \rho _{\textit{full}} = \overline {\rho }+\rho , \quad \overline {\rho } = \rho _0 - bz, \quad b\gt 0. \end{equation}

The domain is $[0,L]^3$ , for a length scale $L$ , with triply periodic boundary conditions, to facilitate the comparison of the results presented here with those of previous turbulence studies and the classical fast singular QG limit (Embid & Majda Reference Embid and Majda1998).

We non-dimensionalise the Boussinesq- $\alpha$ equations by considering a characteristic length scale $L$ (e.g. the domain length) and velocity scale $U$ , which correspond to the advective time scale $L/U$ . Using the Brunt–Väisälä frequency $N=(gb/\rho _0)^{1/2}$ , we define the buoyancy fluctuation scale $B=NU/g$ in order to non-dimensionalise the density fluctuations. Note that throughout this work, $N$ is considered constant. We non-dimensionalise $\boldsymbol{x}'=\boldsymbol{x}/L$ , $t'=t/T$ , $\boldsymbol{v'}=\boldsymbol{v}/U$ , $\rho '=\rho /(\rho _0 B)$ and the pressure via $p'=p/\overline {p}$ , where $\overline {p}$ is a reference pressure (cf. Majda Reference Majda2003). Using these scalings, the Rossby, Froude, Reynolds, Prandtl and Euler numbers arise naturally:

(2.5) \begin{equation} {\textit{Ro}} \overset {{\text{def}}}{=} \frac {U}{\textit{fL}}, \quad {\textit{Fr}} \overset {{\text{def}}}{=} \frac {U}{NL}, \quad Re \overset {{\text{def}}}{=} \frac {\rho _0UL}{\nu }, \quad Pr \overset {{\text{def}}}{=} \frac {\nu }{\rho _0\kappa }, \quad Eu \overset {{\text{def}}}{=} \frac {\overline {p}}{\rho _0 U^2}. \end{equation}

The parameter $F$ , which measures the relative strength of rotation versus stratification,

(2.6) \begin{equation} F\overset {{\text{def}}}{=} \frac {\textit{Fr}}{\textit{Ro}}=\frac {f}{N}, \end{equation}

also arises in the system and must be assumed fixed when taking the QG limit. This parameter is closely related to the Burger number

(2.7) \begin{equation} Bu\overset {{\text{def}}}{=} \frac {Ro^2}{{\textit{Fr}}^2}=\frac {1}{F^2}, \end{equation}

and plays an important role in the dynamics.

Omitting the primes, the non-dimensionalised version of (2.1) is given by

(2.8a) \begin{align} \frac {\partial \boldsymbol{v}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} +\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^T+ \frac {1}{\textit{Ro}} \hat {\boldsymbol{z}} \times \boldsymbol{u} + \boldsymbol{\nabla }\phi + \frac {1}{\textit{Fr}}\rho \hat {\boldsymbol{z}} &= \frac {1}{\textit{Re}}\Delta \boldsymbol{v}+\mathcal{F}, \end{align}

(2.8b) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho - \frac {1}{\textit{Fr}}\boldsymbol{u}\boldsymbol{\cdot }\hat {\boldsymbol{z}} &= \frac {1}{\textit{RePr}}\Delta \rho , \end{align}

(2.8c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0, \end{align}

(2.8d) \begin{align} \boldsymbol{u}&=\mathcal{S}^{-1}\boldsymbol{v}, \end{align}

(2.8e) \begin{align} \boldsymbol{\nabla }\varDelta ^{-1}\left (\frac {1}{\textit{Ro}}\hat {\boldsymbol{z}}\boldsymbol{\cdot }\mathcal{S}^{-1} \boldsymbol{\omega } -\boldsymbol{\nabla }\boldsymbol{\cdot }(\mathcal{S}^{-1}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }(\mathcal{S}^{-1}\boldsymbol{v})^T)-\frac {1}{\textit{Fr}}\frac {\partial \rho }{\partial z}\right )&=\boldsymbol{\nabla }\phi . \end{align}

Here, the vorticity $\boldsymbol{\omega }\overset {{\text{def}}}{=} \boldsymbol{\nabla }\times \boldsymbol{v}$ and $\mathcal{S}$ is actually the non-dimensionalised Helmholtz operator, $(1-L_\alpha \Delta )$ , where $L_\alpha \overset {{\text{def}}}{=} \alpha ^2/L^2$ . Since $\alpha$ is a fraction of the domain length, and we use the domain size as the characteristic length scale, $L_\alpha \in [0,1]$ . The modified pressure $\phi$ and forcing $\mathcal{F}$ are also the non-dimensional versions. The scalings $BgLU^{-2} = (Fr)^{-1}$ and $Eu=(Ro)^{-1}$ have been used in (2.8) (Majda Reference Majda2003, § 7.4). The non-local form of the momentum equations is found by solving the elliptic equation for the modified pressure, i.e. replacing $\boldsymbol{\nabla }\phi$ in (2.8a ) with (2.8e ). The modified pressure thus does not further differentiate the Boussinesq- $\alpha$ equations from the unregularised Boussinesq equations; $\phi$ is a free variable in exactly the same sense that $p$ is.

This system reduces to the usual Boussinesq equations in the case of $\alpha =0$ , so we expect the slow and fast dynamics derived for the LANS- $\alpha$ model to reduce to those found in Embid & Majda (Reference Embid and Majda1998), Babin et al. (Reference Babin, Mahalov, Nicolaenko and Zhou1997b ) and Wingate et al. (Reference Wingate, Embid, Holmes-Cerfon and Taylor2011) when $\alpha =0$ .

The LANS- $\alpha$ equations may be viewed as a Lagrangian-averaged closure formulated within the Euler–Poincaré variational framework (Holm, Marsden & Ratiu Reference Holm, Marsden and Ratiu1998; Holm Reference Holm1999). In common with generalised Lagrangian mean (GLM) theory (Andrews & McIntyre Reference Andrews and McIntyre1978), the circulation velocity $\boldsymbol{v}$ is advected by the transport velocity $\boldsymbol{u}$ , which is identified with the Lagrangian-mean velocity. Closure of the Lagrangian-averaged system requires independently prescribing the covariance of particle displacement fluctuations $\langle \boldsymbol{\xi }\boldsymbol{\xi }\rangle$ , where $\boldsymbol{\xi }$ denotes a rapid fluctuation of a fluid parcel trajectory; this covariance enters through a general operator $\tilde {\Delta } = \boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\xi }\boldsymbol{\xi }\rangle \boldsymbol{\cdot }\boldsymbol{\nabla}$ that governs the relationship between $\boldsymbol{u}$ and $\boldsymbol{v}$ in the full Lagrangian-averaged system (Holm Reference Holm1999). The Helmholtz operator $\mathcal{S}$ is an isotropic approximation of this operator. Specifically, the isotropy condition $\langle \xi ^k \xi ^l\rangle = \alpha ^2 \delta ^{kl}$ , where $\alpha$ is constant if the statistics are homogeneous, reduces $\tilde {\Delta }$ to $\alpha ^2 \Delta$ and yields (2.1d ). In GLM theory, the standard interpretation identifies the circulation velocity with the difference between the Lagrangian-mean velocity and the pseudo-momentum (Andrews & McIntyre Reference Andrews and McIntyre1978; Gilbert & Vanneste Reference Gilbert and Vanneste2018).

Interpreting the LANS- $\alpha$ system as Lagrangian-averaged Navier–Stokes in this way, this closure allows the Helmholtz operator $\mathcal{S}$ to be commuted with the advective operator, a key property used in the derivation of the slow dynamics. While it is by now standard to regularise Navier–Stokes by applying a number of different filters, this particular filter, with the given meaning of $\alpha$ , is especially advantageous due to its commutative properties with space and time derivatives and advection.

3.1. Abstract framework

As in Embid & Majda (Reference Embid and Majda1996, Reference Embid and Majda1998), Babin et al. (Reference Babin, Mahalov, Nicolaenko and Zhou1997b ), Wingate et al. (Reference Wingate, Embid, Holmes-Cerfon and Taylor2011) and Whitehead & Wingate (Reference Whitehead and Wingate2014), we write the dependent variables as a vector:

(3.1) \begin{equation} \boldsymbol{w} \overset {{\text{def}}}{=} \begin{pmatrix} \boldsymbol{v} \\ \rho \end{pmatrix}. \end{equation}

The non-dimensionalised LANS- $\alpha$ equations (2.8) may then be rewritten in the abstract operator form:

(3.2) \begin{equation} \frac {\partial \boldsymbol{w}}{\partial t} + \frac {1}{\textit{Ro}}\mathcal{L}_{\textit{Ro}}\boldsymbol{w} + \frac {1}{\textit{Fr}}\mathcal{L}_{\textit{Fr}}\boldsymbol{w} + \mathcal{B}(\boldsymbol{w},\boldsymbol{w}) = \mathcal{D}\boldsymbol{w} \end{equation}

(omitting any forcing, as in Embid & Majda (Reference Embid and Majda1998)). In (3.2), the operators are defined:

(3.3) \begin{align} \mathcal{L}_{\textit{Ro}} \boldsymbol{w} &= \begin{pmatrix}\mathcal{S}^{-1}\Big (\boldsymbol{v}_H^\perp + \boldsymbol{\nabla} _H \varDelta ^{-1} \left (\hat {\boldsymbol{z}}\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v})\right )\Big ) \\[8pt] \mathcal{S}^{-1}\Big (\dfrac {\partial }{\partial z}\varDelta ^{-1}(\hat {\boldsymbol{z}}\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v}))\Big ) \\[5pt] 0 \end{pmatrix}, \notag\\ \mathcal{L}_{\textit{Fr}}\boldsymbol{w} &= \begin{pmatrix} - \boldsymbol{\nabla} _H \varDelta ^{-1}\left (\dfrac {\partial \rho }{\partial z}\right ) \notag\\[12pt] -\varDelta ^{-1}\left (\dfrac {\partial ^2\rho }{\partial z^2} \right )+ \rho \\[12pt] -\mathcal{S}^{-1}v_3 \end{pmatrix},\notag\\ \mathcal{B}(\boldsymbol{w}, \boldsymbol{w}) &= \begin{pmatrix} \big(\mathcal{S}^{-1} \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\big) \boldsymbol{v} - \boldsymbol{\nabla }\varDelta ^{-1}\Big (\boldsymbol{\nabla }\boldsymbol{\cdot }\Big(\big(\mathcal{S}^{-1} \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\big) \boldsymbol{v}+ \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }(\mathcal{S}^{-1}\boldsymbol{v})^T\Big)\Big ) \\[6pt] (\mathcal{S}^{-1}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }) \rho \end{pmatrix},\notag \\ \mathcal{D}\boldsymbol{w} &= \frac {1}{\textit{Re}}\begin{pmatrix} \Delta \boldsymbol{v} \\[5pt] \dfrac {1}{Pr}\Delta \rho \end{pmatrix}. \end{align}

To write the operators in this form, we made use of the commutative properties of $\mathcal{S}$ . Each of these operators reduces to those in Embid & Majda (Reference Embid and Majda1998), Babin et al. (Reference Babin, Mahalov, Nicolaenko and Zhou1997b ) and Wingate et al. (Reference Wingate, Embid, Holmes-Cerfon and Taylor2011) in the case of $\alpha =0$ .

We decompose the full solution to (3.2), $\boldsymbol{w}$ , as

(3.4) \begin{equation} \boldsymbol{w} = \boldsymbol{w}_F+\boldsymbol{w}_S, \end{equation}

where the subscript $F$ is for ‘fast’ and the subscript $S$ is for ‘slow’. This decomposition is possible because there exists a projection operator $P$ onto the null space of the fast operator (see Appendix A.2), in a triply periodic domain, such that

(3.5) \begin{equation} P\boldsymbol{w}_S = \boldsymbol{w}_S, \quad P\boldsymbol{w}_F = 0. \end{equation}

In the fast QG- $\alpha$ limit, the fast operator is given by $F\mathcal{L}_{\textit{Ro}}+\mathcal{L}_{\textit{Fr}}$ . We will use the decomposition (3.5) to find evolution equations for the components of the flow and energy on and off the slow manifold in the subsequent sections.

The derivation of both the slow dynamics and the fast dynamics involves a spectral analysis of the fast operator; for simplicity, we consider (3.2) with triply periodic boundary conditions.

3.2. Energy conservation laws

The equation for the global integrated total energy is given by (see Appendix C)

(3.6) \begin{align} \frac {1}{2}\frac {\rm d}{{\rm d}t} \int _\varOmega \Big (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{v} + \rho ^2\Big ) \,\text{d}\boldsymbol{x}&= -\frac {1}{\textit{Re}}\int _\varOmega (\boldsymbol{\nabla }\times \boldsymbol{u})\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v})\, \text{d}\boldsymbol{x}- \frac {1}{\textit{RePr}}\int _\varOmega |\boldsymbol{\nabla }\rho |^2 \,\text{d}\boldsymbol{x}\notag\\&\quad +\int _{\varOmega }\boldsymbol{u}\boldsymbol{\cdot }\mathcal{F}\,\text{d}\boldsymbol{x}. \end{align}

Since the LANS- $\alpha$ equations can be written in the abstract operator form (3.2) and the fast operator is skew-Hermitian, the results of Babin et al. (Reference Babin, Mahalov and Nicolaenko1995, Reference Babin, Mahalov and Nicolaenko1996a , Reference Babin, Mahalov and Nicolaenkob , Reference Babin, Mahalov and Nicolaenko1997a , Reference Babin, Mahalov, Nicolaenko and Zhoub ) and Embid & Majda (Reference Embid and Majda1998), together with the energy conservation equation (3.6) (when there is no dissipation or forcing) yield that

(3.7) \begin{equation} \boldsymbol{u}(0)\boldsymbol{\cdot }\boldsymbol{v}(0) + (\rho (0))^2 = \boldsymbol{v}_F(t)\boldsymbol{\cdot }\boldsymbol{u}_F(t) + (\rho _F(t))^2 + \boldsymbol{v}_S(t)\boldsymbol{\cdot }\boldsymbol{u}_S(t)+(\rho _S(t))^2, \end{equation}

where, in the notation of (3.4), $\boldsymbol{w}_F=(\boldsymbol{v}_F,\rho _F)^T$ , $\boldsymbol{w}_S=(\boldsymbol{v}_S,\rho _S)^T$ and an analogous decomposition for $\mathcal{S}^{-1}\boldsymbol{w}$ defines $\boldsymbol{u}_F$ and $\boldsymbol{u}_S$ . The limiting slow dynamics also conserves energy in the absence of dissipation and forcing, i.e. $\boldsymbol{v}_S(t)\boldsymbol{\cdot }\boldsymbol{u}_S(t)+(\rho _S(t))^2 = \boldsymbol{v}_S(0)\boldsymbol{\cdot }\boldsymbol{u}_S(0)+(\rho _S(0))^2$ . Thus there is a constant ratio of total (as in kinetic plus potential) fast energy ( $\boldsymbol{v}_F\boldsymbol{\cdot }\boldsymbol{u}_F+\rho _F^2$ ) to total slow energy ( $\boldsymbol{v}_S\boldsymbol{\cdot }\boldsymbol{u}_S+\rho _S^2$ ) in time, in the limit as $\epsilon \to 0$ ( $\epsilon =Ro$ ) (in the inviscid, unforced case).

This allows the energy of the slow balanced flow and that of the fast waves to be treated as distinct quantities. The energy partition is necessary to investigate the mechanism by which the Lagrangian-averaging parameter $\alpha$ alters the interplay between slow vortical and fast wave modes.

3.3. Slow dynamics, PV- $\alpha$

To find the slow dynamics, we formulate a projection operator onto the null space of the fast linear operator, as in (3.5). A projection is given by (see Appendix A.1)

(3.8) \begin{equation} P_{{\textit{QG}}_\alpha }\boldsymbol{w} = \begin{pmatrix} \boldsymbol{v}_H - F^2\mathcal{S}^{-2} \varDelta _{{\textit{QG}}_\alpha }^{-1}\dfrac {\partial ^2 \boldsymbol{v}_H}{\partial z^2} - \varDelta _{{\textit{QG}}_\alpha }^{-1}\left (\boldsymbol{\nabla} _H(\boldsymbol{\nabla} _H\boldsymbol{\cdot }\boldsymbol{v}_H) + F\mathcal{S}^{-1} \boldsymbol{\nabla} _H^\perp \left (\dfrac {\partial \rho }{\partial z}\right )\right )\\[8pt] 0 \\[3pt] \rho -F \mathcal{S}^{-1}\varDelta _{{\textit{QG}}_\alpha }^{-1}\left (\dfrac {\partial }{\partial z}(\boldsymbol{\nabla} _H\times \boldsymbol{v}_H)\right ) - \varDelta ^{-1}_{{\textit{QG}}_\alpha }\varDelta _H\rho \end{pmatrix}, \end{equation}

where $\varDelta _{{\textit{QG}}_\alpha } = \varDelta _H + F^2\mathcal{S}^{-2}( {\partial ^2}/{\partial z^2})$ , and $\mathcal{S}^{-2}$ denotes $\mathcal{S}^{-1}$ applied twice. The subscript ‘H’ in (3.8) refers to the horizontal portion of a vector, i.e. $\boldsymbol{v}_H\overset {{\text{def}}}{=} (v_1,v_2)$ . Additionally, $\varDelta _H\overset {{\text{def}}}{=} \partial _{xx}+\partial _{yy}$ and $\boldsymbol{\nabla} _H^\perp \overset {{\text{def}}}{=} (-\partial _y, \partial _x )$ . When $\alpha =0$ , since $\mathcal{S}=\mathcal{S}^{-1}=I$ , this correctly reduces to the projection operator in the QG, $\alpha =0$ case in Whitehead & Wingate (Reference Whitehead and Wingate2014).

The (slow) QG- $\alpha$ dynamics are (see Appendix A.3)

(3.9) \begin{equation} \frac {D^{H,\alpha }}{Dt}\left (\varDelta _H\mathcal{S} \phi + F^2 \mathcal{S}^{-1}\frac {\partial ^2}{\partial z^2}\phi \right )=0, \end{equation}

where $( {D^{H,\alpha }}/{Dt}\overset {{\text{def}}}){=} ( {\partial }/{\partial t} )+ \boldsymbol{u}_H\boldsymbol{\cdot }\boldsymbol{\nabla} _H$ and $\phi$ denotes the potential (here, the modified pressure), so that the vertical vorticity $\omega _3= \mathcal{S} \varDelta _H \phi$ . By using the commutative properties of $\mathcal{S}$ , we can see that analogously to the $\alpha =0$ case, the potential vorticity is simply the $\Delta _{{\textit{QG}}_\alpha }$ operator applied to the potential $\phi$ . We can also rewrite (3.9) as

(3.10) \begin{equation} \frac {D^{H,\alpha }}{Dt}\left (\omega _3 - F \mathcal{S}^{-1}\frac {\partial }{\partial z}\rho \right )=0, \end{equation}

which is very similar to the PV conservation equation in the $\alpha =0$ case:

(3.11) \begin{equation} \frac {D^H}{Dt}\left (\omega _3 - F \frac {\partial }{\partial z}\rho \right )=0, \end{equation}

where the vertical vorticity reduces to $\omega _3 =\varDelta _H p$ as usual, since $\mathcal{S}=I$ in the $\alpha =0$ case. The key differences between (the slow portion of) QG- $\alpha$ , (3.10), and the usual QG, (3.11), are that in QG- $\alpha$ , the density is now filtered, the vorticity has a different relationship to the (modified) pressure, now involving the Helmholtz operator, and the material derivative is with respect to the transport velocity. Thus, the dominant, slow QG- $\alpha$ dynamics are regularised compared with the dominant, slow QG $(\alpha =0)$ dynamics (figure 1): an expected outcome given that LANS- $\alpha$ is a regularised model.

3.4. Fast limiting dynamics and modifications of triad interactions by $\alpha$

While $\alpha$ regularises the dynamics, we find that it alters the balance between the slow and fast dynamics, amplifying the role of the fast dynamics. To show this, we analyse the triad interactions within the LANS- $\alpha$ equations and compare with those of the unregularised equations. Although LANS- $\alpha$ is motivated by Lagrangian averaging, in this paper we interpret $\alpha$ primarily as a scale-selective regularisation: it modifies dynamics at wavelengths comparable to, or smaller than, $\alpha$ , while leaving longer-wave motions essentially unchanged in the resolved fields. A useful parallel is Bretherton’s (Reference Bretherton1969) calculation of the mean flow induced by a short internal-gravity wave packet, where averaging over the fast oscillatory wave phase removes the short-wave oscillations but, in the two-dimensional case, leaves a quasi-hydrostatic radiation field of long internal waves; to leading order, the packet propagates with the group velocity while the total horizontal pseudo-momentum associated with the packet is conserved. Hence, averaging over fast wave scales need not remove all wave motion; in the $\alpha$ -regularised equations, long-wave inertia–gravity modes likewise remain, so it is natural to consider how $\alpha$ shifts their dispersion relations and the associated resonance conditions.

The dispersion relations for the modes associated with the fast QG- $\alpha$ limit are found by seeking eigenfunctions of the fast operator $\mathcal{L}_F$ of the form

(3.12) \begin{equation} \boldsymbol{w} = \exp (i\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x} - i\omega (\boldsymbol{k})t)\boldsymbol{r}, \end{equation}

where $\boldsymbol{r}$ is an associated eigenvector, $\omega (\boldsymbol{k})$ is an associated eigenfrequency and $\boldsymbol{k}=(k_1,k_2,k_3)$ is a wavenumber. In this context, $\omega (\boldsymbol{k})$ will be either zero or the classic inertia–gravity wave dispersion relation. We also denote

(3.13) \begin{equation} |\boldsymbol{k}|\overset {{\text{def}}}{=} \left(k_1^2+k_2^2+k_3^2\right)^{{1}/{2}} \quad \text{and} \quad |\boldsymbol{k}_H|\overset {{\text{def}}}{=}\left(k_1^2+k_2^2\right)^{ {1}/{2}}. \end{equation}

Analytic formulas for the eigenfrequencies are required to characterise all possible three-wave resonances. In the usual fast QG limit, the eigenfrequencies are given by (Embid & Majda Reference Embid and Majda1998)

(3.14) \begin{equation} \omega ^{(0)}(\boldsymbol{k})=0 \quad \text{(double root)} \quad \text{ and } \quad \omega ^{(\pm 1)}(\boldsymbol{k}) = \pm \frac {\left(|\boldsymbol{k}_H|^2+F^2k_3^2\right)^{ {1}/{2}}}{|\boldsymbol{k}|}. \end{equation}

Physically, the zero-frequency (slow) modes are zero-frequency Rossby waves, while the modes with frequency $\omega ^{(\pm 1)}$ are (dispersive) inertia–gravity waves. The parameter $F=f/N$ appears in the dispersion relation (3.14) and sets the relative importance of vertical versus horizontal structure. For instance, when stratification is stronger than rotation ( $N\gt f)$ , the impact of vertical variations on the inertia–gravity wave frequency is reduced.

Table 1.

Summary of the four types of relevant resonant and near-resonant interactions. Note that Fourier coefficients for fast–fast–slow interactions are zero (Embid & Majda Reference Embid and Majda1998), which does not change when $\alpha \neq 0$ . For wavenumbers $\boldsymbol{k}+\boldsymbol{p}=\boldsymbol{q}$ , the eigenfrequency condition for (all but slow–slow–slow) resonances to occur changes with the parameter $\alpha$ . The resonance types are listed in order of their frequency and importance for energy transfer in the case $\alpha =0$ . Slow–slow–slow interactions are the most numerous but are unaffected by $\alpha$ . Fast–fast–fast interactions have been less considered because they occur only sparsely, though this changes when $\alpha \neq 0$ . For more details, see Appendix B.2.

Table 1. Long description

Table 1 has four columns: resonance type, α = 0 condition, α ≠ 0 condition, and role in energy evolution. The four rows list slow-slow-slow, slow-fast-fast, fast-fast-fast, and slow-slow-fast interactions. The table shows that slow-slow-slow interactions occur for all wavenumbers and are unaffected by α. Slow-fast-fast and fast-fast-fast resonance conditions change when α is nonzero, because the dispersion relation is modified. Exact slow-slow-fast resonances remain absent, although near-resonance locations may be distorted.

In the fast QG- $\alpha$ limit, the eigenfrequencies are modified by the non-dimensional Helmholtz operator in spectral space, $(1+L_\alpha |\boldsymbol{k}|^2)$ , where $L_\alpha =\alpha ^2/L^2$ is the non-dimensionalised version of $\alpha ^2$ . The eigenfrequencies when $\alpha \neq 0$ ( $L_\alpha \neq 0$ ) are given by

(3.15) \begin{align} \omega ^{(0)}(\boldsymbol{k})=0 \quad \text{(double root)} \quad \text{ and } \quad \omega ^{(\pm 1)}(\boldsymbol{k}) = \pm \frac {\left(\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}_H|^2+F^2k_3^2\right)^{ {1}/{2}}}{\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}|}. \end{align}

Smaller $L_\alpha$ corresponds to $\alpha$ -regularised dynamics that is closer to the unregularised dynamics. In the limit $L_\alpha =0$ , the regularisation vanishes; (3.15) then reduces to (3.14).

The set of resonant triads is described by

(3.16) \begin{equation} R_{\beta , \boldsymbol{q}}^\alpha \overset {{\text{def}}}{=} \left\{\big(\boldsymbol{k}, \boldsymbol{p}, \beta ', \beta ''\big)\;|\; \boldsymbol{k}+\boldsymbol{p}=\boldsymbol{q}, \; \omega ^{(\beta ')}(\boldsymbol{k})+\omega ^{(\beta '')}(\boldsymbol{p}) = \omega ^{(\beta )}(\boldsymbol{q})\right\}\!, \end{equation}

where the superscript $\alpha$ emphasises that the set of possible interactions now depends on $\alpha$ due to (3.15). The relevant classes of triads are described in table 1. Of these, we consider the slow–fast–fast and fast–fast–fast triads, as the set of possible slow–slow–slow resonances does not change with $\alpha$ and there are no slow–slow–fast resonances or near-resonances when both $\alpha =0$ and $\alpha \neq 0$ (for conservative tolerance levels, see supplementary material, figures S1–S3 and S7, available at https://doi.org/10.1017/jfm.2026.11720.).

Slow–fast–fast interactions, when $F\approx 1$ ( ${\textit{Bu}}\approx 1$ ), have been seen to result in a rapid downscale energy transfer, and have been described as catalytic interactions responsible for the emergence of a geostrophically adjusted state (Bartello Reference Bartello1995). When the eigenfrequencies are given by (3.14), the condition for the pair of wavenumbers $\boldsymbol{k}$ and $\boldsymbol{p}$ , with $\boldsymbol{k}+\boldsymbol{p}=\boldsymbol{q}$ , to be an element of $R_{\beta ,\boldsymbol{q}}^0$ (with $\beta '=\pm 1, \beta ''=0, \beta =\pm 1$ ), as defined in (3.16), is

(3.17) \begin{equation} \pm \frac {\left(|\boldsymbol{k}_H|^2+F^2k_3^2\right)^{{1}/{2}}}{|\boldsymbol{k}|} = \pm \frac {\left(|\boldsymbol{q}_H|^2+F^2q_3^2\right)^{{1}/{2}}}{|\boldsymbol{q}|}, \end{equation}

or similarly with $\boldsymbol{p}$ in place of $\boldsymbol{k}$ , interchanging the roles of $\beta '$ and $\beta ''$ . The condition (3.17) produces a resonance for any $\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}$ when $F=1$ or a near-resonance when $F=O(1)$ . When the eigenfrequencies are instead given by (3.15), the condition is

(3.18) \begin{equation} \pm \frac {\left(\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}_H|^2+F^2 k_3^2\right)^{{1}/{2}}}{\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}|} = \pm \frac {\left(\left(1+L_\alpha |\boldsymbol{q}|^2\right)|\boldsymbol{q}_H|^2+F^2 q_3^2\right)^{{1}/{2}}}{\left(1+L_\alpha |\boldsymbol{q}|^2\right)|\boldsymbol{q}|}, \end{equation}

or similarly with $\boldsymbol{p}$ upon interchanging the roles of $\beta '$ and $\beta ''$ .

Due to the additional nonlinearity of (3.18) as compared with (3.17), it is challenging to reason analytically exactly how the number and distribution of the elements of $R_{\beta ,\boldsymbol{q}}^\alpha$ change with $\alpha$ (or $L_\alpha$ ), but it is possible to do so numerically. Our results show that for $F=1$ ( ${\textit{Bu}}=1$ ), the introduction of a non-zero $\alpha$ drastically curtails the number of slow–fast–fast resonant triads, including the cross-scale interactions responsible for transferring energy between large-scale and small-scale dynamics (figure 2). This reduction disrupts the primary channel for downscale energy transfer. Further increasing $L_\alpha$ leads to an increase in the number of resonances, though it does not return to the number when $L_\alpha =0$ for reasonable values of $L_\alpha$ (given that in modelling studies, typically $\alpha$ is only a few grid points in length, i.e. $\alpha = 2\Delta x$ ). For the regimes in which this energy pathway is already less active, $F=1/2$ ( ${\textit{Bu}}=4$ ) and $F=2$ ( ${\textit{Bu}}=1/4$ ), the introduction of non-zero $\alpha$ has a different effect, only increasing the number of resonances and thus altering rather than shutting down energy exchange via slow–fast–fast resonances.

Figure 2.

Regularisation by $\alpha$ ( $L_\alpha =\alpha ^2/L^2$ ) reshapes the landscape of resonant interactions, superseding the ratio $F=f/N$ (thus also the Burger number, ${\textit{Bu}}=N^2/f^2$ ) as the dominant parameter determining the relative amounts of different resonances, once $L_\alpha$ is sufficiently large. As $L_\alpha$ increases, the counts of both slow–fast–fast and fast–fast–fast resonant and near-resonant interactions converge to values largely independent of $F$ . All counts were computed on a grid of $256^3$ using a tolerance of $10^{-5}$ for numerically checking the resonance conditions (see supplementary materials for the robustness across different tolerance levels). For each line style (fixed $F$ ), both slow–fast–fast (black) and fast–fast–fast (grey) triad counts are computed using that same $F$ in the dispersion relation for all three modes. The slow–fast–fast curves for different $F$ are indistinguishable, given sufficiently large $L_\alpha$ , on the logarithmic $y$ axis and appear as a single curve. The fast–fast–fast curves are similarly close, for all $L_\alpha$ , and appear as a single (overlapping) curve. Given the large number of wavenumbers associated with the $256^3$ grid, it was possible for the resonance conditions to be satisfied a large number of times (up to trillions, as shown here). To focus on general three-dimensional interactions, triads involving purely horizontal or purely vertical wavenumbers were excluded from the count.

Figure 2. Long description

Figure 2 is a line plot with L_α on the horizontal axis and number of interactions on the logarithmic vertical axis. Black markers and curves represent slow-fast-fast triads, while grey markers and curves represent fast-fast-fast triads. Different line styles represent F = 1/2, F = 1, and F = 2. At L_α = 0, the slow-fast-fast counts depend strongly on F, while fast-fast-fast counts are absent or much smaller. As L_α increases, the fast-fast-fast counts rise sharply, and the counts for both interaction classes become much less dependent on F.

Figure 3.

A detailed view of the change in the number of resonant interactions for small values of $L_\alpha$ , for three values of the rotation-to-stratification ratio $F=f/N$ ( $F=1/2,1,2$ ), shows that the number of fast–fast–fast triads appears to grow continuously with $L_\alpha$ , rather than emerging in a discontinuous jump. Concurrently, the count of slow–fast–fast interactions exhibits an inflection point with $L_\alpha$ . These counts are computed in the same manner as those in figure 2, but with respect to a range of smaller $L_\alpha$ values.

Figure 3. Long description

Figure 3 contains three narrow plots arranged horizontally and labelled F = 1/2, F = 1, and F = 2. Each panel has L_α on the horizontal axis and number of interactions on the logarithmic vertical axis. Black curves show slow-fast-fast counts, and grey curves show fast-fast-fast counts. The grey fast-fast-fast curves rise from zero or near zero as L_α increases. The black slow-fast-fast curves bend downward and then upward, showing an inflection in the count at small L_α values.

Fast–fast–fast exact resonant interactions do not occur when $1/2\leq F\leq 2$ in the $\alpha =0$ case, and resonant and near-resonant interactions occur only rarely otherwise. The condition is

(3.19) \begin{equation} \pm \frac {\left(|\boldsymbol{k}_H|^2+F^2k_3^2\right)^{ {1}/{2}}}{|\boldsymbol{k}|} \pm \frac {\left(|\boldsymbol{p}_H|^2+F^2p_3^2\right)^{ {1}/{2}}}{|\boldsymbol{p}|}= \pm \frac {\left(|\boldsymbol{q}_H|^2+F^2q_3^2\right)^{ {1}/{2}}}{|\boldsymbol{q}|}. \end{equation}

When $\alpha \neq 0$ , (3.19) becomes instead

(3.20) \begin{align} \pm \frac {\left(\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}_H|^2+F^2 k_3^2\right)^{ {1}/{2}}}{\left(1+L_\alpha |\boldsymbol{k}|^2\right)|\boldsymbol{k}|} &\pm \frac {\left(\left(1+L_\alpha |\boldsymbol{p}|^2\right)|\boldsymbol{p}_H|^2+F^2 p_3^2\right)^{ {1}/{2}}}{\left(1+L_\alpha |\boldsymbol{p}|^2\right)|\boldsymbol{p}|} \notag\\&= \pm \frac {\left(\left(1+L_\alpha |\boldsymbol{q}|^2\right)|\boldsymbol{q}_H|^2+F^2 q_3^2\right)^{ {1}/{2}}}{\left(1+L_\alpha |\boldsymbol{q}|^2\right)|\boldsymbol{q}|}. \end{align}

The new nonlinearity and dependence on wavenumber in (3.20) allows the number of fast–fast–fast resonances to significantly increase, even in cases where it is usually zero without $\alpha$ -regularisation. When $\alpha \neq 0$ , the number of fast–fast–fast resonant interactions becomes comparable to the number of slow–fast–fast resonant interactions (figure 2) for every value of $F$ ( $F={1}/{2},1,2$ ). This stands in stark contrast to the usual fast QG limit, where fast–fast–fast triads are typically sparse or non-existent. The emergence of a large number of fast–fast–fast resonances allows energy to be efficiently exchanged and redistributed purely among the fast modes, independent of the slow vortical flow. A detailed view of this emergence (figure 3) confirms that the number of fast–fast–fast triads appears to grow continuously from $L_\alpha =0$ rather than spiking discontinuously; this is especially evident in the $F=1/2$ ( ${\textit{Bu}}=4$ ) and $F=2$ ( ${\textit{Bu}}=1/4$ ) plots.

Perhaps most significantly, as $\alpha$ becomes sufficiently large, the numbers of both slow–fast–fast triads and fast–fast–fast triads converge to values largely independent of the rotation-to-stratification ratio $F$ , indicating that $\alpha$ supersedes $F$ as the dominant parameter shaping the potential for resonant interactions. The overall effect of $\alpha$ on the fast dynamics is not to simply regularise but to restructure the energy pathways.

4.1. Description of the numerical method

For all numerical experiments, we use Dedalus, a pseudo-spectral solver in Python (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The non-dimensional LANS- $\alpha$ equations (2.8a – c ) are solved within the triply periodic domain $[0,1]^3$ , discretised on a uniform grid of $N^3=256^3$ points with a $3/2$ dealiasing rule. The system is initialised from a state of rest. For time stepping, we employ an explicit second-order, two-stage Runge–Kutta scheme (RK222). This scheme was selected for computational efficiency and validation runs confirmed its results were comparable to those from the third-order, four-stage scheme (RK443). The time step is adjusted dynamically by a Courant–Friedrichs–Lewy condition determined by the maximum total frequency corresponding to the unregularised velocity $\boldsymbol{v}$ , across the grid (see Burns et al. (Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) for more details).

The physical set-up of our experiments is similar to that of Smith & Waleffe (Reference Smith and Waleffe2002). Energy is injected by selecting a divergence-free forcing $\mathcal{F}$ in the momentum equation whose spectral shape is inspired by Smith & Waleffe (Reference Smith and Waleffe2002), but which is constant in time (as in Wingate et al. Reference Wingate, Embid, Holmes-Cerfon and Taylor2011). An isotropic amplitude distribution is prescribed via the scalar spectrum $f(|\boldsymbol{k}|)$ , which is a Gaussian in wavenumber space, centred at the forcing wavenumber $k_f$ :

(4.1) \begin{equation} f(|\boldsymbol{k}|) = \frac {\epsilon _f}{(2\pi )^{{1}/{2}}s}\exp \left (\frac {-(|\boldsymbol{k}| -k_f)^2}{2s^2}\right )\!. \end{equation}

Specifically, the coefficients are normalised so that the shell-wise mean-square amplitude follows $f(|\boldsymbol{k}|)$ , so $\epsilon _f$ is referred to as the forcing-amplitude parameter (it is not generally a realised energy injection rate, except under specific scaling choices). The forcing-amplitude parameter is $\epsilon _f=1$ and the standard deviation is $s=1$ for all simulations. Since we are interested in the QG limit, we use low wavenumber forcing, $k_f=3(2\pi )$ . The choice of a time-independent forcing is motivated primarily by the need for grid-scale convergence studies. One issue is that the LANS- $\alpha$ model presents known stability challenges, and verifying grid-scale convergence at each value of $L_\alpha$ is essential. Another issue is that $\alpha$ itself modifies the effective dissipation range (Chen et al. Reference Chen, Holm, Margolin and Zhang1999), so that the resolution requirements change with $L_\alpha$ . Together, these make convergence verification both essential and non-trivial. A stochastic-in-time forcing as in Smith & Waleffe (Reference Smith and Waleffe2002) and Sukhatme & Smith (Reference Sukhatme and Smith2008) compounds both difficulties, as random realisations can push the model outside its region of stability unpredictably. These considerations lead us to choose forcing and viscosity for the study in this paper that are fixed in time. This forcing defines a characteristic time scale:

(4.2) \begin{equation} \tau = \left(\epsilon _f k_f^2\right)^{-{1}/{3}}. \end{equation}

All simulations are integrated to a final non-dimensional time of $t_f=1$ , corresponding to approximately $7\tau$ . In the simulations, we are particularly interested in whether $\alpha$ modifies the first crossover time (as in Sukhatme & Smith Reference Sukhatme and Smith2008), which is the first time at which the energy in the geostrophic (or vortical) modes exceeds that in the wave modes. Comparing the crossover times for different values of $\alpha$ , all else equal, provides a measure of how rapidly the partition of energy reorganises and thus how $\alpha$ modifies the relative effectiveness of wave-sustaining pathways. This will also show how the $\alpha$ -modification of the triad interactions (§ 3.4) manifests in the energy evolution. We choose the final time of $t_f=1$ for the simulations in accordance with this goal of observing the initial development of the turbulent flow involving the transfer of energy between vortical and wave modes.

In (2.8a ), we use a hyperviscosity and set $\mathcal{D}_\nu \boldsymbol{v} = \nu (-1)^{n+1}\boldsymbol{\nabla} ^{2n} \boldsymbol{v}$ and $\mathcal{D}_\kappa \rho = \kappa (-1)^{n+1} \boldsymbol{\nabla} ^{2n} \rho$ with $n=8$ , as in Smith & Waleffe (Reference Smith and Waleffe2002), except the coefficients $\nu$ and $\kappa$ are constant in time and defined based on the dealiased grid cutoff wavenumber, to allow for more direct comparisons across simulations of different grid sizes.

We investigate the system’s behaviour for three values of $F$ : 1 (equally strong rotation and stratification), 2 (stronger rotation) and 1/2 (stronger stratification), achieved by setting $(Ro, Fr) = (0.1, 0.1)$ , $(0.1, 0.2)$ and $(0.2, 0.1)$ , respectively. In all cases, both $Ro$ and ${\textit{Fr}}$ are small in order to reflect the QG limit. For each $F$ , we systematically vary the non-dimensional regularisation parameter across five values: $L_\alpha = 0, 0.01, 0.1, 0.15, 0.2$ .

A representative snapshot from the simulations is shown in figure 4, which also illustrates how the distinction between the velocities $\boldsymbol{v}$ and $\boldsymbol{u}$ manifests under LANS- $\alpha$ regularisation.

Figure 4.

Representative snapshots of horizontal velocity fields from the LANS- $\alpha$ simulations. Panels show the first horizontal component over a horizontal plane at $z=0.5$ and $t=0.5$ : (a) the unregularised $x$ -velocity $v_1$ for $L_\alpha =0$ , (b) the unregularised $x$ -velocity $v_1$ for $L_\alpha =0.1$ and (c) the associated regularised $x$ -velocity $u_1$ for $L_\alpha =0.1$ . In these snapshots, the $L_\alpha =0.1$ field exhibits less pronounced fine-scale structure in $v_1$ than the $L_\alpha =0$ case, and $u_1$ shows reduced small-scale contrast and smaller spatial variability in magnitude. The colourbar indicates the value of the plotted velocity component (m s $^{-1}$ ), and a common colour scale is used across all panels.

Figure 4. Long description

Figure 4 has three square greyscale heat-map panels with x and y axes from 0 to 1 and a shared horizontal colorbar. Panels (a) and (b) are on the top row, and panel (c) is centered below them. Panel (a) shows v1 for L_α = 0, panel (b) shows v1 for L_α = 0.1, and panel (c) shows u1 for L_α = 0.1. Dark and light regions indicate negative and positive values of the velocity component. Compared with panel (a), panel (b) has less fine-scale texture, and panel (c) is smoother and lower-contrast than panel (b).

Figure 5.

The regularisation parameter $L_\alpha$ extends the persistence of wave energy by delaying the dominance of vortical modes, as shown here for runs with $F=1$ ( ${\textit{Bu}}=1$ ). (a) panel shows the decomposition of domain-integrated total energy (kinetic plus potential) into slow vortical energy (dashed) and wave energy (solid) for increasing $L_\alpha$ . (b) panel quantifies this delay, showing that the crossover time (when vortical energy surpasses wave energy) is a monotonically increasing function of $L_\alpha$ .

Figure 5. Long description

Figure 5 has two panels for the F = 1 simulations. The larger left panel plots total energy versus time, with solid curves for wave energy and dashed curves for vortical energy; markers distinguish L_α values from 0 to 0.2. The narrow right panel plots crossover time versus L_α, where crossover means the first time vortical energy exceeds wave energy. The right-panel curve rises with L_α, showing that increasing L_α delays the crossover.

Figure 6.

Energy decomposition for $F=1$ ( ${\textit{Bu}}=1$ ) showing the proportions of the domain-integrated total energy lying in slow vortical modes (dashed) and fast wave modes (solid) over time, corresponding to the magnitudes in figure 5. The crossover time is delayed with increasing $L_\alpha$ .

Figure 6. Long description

Figure 6 is a single panel line plot for the F = 1 simulations. The horizontal axis is time, and the vertical axis is proportion of total energy. Solid curves show wave-energy fractions, and dashed curves show vortical-energy fractions. The outermost dashed and solid lines form a shape akin to an hourglass open on both sides, with the leftmost half of the hourglass narrow compared to the rightmost half. Marker styles distinguish L_α values from 0 to 0.2. For small L_α, the wave fraction decreases quickly while the vortical fraction increases quickly. For larger L_α, both changes are slower, so wave energy remains a larger share for longer.

4.2. Evolution of slow vortical and fast wave energy

The numerical simulations illustrate how the modifications to the set of resonant triads analysed in § 3.4 manifest in the flow’s energy evolution. To see this, we track the partition of total energy between the vortical modes of the slow dynamics and the fast wave modes. Note that the total energy for the $L_\alpha \neq 0$ cases is given by $( {1}/{2})\int _\varOmega (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{v} + \rho ^2) \,\text{d}\boldsymbol{x}$ as opposed to $({1}/{2})\int _\varOmega (\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v} + \rho ^2) \,\text{d}\boldsymbol{x}$ , as the former is the conserved quantity in (3.6) when there is no dissipation or forcing. The former is smaller in magnitude than the latter (figure 4), especially for large $L_\alpha$ , because $\boldsymbol{u}$ is modified according to $\boldsymbol{v} = (1-L_\alpha \Delta )\boldsymbol{u}$ (the energy $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}$ for the $L_\alpha \neq 0$ cases better illustrates how the LANS- $\alpha$ simulations are indeed more energetic overall).

Figure 7.

The decomposition of domain-integrated total (kinetic plus potential) energy into slow vortical energy (dashed) and fast wave energy (solid) for increasing $L_\alpha$ , for $F=1/2$ ( ${\textit{Bu}}=4$ ). The right-hand panel quantifies this delay, showing that the first crossover time (when vortical energy first surpasses wave energy) is a monotonically increasing function of $L_\alpha$ .

Figure 7. Long description

Figure 7 has two panels for the F = 1/2 simulations. The larger left panel plots total energy versus time, with solid curves for fast wave energy and dashed curves for slow vortical energy; markers distinguish increasing L_α values. The narrow right panel plots first crossover time versus L_α. The crossover-time curve increases with L_α, indicating that regularization delays the transition to vortical dominance.

Figure 8.

Energy decomposition for $F=1/2$ ( ${\textit{Bu}}=4$ ) showing the proportions of the domain-integrated total energy lying in slow vortical modes (dashed) and fast wave modes (solid) over time, corresponding to the magnitudes in figure 7.

Figure 8. Long description

Figure 8 is a single panel line plot for the F = 1/2 simulations. The horizontal axis is time, and the vertical axis is proportion of total energy. Solid curves show wave-mode fractions, and dashed curves show vortical-mode fractions. The outermost dashed and solid lines form a shape akin to an hourglass open on both sides, with the leftmost half of the hourglass narrow compared to the rightmost half, similarly to that in Figure 6 except that both halves of the hourglass are narrower than those of Figure 6. Marker styles distinguish L_α values from 0 to 0.2. As L_α increases, the decrease in wave-energy fraction and increase in vortical-energy fraction occur more slowly.

The decomposition is computed by performing an eigenvalue decomposition of the fields outputted by the simulations, as in Smith & Waleffe (Reference Smith and Waleffe2002), according to the eigenvalues in § 3.4. This method is equivalent to applying the projection operator described in § 3.3.

The most immediate effect of the changes in the fast dynamics, as shown in the total energy plots for each value of $F$ (figures 5, 7 and 9), is a consistent delay in the evolution towards a state dominated by slow vortical modes. For all values of $F$ considered, increasing $L_\alpha$ systematically extends the crossover time at which the energy in vortical modes surpasses that in the wave modes (right-hand panels of figures 5, 7 and 9). This directly corresponds to a delay of geostrophic adjustment, and by hindering the flow’s relaxation towards the balanced state, more energy is retained in the wave modes for a longer duration.

Figure 9.

The decomposition of domain-integrated total (kinetic plus potential) energy into slow vortical energy (dashed) and fast wave energy (solid) for increasing $L_\alpha$ , for $F=2$ ( ${\textit{Bu}}=1/4$ ). The right-hand panel shows that in this case, the crossover time (when vortical energy surpasses wave energy) is no longer a monotonically increasing function of $L_\alpha$ and instead has an inflection point. Simulations for additional values of $L_\alpha$ ( $0.025,\ 0.05,\ 0.075$ ) were performed for the $F=2$ case to more clearly illustrate this unique behaviour. When $L_\alpha$ exceeds 0.05, there is no crossover within the simulation period (non-dimensional time $t_f=1$ ), which is why the crossover times for those $L_\alpha$ values exceed the bounds of the plot. In all other simulations, the crossover time was achieved well within the simulation time, and there was no inflection point.

Figure 9. Long description

Figure 9 has two panels for the F = 2 simulations. The larger left panel plots total energy versus time, with solid curves for fast wave energy and dashed curves for slow vortical energy; markers distinguish L_α values from 0 to 0.2, including additional small L_α values. The narrow right panel plots crossover time versus L_α. Unlike Figures 5 and 7, the right-panel curve first decreases and then increases. For L_α values above about 0.05, the figure indicates that crossover occurs after the final simulated time.

The emergence of fast–fast–fast resonant triads offers a potential mechanism to explain this delay. As shown in § 3.4, the introduction of a non-zero $L_\alpha$ creates a dense network of these interactions, which are sparse or non-existent in the $L_\alpha =0$ case. This new resonant pathway allows for the efficient exchange and redistribution of energy purely among the wave modes. Such a self-sustaining wave field, where energy is actively transferred among fast modes rather than being propagated away, is consistent with the observed delay in the system’s relaxation to a geostrophic mode-dominated state. The creation of this fast–fast–fast pathway is a universal outcome of the regularisation across all values of $F$ , accompanied by a reorganisation of slow–fast–fast triads whose number is sharply reduced for $F=1$ ( ${\textit{Bu}}=1$ ) but increased for $F \neq 1$ . For the ${\textit{Bu}}=1$ case, where geostrophic adjustment is extremely efficient in the absence of regularisation, the disruption of this slow–fast–fast pathway likely plays an important, complementary role in producing the observed delay (figures 5 and 6).

Figure 10.

Energy decomposition for $F=2$ ( ${\textit{Bu}}=1/4$ ) showing the proportions of domain-integrated total energy lying in slow vortical modes (dashed) and fast wave modes (solid) over time, corresponding to the magnitudes in figure 9.

Figure 10. Long description

Figure 10 is a single panel line plot for the F = 2 simulations. The horizontal axis is time, and the vertical axis is proportion of total energy. Solid curves show wave-mode fractions, and dashed curves show vortical-mode fractions. Marker styles distinguish several L_α values from 0 to 0.2. The curves show a more complicated response than in Figures 6 and 8: some wave fractions decrease and then level off or rise, while the corresponding vortical fractions increase more slowly or non-monotonically.

Although the regularisation parameter $L_\alpha$ becomes the dominant factor (as opposed to $F$ ) shaping the fast dynamics (e.g. figure 2), the system’s energy evolution still exhibits clear dependencies on the relative strengths of rotation and stratification. In the unregularised case, the dynamics varies considerably, as expected: the $F=2$ flow (figure 10) is most persistently wave-dominated, while the $F=1/2$ case (figure 8) becomes dominated by the slow modes most rapidly. Yet, the response to the regularisation is also not uniform. For $F=1/2$ and $F=1$ , the crossover time is a monotonically increasing function of $L_\alpha$ (figures 7 and 5). The $F=2$ case, however, exhibits a non-monotonic response, where a small amount of regularisation ( $L_\alpha =0.01$ ) actually accelerates the onset of vortical dominance before larger values delay it.

The regularisation parameter $L_\alpha$ acts as a second control on the system’s energetics, in the sense that it governs the evolution of the partition of energy in a manner analogous to $F$ (and ${\textit{Bu}}$ ). While higher values of both parameters ultimately limit the maximum proportion of vortical energy (see Sukhatme & Smith (Reference Sukhatme and Smith2008) for a comparison of the impact of many different values of $F$ in the $\alpha =0$ case), their effects on the initial distribution of wave and vortical energy are opposite. As shown in figures 6, 8 and 10, higher values of $L_\alpha$ lead to a smaller initial portion of wave energy – the opposite effect from that of $F$ .

Over longer time scales, the growth in total energy is primarily driven by the slow vortical modes. This is evident from the steeper slope of the vortical energy curves after the initial adjustment period. However, the wave energy does not simply saturate or decay after the crossover. While in some cases it does appear to stagnate for brief periods, it generally continues to increase throughout the simulations, albeit at a slower rate than its vortical counterpart.

These results show that $\alpha$ acts as a second control parameter for the system’s energetic pathways, with effects that complement and sometimes oppose those of the rotation-to-stratification ratio, delaying geostrophic adjustment and sustaining energy growth in both wave and vortical components.