1 Introduction↩︎

The Lagrangian-Averaged Navier-Stokes-\(\alpha\) (LANS\(-\alpha\)) equations constitute a turbulence closure model that has achieved at lower resolution what other more commonly used closure schemes require higher resolution to accomplish, in a number of contexts [1]–[6]. Notably, LANS-\(\alpha\) produces more energetic simulations [7], [8]. This improves heat and salinity transport resolution [5] and turbulence statistics [9] in ocean models and representations of eddy momentum flux in atmosphere models [10]. The \(\alpha\)-regularization also lowers the critical wavenumber for baroclinic instability [11], improving its ability to model the transfer of potential to kinetic energy on coarser grids.

The LANS-\(\alpha\) model is fundamentally different from classical closure schemes: it regularizes nonlinear advection rather than introducing an eddy diffusivity. This functions to filter the nonlinear dynamics at length scales smaller than the namesake parameter \(\alpha\). The effectiveness of the LANS-\(\alpha\) model, as well as its key physical and mathematical properties — such as a Kelvin circulation theorem [12], a Kármán-Howarth theorem [13], a conserved energy [14], bounded enstrophy [15], well-posedness [16], and the existence of a finite-dimensional global attractor [17], [18] — are perhaps unsurprising, given that its formulation emerged independently from several lines of inquiry. For instance, while its mathematical roots trace back to Leray’s seminal work, in which regularizing Navier-Stokes is the key idea needed to prove the existence of weak solutions [19], it more directly originates from the Gjaja-Holm Wave-Mean Flow Interaction (WMFI) equations [20]. The LANS-\(\alpha\) equations also share mathematical connections with closure schemes that parameterize eddy stress based on deformation [21]–[23]. As in those methods, numerical instabilities can sometimes arise when implementing the LANS-\(\alpha\) model, limiting its broader application [6], [8], [23]. Both the significant advantages and the primary obstacle of the LANS-\(\alpha\) model lead to the central question of this work: can the model’s energetic enhancement be understood through the parameter \(\alpha\)’s role in modifying the underlying wave dynamics and energy transfers?

While there are partial explanations for why the LANS-\(\alpha\) model is more energetic — varying \(\alpha\) modifies the dissipation range [1], a property shared by related \(\alpha\)-regularization models [24], [25], and introduces an effective Rossby deformation radius, enabling baroclinic instability to be resolved at coarser resolution [11] — further insight into the processes driving this increased energy would benefit both the theoretical understanding of the model and its implementation. To this end, we investigate the energetics of the LANS-\(\alpha\) model within the quasi-geostrophic (QG) limit.

The QG equations, originally derived by [26], describe the large-scale, balanced dynamics primarily driven by planetary rotation. The three-dimensional, continuously stratified system considered here is further constrained by strong background stratification. By design, the QG equations effectively filter out fast, small-scale motions like inertia-gravity waves, isolating the dominant “balanced" flow. The first derivation of the QG limiting dynamics for the Boussinesq-\(\alpha\) equations appears in [27]. We extend [27] by presenting the first fast singular limit derivation of QG-\(\alpha\) for the Boussinesq-\(\alpha\) model. In our work, the leading-order asymptotic solution is the conservation of PV-\(\alpha\) modulated by the leading-order fast dynamics, which is found by deriving an operator to project the full dynamics onto the null space of the fast operator. Using this framework, we show there exists a rich dynamical structure involving interactions purely among fast waves (fast-fast-fast dynamics) that is absent in the unregularized (\(\alpha = 0\)) case (Figure 1). Further, the prognostic variable PV-\(\alpha\) derived in this limit is intrinsically useful for future studies of LANS-\(\alpha\); it is the central dynamical variable of the balanced system and a quantity fundamental to understanding large-scale geophysical flows (e.g., [26]; [28]; [29]; [30]; [14]; [31]). In other words, we derive QG-\(\alpha\) by formulating a projection operator onto the balanced flow. This projection serves a dual purpose: it also provides the foundation for our primary goal of analyzing the energetics of LANS-\(\alpha\) by enabling the partitioning of the system’s energy between the balanced and fast wave dynamics.

The complete elimination of fast dynamics without affecting the evolution of the mean flow, as is intended in Charney’s QG equations, has been described as the existence of a “slow manifold," a concept furthered by [32] and [33]. Ideally, a reduced equation resulting from an asymptotic expansion in physical parameters of interest would correspond to a slow manifold: attracting other dynamics toward balanced states, representing dominant dynamics in a lower-dimensional subspace, and ensuring that trajectories initialized on this subspace remain there. In reality, this ideal is rarely achieved, and the”slow manifolds” that can be obtained are better described as “fuzzy" slow manifolds [34]. Even if an original set of equations permits a corresponding slow manifold in the mathematical sense, it may not serve as the desired physical slow manifold, e.g., in the QG limit, one where gravity wave activity is entirely absent [35]–[37]. Indeed, the QG equations do not represent a perfectly invariant subspace, since even perfectly balanced initial conditions can generate exponentially small inertia-gravity waves [38]. Furthermore, even though fast and slow dynamics are decoupled at leading order, they are connected in subtle but dynamically significant ways. Resonant interactions allow slow (also called potential vorticity or geostrophic) modes to affect the distribution of fast waves [39]. When forcing is not restricted to the slow dynamics, a common scenario in realistic settings, fast dynamics can significantly impact the distribution of energy throughout the system [40]–[42]. A robust analysis of the energetics of the LANS-\(\alpha\) model thus requires not only deriving the slow QG-\(\alpha\) equations but also accounting for the influence of the fast dynamics.

To do so, we depart from methods that strictly isolate a slow manifold and instead take a fast singular limit [43]–[48]. By taking the same parameter limit as traditional QG – low Rossby number (\(Ro\to 0\), geostrophic balance) and low Froude number (\(Fr\to 0\), hydrostatic balance), while the ratio between them (\(F= Fr/Ro=f/N\)) remains fixed – this approach recovers the standard QG equations for the slow dynamics while retaining fast dynamics at leading order. Applying this process to the LANS-\(\alpha\) equations yields a system we term the “fast QG-\(\alpha\) dynamics" (or fast QG-\(\alpha\) limit). The fast QG-\(\alpha\) dynamics decompose into the slow dynamics, the balanced geostrophic flow corresponding to the traditional QG equations (but with \(\alpha\) regularization), and the fast dynamics, the inertia-gravity waves filtered out in the traditional limit. The fast dynamics are decoupled from the usual slow dynamics in the sense that the leading-order fast dynamics cannot influence the slow dynamics because these fast, oscillatory effects cancel out when averaged across many wave periods. In this context, the fluctuations about the balanced QG-\(\alpha\) state are the \(O(1)\) fast dynamics together with all \(O(\epsilon)\) dynamics, where the small parameter \(\epsilon\) is taken to be the Rossby number. Physically, the frequency of the fast-wave oscillations is inversely proportional to \(\epsilon\). This approach [46], [47] is rigorously justified by Schochet’s [49] method of cancellation of oscillations for hyperbolic PDEs and may be thought of as averaging over fast waves in geophysical flows.

The three-dimensional Boussinesq (or in this work, Boussinesq-\(\alpha\)) equations admit two kinds of eigenfrequencies: slow (geostrophic) modes (vortical modes) which have zero frequency for all wavenumbers \(\boldsymbol{k}\) and fast (ageostrophic) modes (dispersive waves). While potential vorticity in the fast QG (or QG-\(\alpha\)) limit is determined solely by slow modes, in other balanced limits such as low Rossby/finite Froude or low Froude/finite Rossby, waves associated with the slower physical process can contribute to the slow dynamics [47], [50], [51].

To understand what drives the LANS-\(\alpha\) model’s distinctive energetics, we analyze three-wave resonant and near-resonant interactions. Resonant interactions between slow and fast modes, i.e., interactions between geostrophic turbulence and waves, are a primary mechanism for geostrophic adjustment via a direct cascade of wave energy [52]. These resonant interactions are highly dependent on the Burger number (\(Bu=Ro^2/Fr^2\)), which is related to the ratio \(F=f/N\) via \(Bu=1/F^2\). For example, [50] show analytically that as stratification increases compared to rotation, ageostrophic energy cascades are “unfrozen," allowing nonlinear geostrophic adjustment to occur. With small-scale forcing, geostrophic dynamics dominate at large scales for \(1/2\leq F\leq 2\) (or equivalently, \(1/4\leq Bu\leq 4\)), a range with no triad interactions exclusively involving fast waves [40]. With large-scale forcing, the system transitions from being vortical-energy dominated for \(F\leq 1\) (\(Bu\geq 1\)) to being wave-energy dominated for \(F>1\) (\(Bu<1\)) [53]. Near-resonant interactions, while more challenging to analyze [54] and weak for \(F\sim O(1)\) [55], also play an important role in the transfer of energy to slow wave modes [41].

Our numerical simulations reveal that the smoothing parameter \(\alpha\) plays a role analogous to that of the rotation-to-stratification ratio \(F\) and affects wave-vortex energy ratios in a similar manner. However, a key difference emerges when varying \(\alpha\) as opposed to \(F\) (or \(Bu\)): increasing \(\alpha\) markedly extends the persistence of wave energy, thereby slowing the typical transition to vortical energy dominance. A likely mechanism for this is found in how \(\alpha\) reshapes three-wave resonant interactions. In the standard fast QG limit, resonant interactions in which two fast waves produce another fast wave are nonexistent for many Burger numbers and sparse otherwise. In contrast, by numerically tabulating how often our analytically-derived eigenfrequencies satisfy the condition for resonance interactions, we find that \(\alpha\)-regularization drastically increases the number of these interactions (as illustrated in Figure 1). While \(\alpha\) was previously known to modify the dissipation range, these alterations to the triad interactions reveal that the introduction of \(\alpha\) can both cause energy to accumulate within the waves and decrease resonant interactions generally responsible for cross-scale energy transfer, providing a mechanism by which energy may be prevented from reaching dissipative scales and thus accumulate within the system. Collectively, our results indicate that LANS-\(\alpha\) is likely more energetic than other models because of the way \(\alpha\)-regularization modulates fast waves, affecting their persistence and ability to satisfy resonance conditions that drive energy transfer among the slow and fast components of the flow.

The remainder of this paper is organized as follows. In \(\S 2\), we introduce the LANS-\(\alpha\) equations, including their non-dimensional, non-local form, and interpretation as Lagrangian-averaged Navier-Stokes. The derivation of the fast QG limit for LANS-\(\alpha\) is presented in \(\S 3\); this section includes energy conservation laws, the slow dynamics (conservation of PV-\(\alpha\)), the fast limiting dynamics (\(O(1)\) wave dynamics), and an analysis of how the parameter \(\alpha\) modifies three-wave resonance interactions. In \(\S 4\), we present numerical simulations focusing on the evolution of slow and wave energies in the QG limit and how it depends on \(\alpha\). Finally, in \(\S 5\), we discuss the main findings and their broader implications.

2 The LANS-\(\alpha\) equations↩︎

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 [12], [14]: \[\tag{1} \begin{align} \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{u} \cdot \nabla \boldsymbol{v} + \boldsymbol{v}\cdot \nabla \boldsymbol{u}^T + f\hat{\boldsymbol{z}}\times \boldsymbol{u} + \nabla \phi+\frac{\rho}{\rho_0}g\hat{\boldsymbol{z}} &= \nu \nabla^{2n} \boldsymbol{v}, \tag{2} \\ \frac{\partial \rho}{\partial t} + \boldsymbol{u} \cdot \nabla \rho -\boldsymbol{u}\cdot b\hat{\boldsymbol{z}} &=\kappa \nabla^{2n} \rho, \tag{3} \\ \nabla \cdot \boldsymbol{u} &= 0, \\ \boldsymbol{v} &= \boldsymbol{u} - \alpha^2 \Delta \boldsymbol{u}, \tag{4} \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, \(\nu\) is kinematic viscosity, \(b\) is the strength of the underlying density stratification, \(\kappa\) is buoyancy diffusivity, and \(\alpha\in \mathbb{R}\) is a parameter with units of length. While this length scale \(\alpha\) can be interpreted in many ways – such as the typical deviation of a particle trajectory from its time-averaged path [2] or a filter width [56] – we adopt the interpretation that it is the length scale below which the dynamics are regularized. This regularization is accomplished through the introduction of a new velocity field, \(\boldsymbol{u}\), as defined in 4 . In this framework, the smoothed field \(\boldsymbol{u}\) advects the momentum \(\boldsymbol{v}\), mitigating the development of sharp gradients at scales smaller than \(\alpha\).

The difference between the Boussinesq-\(\alpha\) equations and the usual Boussinesq equations is essentially due to the introduction of the smoothed velocity \(\boldsymbol{u}=(u_1, u_2, u_3)\). Denoting the Helmholtz operator as \[\label{def:helmholtz-operator} \mathcal{S}\overset{\text{def}}{=}(1-\alpha^2\Delta )\tag{5}\] (\(\mathcal{S}\) for “smoothing"), the smoothed velocity \(\boldsymbol{u}\) is constructed via \(\mathcal{S}^{-1}\boldsymbol{v}\) 4 . The smoothed velocity contributes an additional advective term, \(\boldsymbol{v} \cdot \nabla \boldsymbol{u}^T\) (also written as \((\nabla \boldsymbol{u})^T \boldsymbol{v}\) or \(v_j\partial_i u_j\)), in 2 and necessitates a modification to the pressure. The usual pressure \(p\) has been replaced with the modified pressure \(\phi\) in 2 , where \[\label{eq:modified-pressure} \phi \overset{\text{def}}{=}\frac{p}{\rho_0} - \frac{1}{2}|\boldsymbol{u}|^2-\frac{\alpha^2}{2}|\nabla \boldsymbol{u}|^2.\tag{6}\] From the definition of \(\mathcal{S}\) 5 , 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 \(\Delta\) is the scalar Laplacian rather than the vector Laplacian. Using the Helmholtz operator to smooth the advecting velocity does not introduce dissipative effects, since 2 and 3 are time-reversible when \(\nu=\kappa=0\).

In 2 4 , the density has been decomposed as \[\rho_{\text{full}} = \overline{\rho}+\rho, \quad \overline{\rho} = p_0 - bz, \quad b>0.\]

The dissipation terms in 2 are presented in a general form, and when \(n=1\), the dissipation is of the usual Laplacian type. 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 to those of previous turbulence studies and the classical fast singular QG limit [47].

We non-dimensionalize 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\). Note that throughout this work, \(N\) is considered constant. We non-dimensionalize \(\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 (c.f. [57]). Using these scalings, the Rossby, Froude, Reynolds, Prandtl, and Euler numbers arise naturally: \[Ro \overset{\text{def}}{=}\frac{U}{fL}, \quad 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}.\] The parameter \(F\), which measures the relative strength of rotation versus stratification, \[\label{def:F} F\overset{\text{def}}{=}\frac{Fr}{Ro}=\frac{f}{N},\tag{7}\] also arises in the system and must be assumed fixed when taking the QG limit. This parameter is closely related to the Burger number \[\label{def:Bu} Bu\overset{\text{def}}{=}\frac{Ro^2}{Fr^2}=\frac{1}{F^2},\tag{8}\] and will play an important role in the dynamics.

Omitting the primes, the non-dimensionalized version of 1 is given by \[\tag{9} \begin{align} \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{u}\cdot \nabla \boldsymbol{v} +\boldsymbol{v} \cdot \nabla \boldsymbol{u}^T+ \frac{1}{Ro} \hat{\boldsymbol{z}} \times \boldsymbol{u} + \nabla \phi + \Gamma\rho\hat{\boldsymbol{z}} &= \frac{1}{Re}\nabla^{2n} \boldsymbol{v},\tag{10} \\ \frac{\partial \rho}{\partial t} + \boldsymbol{u}\cdot\nabla\rho - \frac{1}{Fr}\boldsymbol{u}\cdot\hat{\boldsymbol{z}} &= \frac{1}{RePr}\nabla^{2n} \rho, \tag{11}\\ \nabla \cdot \boldsymbol{u} &= 0, \tag{12} \\ \boldsymbol{u}&=\mathcal{S}^{-1}\boldsymbol{v}. \end{align}\] Here, \(\mathcal{S}\) is actually the non-dimensionalized Helmholtz operator, \((1-L_\alpha^2 \Delta)\), where \(L_\alpha \overset{\text{def}}{=}\alpha/L\). 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\) is also the non-dimensional version. The constant \(\Gamma\overset{\text{def}}{=}BgLU^{-2}\) but we will use the scaling \(\Gamma = Fr^{-1}\). We will consider \(n=1\) in the diffusivity for simplicity.

We obtain the non-local form of 9 by solving the elliptic equation for the modified pressure, so that \[\begin{align} \label{eq:lans-nondim-2} \nabla \phi = \nabla\Delta^{-1}\left(\frac{1}{Ro}\hat{\boldsymbol{z}}\cdot \mathcal{S}^{-1} \boldsymbol{\omega} -\nabla\cdot (\mathcal{S}^{-1}\boldsymbol{v}\cdot \nabla \boldsymbol{v} + \boldsymbol{v} \cdot \nabla (\mathcal{S}^{-1}\boldsymbol{v})^T)-\frac{1}{Fr}\rho\boldsymbol{\hat{z}}\right). \end{align}\tag{13}\] The modified pressure thus does not further differentiate the Bousinessq-\(\alpha\) equations from the unregularized Boussinesq equations; \(\phi\) is a free variable in exactly the same sense \(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 [47], [50], and [51] upon setting \(\alpha=0\).

2.2 Interpretation as Lagrangian-Averaged Navier-Stokes↩︎

The definition of the smoothed velocity 4 indicates that \(\boldsymbol{u}\) should be interpreted as the Lagrangian average of \(\boldsymbol{v}\) to justify that this system is “Lagrangian-averaged" Navier-Stokes. It must then be the case that the Helmholtz operator, applied to a Lagrangian-averaged field, yields an Eulerian-averaged field, and vice versa with the inverse. When defined as the Helmholtz operator dependent on the scalar \(\alpha\), \(\mathcal{S}\) is an approximation to the true operator for which this is (asymptotically) true. This operator, the dynamical Helmholtz operator, is given by \((1-\tilde{\Delta})\) where \(\tilde{\Delta} = \nabla \cdot \langle \boldsymbol{\xi}\boldsymbol{\xi}\rangle \cdot \nabla\) [14]. Here, \(\boldsymbol{\xi}\) represents a displacement fluctuation of a particle trajectory \(\boldsymbol{X}\), so that \(\langle \boldsymbol{\xi}\boldsymbol{\xi}\rangle\) is the covariance of the displacement fluctuation. It can be seen to \(o(|\boldsymbol{\xi}|^2)\) that the quantity \[(\nabla \cdot\langle \boldsymbol{\xi}\boldsymbol{\xi}\rangle \cdot \nabla) \langle \boldsymbol{U}\rangle^L,\] where \(\langle \boldsymbol{U}\rangle^L\) is a Lagrangian mean velocity, is the Stokes drift. Thus, the dynamical Helmholtz operator \((1-\tilde{\Delta})\) can be interpreted as mapping Lagrangian-averaged fields to their Eulerian-averaged counterparts.

However, the LANS-\(\alpha\) equations do not involve the dynamical Helmholtz operator; \(\tilde{\Delta}\) is replaced with \(\alpha^2 \Delta\). This assumes the isotropy condition \(\langle \xi^k \xi^l\rangle = \alpha^2 \delta^{kl}\), where \(\alpha\) is constant if the statistics are homogeneous. Thus, by assuming statistically stationary and isotropic turbulence (in a sense), LANS-\(\alpha\) can be interpreted as Lagrangian-averaged Navier-Stokes. This interpretation allows the dynamical Helmholtz operator to be commuted with the advective operator, a key property used in the derivation the slow dynamics. While it is by now standard to regularize 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 The fast quasi-geostrophic limit of the LANS-\(\alpha\) equations↩︎

3.1 Abstract framework↩︎

As in [46], [47], [50], [51], and [42], we write the dependent variables as a vector, \[\boldsymbol{w} \overset{\text{def}}{=}\begin{pmatrix} \boldsymbol{v} \\ \rho\end{pmatrix}.\] The non-dimensionalized LANS-\(\alpha\) equations 9 may then be rewritten in the abstract operator form \[\label{eq:abstract-operator-form} \frac{\partial \boldsymbol{w}}{\partial t} + \frac{1}{Ro}\mathcal{L}_{Ro}\boldsymbol{w} + \frac{1}{Fr}\mathcal{L}_{Fr}\boldsymbol{w} + \mathcal{B}(\boldsymbol{w},\boldsymbol{w}) = \mathcal{D}\boldsymbol{w}.\tag{14}\] In 14 , the operators are defined: \[\label{eq:lans-abstract-operators} \begin{align} \mathcal{L}_{Ro} \boldsymbol{w} &= \begin{pmatrix}\mathcal{S}^{-1}\Big(\boldsymbol{v}_H^\perp+ \nabla_H \Delta^{-1} \left(\hat{z}\cdot(\nabla \times \boldsymbol{v})\right)\Big) \\ \mathcal{S}^{-1}\Big(\frac{\partial}{\partial z}\Delta^{-1}(\hat{z}\cdot(\nabla \times \boldsymbol{v}))\Big) \\ 0 \end{pmatrix}, \\ \mathcal{L}_{Fr}\boldsymbol{w} &= \begin{pmatrix} - \nabla_H \Delta^{-1}\left(\frac{\partial \rho}{\partial z}\right) \\ -\Delta^{-1}\left(\frac{\partial^2\rho}{\partial z^2} \right)+ \rho \\ -\mathcal{S}^{-1}w \end{pmatrix},\\ \mathcal{B}(\boldsymbol{w}, \boldsymbol{w}) &= \begin{pmatrix} (\mathcal{S}^{-1} \boldsymbol{v} \cdot \nabla) \boldsymbol{v} - \nabla \Delta^{-1}\Big(\nabla \cdot ((\mathcal{S}^{-1} \boldsymbol{v} \cdot \nabla) \boldsymbol{v}+ \boldsymbol{v} \cdot \nabla (\mathcal{S}^{-1}\boldsymbol{v})^T)\Big) \\ (\mathcal{S}^{-1}\boldsymbol{v}\cdot \nabla) \rho \end{pmatrix}, \\ \mathcal{D}\boldsymbol{w} &= \frac{1}{Re}\begin{pmatrix} \Delta \boldsymbol{v} \\ \frac{1}{Pr}\Delta \rho\end{pmatrix}. \end{align}\tag{15}\] 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 [47], [50], and [51] in the case of \(\alpha=0\).

We decompose the full solution to 14 , \(\boldsymbol{w}\), as \[\label{eq:slow-fast-decomp} \boldsymbol{w} = \boldsymbol{w}_F+\boldsymbol{w}_S,\tag{16}\] 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 6.2), in a triply periodic domain, such that \[\label{eq:fast-slow-decomp-projection} P\boldsymbol{w}_S = \boldsymbol{w}_S, \quad P\boldsymbol{w}_F = 0.\tag{17}\] In the fast QG-\(\alpha\) limit, the fast operator is given by \(F\mathcal{L}_{Ro}+\mathcal{L}_{Fr}\). We will use the decomposition 17 to find evolution equations for the components of the flow and energy on and off the slow manifold in the proceeding sections.

Both the derivation of the slow dynamics and the fast dynamics involves a spectral analysis of the fast operator; for simplicity, we consider 14 with triply periodic boundary conditions.

3.2 Energy conservation laws↩︎

The equation for the global integrated total energy is given by (see Appendix 8) \[\label{eq:total-energy-conservation} \frac{1}{2}\frac{d}{dt} \int_\Omega \Big(\boldsymbol{u}\cdot \boldsymbol{v} + \rho^2\Big) \,\text{d}\boldsymbol{x}= -\frac{1}{Re}\int_\Omega (\nabla\times \boldsymbol{u})\cdot(\nabla \times \boldsymbol{v}) - \frac{1}{RePr}(\nabla \rho \cdot \nabla\rho) \,\text{d}\boldsymbol{x}.\tag{18}\]

Since the LANS-\(\alpha\) equations can be written in the abstract operator form 14 and the fast operator is skew-Hermitian, the results of [47] and [45], [50], [58]–[60] together with the energy conservation equation 18 yield that \[\label{eq:energy-decomposition} \boldsymbol{u}(0)\cdot\boldsymbol{v}(0) + \rho(0)^2 = \boldsymbol{v}_F\cdot \boldsymbol{u}_F + \rho_F^2 + \boldsymbol{v}_S\cdot\boldsymbol{u}_S+\rho_S^2,\tag{19}\] where in the notation of 16 , \(\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\). Thus the ratio of total (as in, kinetic plus potential) fast and slow energy is conserved in the limit as \(\epsilon\to 0\) (\(\epsilon=Ro\)).

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 17 . A projection is given by (see Appendix 6.1) \[\label{qg-proj-alpha-neq0} P_{\text{QG}_\alpha}\boldsymbol{w} = \begin{pmatrix} \boldsymbol{v}_H - F^2\mathcal{S}^{-2} \Delta_{\text{QG}_\alpha}^{-1}\frac{\partial^2 \boldsymbol{v}_H}{\partial z^2} - \Delta_{\text{QG}_\alpha}^{-1}\left(\nabla_H(\nabla_H\cdot \boldsymbol{v}_H) + F\mathcal{S}^{-1} \nabla_H^\perp\left(\frac{\partial \rho}{\partial z}\right)\right)\\ 0 \\ \rho -F \mathcal{S}^{-1}\Delta_{\text{QG}_\alpha}^{-1}\left(\frac{\partial}{\partial z}(\nabla_H\times \boldsymbol{v}_H)\right) - \Delta^{-1}_{\text{QG}_\alpha}\Delta_H\rho\end{pmatrix},\tag{20}\] where \(\Delta_{\text{QG}_\alpha} = \Delta_H + F^2\mathcal{S}^{-2}\frac{\partial^2}{\partial z^2}\), and \(\mathcal{S}^{-2}\) denotes \(\mathcal{S}^{-1}\) applied twice. When \(\alpha=0\), since \(\mathcal{S}=\mathcal{S}^{-1}=I\), this correctly reduces to the projection operator in the QG, \(\alpha=0\) case in [42].

The (slow) QG-\(\alpha\) dynamics are (see Appendix 6.3) \[\label{eq:qg-alpha-v1} \frac{D^{H,\alpha}}{Dt}\left(\Delta_H\mathcal{S} \phi + F^2 \mathcal{S}^{-1}\frac{\partial^2}{\partial z^2}\phi\right)=0,\tag{21}\] where \(\frac{D^{H,\alpha}}{Dt}\overset{\text{def}}{=}\frac{\partial}{\partial t} + \boldsymbol{u}_H\cdot\nabla\), and \(\phi\) denotes the potential (here, the modified pressure), so that the vertical vorticity \(\omega_3= \mathcal{S} \Delta_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_{\text{QG}_\alpha}\) operator applied to the potential \(\phi\). We can also rewrite 21 as \[\label{eq:qg-alpha-v2} \frac{D^{H,\alpha}}{Dt}\left(\omega_3 - F \mathcal{S}^{-1}\frac{\partial}{\partial z}\rho\right)=0,\tag{22}\] which is very similar to the PV conservation equation in the \(\alpha=0\) case, \[\label{eq:qg-alphazero} \frac{D^H}{Dt}\left(\omega_3 - F \frac{\partial}{\partial z}\rho\right)=0,\tag{23}\] where the vertical vorticity reduces to \(\omega_3 =\Delta_H p\) as usual, since \(\mathcal{S}=I\) in the \(\alpha=0\) case. The key differences between (the slow portion of) QG-\(\alpha\), 22 , and the usual QG, 23 , are that in QG-\(\alpha\), the density is now smoothed, the vorticity has a different relationship to the (modified) pressure, now involving the Helmholtz operator, and the material derivative is with respect to the smoothed velocity. Thus, the dominant, slow QG-\(\alpha\) dynamics are regularized compared to the dominant, slow QG \((\alpha=0)\) dynamics (Figure 1): an expected outcome given that LANS-\(\alpha\) is a regularized model.

3.4 Fast limiting dynamics and modifications of triad interactions by \(\alpha\)↩︎

While \(\alpha\) regularizes 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 analyze the triad interactions within the LANS-\(\alpha\) equations and compare to those of the unregularized equations. 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 \[\label{def:eigenfunctions} \boldsymbol{w} = \exp(i\boldsymbol{k}\cdot \boldsymbol{x} - i\omega(\boldsymbol{k})t)\boldsymbol{r},\tag{24}\] 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 \[\label{def:wavenumber-mag} |\boldsymbol{k}|\overset{\text{def}}{=}(k_1^2+k_2^2+k_3^2)^{\frac{1}{2}} \quad \text{and} \quad |\boldsymbol{k}_H|\overset{\text{def}}{=}(k_1^2+k_2^2)^{\frac{1}{2}}.\tag{25}\] Analytic formulas for the eigenfrequencies are required to characterize all possible three-wave resonances. In the usual fast QG limit, the eigenfrequencies are given by [47]: \[\label{eq:eigenvalues-QG} \omega^{(0)}(\boldsymbol{k})=0 \quad \text{(double root)} \quad \text{ and } \quad \omega^{(\pm 1)}(\boldsymbol{k}) = \pm \frac{(|\boldsymbol{k}_H|^2+F^2k_3^2)^{\frac{1}{2}}}{|\boldsymbol{k}|}.\tag{26}\] 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 26 and sets the relative importance of vertical versus horizontal structure. For instance, when stratification is stronger than rotation (\(N>f)\), the impact of vertical variations on the inertia-gravity wave frequency is reduced.

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/L\) is the non-dimensionalized version of \(\alpha\). The eigenfrequencies when \(\alpha \neq 0\) (\(L_\alpha \neq 0\)) are given by \[\label{eq:eigenvalues-QG-alpha} \omega^{(0)}(\boldsymbol{k})=0 \quad \text{(double root)} \quad \text{ and } \quad \omega^{(\pm 1)}(\boldsymbol{k}) = \pm \frac{((1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}_H|^2+F^2k_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}|}.\tag{27}\] The smaller \(L_\alpha\) is, the closer the \(\alpha\)-regularized dynamics are to the unregularized dynamics. When \(L_\alpha=0\), the \(\alpha\)-regularized dynamics reduce to the usual dynamics; notice that 27 reduces to 26 when \(L_\alpha=0\).

The set of resonant triads is described by \[\label{def:resonant-set} R_{\beta, \boldsymbol{q}}^\alpha \overset{\text{def}}{=}\{(\boldsymbol{k}, \boldsymbol{p}, \beta', \beta'')\;|\; \boldsymbol{k}+\boldsymbol{p}=\boldsymbol{q}, \; \omega^{(\beta')}(\boldsymbol{k})+\omega^{(\beta'')}(\boldsymbol{p}) = \omega^{(\beta)}(\boldsymbol{q})\},\tag{28}\] where the superscript \(\alpha\) emphasizes that the set of possible interactions now depends on \(\alpha\) due to 27 . 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 exact slow-slow-fast resonances both when \(\alpha=0\) and \(\alpha\neq 0\).

Slow-fast-fast interactions, when \(F\approx 1\) (\(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 [52]. When the eigenfruencies are given by 26 , the condition for wavenumbers \(\boldsymbol{k}\) and \(\boldsymbol{p}\) such that \(\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 28 , is \[\label{eq:slow-fast-fast-cond-alphazero} \pm \frac{(|\boldsymbol{k}_H|^2+F^2k_3^2)^{\frac{1}{2}}}{|\boldsymbol{k}|} = \pm \frac{(|\boldsymbol{q}_H|^2+F^2q_3^2)^{\frac{1}{2}}}{|\boldsymbol{q}|},\tag{29}\] or similarly with \(\boldsymbol{p}\) in place of \(\boldsymbol{k}\), interchanging the roles of \(\beta'\) and \(\beta''\). The condition 29 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 27 , the condition is \[\label{eq:slow-fast-fast-cond-alpha} \pm \frac{((1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}_H|^2+F^2 k_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}|} = \pm \frac{((1+L_\alpha|\boldsymbol{q}|^2)|\boldsymbol{q}_H|^2+F^2 q_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{q}|^2)|\boldsymbol{q}|},\tag{30}\] or similarly with \(\boldsymbol{p}\) upon interchanging the roles of \(\beta'\) and \(\beta''\).

Due to the additional nonlinearity of 30 as compared to 29 , it is challenging to reason analytically exactly how the number and distribution of the elements of \(R_{\beta,\boldsymbol{q}}^\alpha\) changes with \(\alpha\) (or \(L_\alpha\)), but it is possible to do so numerically. Our results show that for \(F=1\) (\(Bu=1\)), the introduction of a nonzero \(\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 modeling 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\) (\(Bu=4\)) and \(F=2\) (\(Bu=1/4\)), the introduction of nonzero \(\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.

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 \[\label{eq:fast-fast-fast-cond-alphazero} \pm \frac{(|\boldsymbol{k}_H|^2+F^2k_3^2)^{\frac{1}{2}}}{|\boldsymbol{k}|} \pm \frac{(|\boldsymbol{p}_H|^2+F^2p_3^2)^{\frac{1}{2}}}{|\boldsymbol{p}|}= \pm \frac{(|\boldsymbol{q}_H|^2+F^2q_3^2)^{\frac{1}{2}}}{|\boldsymbol{q}|}.\tag{31}\] When \(\alpha \neq 0\), 31 becomes instead \[\label{eq:fast-fast-fast-cond-alpha} \begin{align} \pm \frac{((1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}_H|^2+F^2 k_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{k}|^2)|\boldsymbol{k}|} \pm &\frac{((1+L_\alpha|\boldsymbol{p}|^2)|\boldsymbol{p}_H|^2+F^2 p_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{p}|^2)|\boldsymbol{p}|} \\&= \pm \frac{((1+L_\alpha|\boldsymbol{q}|^2)|\boldsymbol{q}_H|^2+F^2 q_3^2)^{\frac{1}{2}}}{(1+L_\alpha|\boldsymbol{q}|^2)|\boldsymbol{q}|}. \end{align}\tag{32}\]

The new nonlinearity and dependence on wavenumber in 32 allows the number of fast-fast-fast resonances to significantly increase, even in cases where it is usually zero without \(\alpha\)-regularization. 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). This stands in stark contrast to the usual fast QG limit, where fast-fast-fast triads are typically sparse or nonexistent. 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\) (\(Bu=4\)) and \(F=2\) (\(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 smooth but to restructure the energy pathways.

4 Numerical simulations↩︎

4.1 Description of the numerical method↩︎

For all numerical experiments, we use Dedalus, a pseudo-spectral solver in Python [61]. The non-dimensional LANS-\(\alpha\) equations 10 12 are solved within the triply periodic domain \([0,1]^3\), discretized on a uniform grid of \(N^3=256^3\) points with a \(3/2\) dealiasing rule. The system is initialized 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 CFL condition.

The physical setup of our experiments is similar to that of [40]. Energy is injected into the system via a stochastic forcing that is constant in time. The forcing spectrum, \(F(k)\), is a Gaussian in wavenumber space, centered at the forcing wavenumber \(k_f\): \[F(k) = \frac{\epsilon_f}{(2\pi)^{\frac{1}{2}}s}\exp\left(\frac{-(k-k_f)^2}{2s^2}\right).\] The energy injection rate is \(\epsilon_f=1\) and the standard deviation is \(s=1\) for all simulations. Since we are interested in the quasi-geostrophic limit, we use low wavenumber forcing, \(k_f=3(2\pi)\). This forcing defines a characteristic eddy turnover time, \[\tau = (\epsilon_f k_f^2)^{-\frac{1}{3}}.\] All simulations are integrated to a final non-dimensional time of \(t=1\), corresponding to approximately \(7\tau\). This duration is sufficient for observing the initial development of the turbulent flow involving the transfer of energy between vortical and wave modes. In 10 , we have hyperviscosity with \(n=8\), as in [40], except the coefficient is constant in time and defined based on the dealiased grid cutoff wavenumber, in order to allow for more direct comparisons across simulations of different grid sizes.

We investigate the system’s behavior 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 \(Fr\) are small in order to reflect the QG limit. For each \(F\), we systematically vary the nondimensional regularization parameter across five values: \(L_\alpha = 0, 0.01, 0.1, 0.15, 0.2\).

4.2 Evolution of slow vortical and fast wave energy↩︎

The numerical simulations illustrate how the modifications to the set of resonant triads analyzed in \(\S\)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 \(\frac{1}{2}\int_\Omega (\boldsymbol{u}\cdot \boldsymbol{v} + \rho^2) \,\text{d}\boldsymbol{x}\) as opposed to \(\frac{1}{2}\int_\Omega (\boldsymbol{v}\cdot \boldsymbol{v} + \rho^2) \,\text{d}\boldsymbol{x}\), as the former is the conserved quantity in 18 . The former is smaller in magnitude than the latter, especially for large \(L_\alpha\), because \(\boldsymbol{u}\) is modified according to \(\boldsymbol{v} = (1-L_\alpha^2\Delta)\boldsymbol{u}\) (the energy \(\boldsymbol{v}\cdot\boldsymbol{v}\) for the \(L_\alpha\neq 0\) cases better illustrate how the LANS-\(\alpha\) simulations are indeed more energetic overall).

The decomposition is computed by performing an eigenvalue decomposition of the fields outputted by the simulations, as in [40], according to the eigenvalues in \(\S\)3.4. This method is equivalent to applying the projection operator described in \(\S\)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 4, 6, and 8), 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 4, 6, and 8). 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.

The emergence of fast-fast-fast resonant triads offers a potential mechanism to explain this delay. As shown in \(\S\)3.4, the introduction of a nonzero \(L_\alpha\) creates a dense network of these interactions, which are sparse or nonexistent 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 regularization across all values of \(F\), contrasting with the disruption of slow-fast-fast triads that is expected to be most impactful only in the \(F=1\) (\(Bu=1\)) case. For the \(Bu=1\) case, where geostrophic adjustment is extremely efficient in the absence of regularization, the disruption of this slow-fast-fast pathway likely plays an important, complementary role in producing the observed delay (Figures 4 and 5).

Although the regularization 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 unregularized case, the dynamics vary considerably, as expected: the \(F=2\) flow (Figure 9) is most persistently wave-dominated, while the \(F=1/2\) case (Figure 7) becomes dominated by the slow modes most rapidly. Yet, the response to the regularization is also not uniform. For \(F=1/2\) and \(F=1\), the crossover time is a monotonically increasing function of \(L_\alpha\) (Figures 6 and 4). The \(F=2\) case, however, exhibits a non-monotonic response, where a small amount of regularization (\(L_\alpha=0.01\)) actually accelerates the onset of vortical dominance before larger values delay it.

The regularization 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 \(Bu\)). While higher values of both parameters ultimately limit the maximum proportion of vortical energy (see [53] 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 5, 7, and 9, higher values of \(L_\alpha\) lead to a smaller initial portion of wave energy – the opposite effect from that of \(F\).

Over longer timescales, 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.

5 Conclusions↩︎

In this paper, we derive the fast QG limit for the LANS-\(\alpha\) model to analyze its energy evolution. In the absence of forcing and dissipation, the slow dynamics amount to the conservation of PV-\(\alpha\). We find that \(\alpha\) has the expected smoothing effect on these dominant, slow dynamics, but also that it amplifies the role of the fast wave dynamics. The regularization parameter \(\alpha\) alters the frequencies of the waves and, in turn, the conditions for three-wave resonances, which play a key role in the distribution and transfer of energy in the system. By numerically testing integer wavenumbers, we find that introducing \(\alpha\) allows for a rapid increase in the number of resonant triads where two fast waves produce a third fast wave. We also find that the number of slow-fast-fast resonant interactions greatly decreases in the \(F=1\) case, unlike the creation of a large number of fast-fast-fast resonant interactions, which is universal across different values of \(F\) (Figures 2 and 3). Since the fast and slow dynamics are decoupled, we derive an energy conservation equation partitioned into slow and fast energies. Numerical simulations of the LANS-\(\alpha\) system for three values of the rotation-to-stratification ratio \(F\) and a range of \(\alpha\) (\(L_\alpha\)) values produce fields that, once decomposed into slow and fast portions, exhibit energy dynamics consistent with this modification to the resonances.

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 \(F\), indicating that \(\alpha\) supersedes \(F\) (and the Burger number) as the dominant parameter shaping the potential for resonant interactions. This reshaping both fundamentally alters the pathway for geostrophic adjustment and creates a rich wave-wave interaction regime. Indeed, the inflection point seen in the slow-fast-fast triad count (Figure 3) (and in the crossover time in Figure 8) may hint at a critical threshold for this behavior; exploring whether this value of \(L_\alpha\) separates a stable regime from an unstable one, and its dependence on \(F\), merits further investigation.

Beyond its implications for the physics of the LANS-\(\alpha\) model, this surge of fast wave interactions points to a potential physical origin for the model’s numerical instabilities. We propose that such a dense field of interacting waves poses a challenge to numerical schemes because, for example, numerical discretizations can cause dispersion relations to produce both physical modes and computational modes, and the latter can lead to numerical instabilities [62]. A model whose physics permits many more interactions between the physical modes may, in its numerical implementation, also provide more opportunities for the computational modes to interact with each other and the physical modes to trigger numerical instability. This supports the hypothesis that the origin for the known numerical instabilities of the LANS-\(\alpha\) and related models is nonlinear, consistent with linear analyses indicating the model is stable (e.g., [6]). In fact, a (linear) stability analysis of the shallow water LANS-\(\alpha\) model indicates that its numerical implementation should be more stable than that of its unregularized counterpart in the sense that it permits a larger maximum allowable timestep as \(\alpha\) increases [63]. However, such linear analyses do not account for how numerical errors in individual waves can combine via triad interactions to further compromise the accuracy of the nonlinear dynamics, which can be a significant source of error depending on the choice of time-stepping scheme [64]. Our results on enhanced wave-wave interactions indicate that future work, when using the LANS-\(\alpha\) model or related models [21], [22], could not only consider the inherent stability of the computational modes (e.g., the thorough analysis in [65]), but also how they interact with other computational or physical modes in the context of fast-fast-fast interactions.

Ultimately, this analysis demonstrates that by reshaping the conditions for resonance, the LANS-\(\alpha\) model permits a previously inaccessible regime of fast wave interactions which in turn affect the model’s energy evolution, though the full implications of these dynamics and how they change in other relevant parameter regimes merit further study.

References↩︎