Turbulence and Rossby Wave Dynamics with Realizable Eddy Damped Markovian Anisotropic Closure

Abstract

The theoretical basis for the Eddy Damped Markovian Anisotropic Closure (EDMAC) is formulated for two-dimensional anisotropic turbulence interacting with Rossby waves in the presence of advection by a large-scale mean flow. The EDMAC is as computationally efficient as the Eddy Damped Quasi Normal Markovian (EDQNM) closure, but, in contrast, is realizable in the presence of transient waves. The EDMAC is arrived at through systematic simplification of a generalization of the non-Markovian Direct Interaction Approximation (DIA) closure that has its origin in renormalized perturbation theory. Markovian Anisotropic Closures (MACs) are obtained from the DIA by using three variants of the Fluctuation Dissipation Theorem (FDT) with the information in the time history integrals instead carried by Markovian differential equations for two relaxation functions. One of the MACs is simplified to the EDMAC with analytical relaxation functions and high numerical efficiency, like the EDQNM. Sufficient conditions for the EDMAC to be realizable in the presence of Rossby waves are determined. Examples of the numerical integration of the EDMAC compared with the EDQNM are presented for two-dimensional isotropic and anisotropic turbulence, at moderate Reynolds numbers, possibly interacting with Rossby waves and large-scale mean flow. The generalization of the EDMAC for the statistical dynamics of other physical systems to higher dimension and higher order nonlinearity is considered.

1. Introduction

Orszag [1] proposed the Eddy Damped Quasi Normal Markovian (EDQNM) model as an efficient and realizable closure for three-dimensional (3D) isotropic turbulence. The first numerical study of the EDQNM was made by Leith [2], who formulated and solved it for two-dimensional (2D) isotropic turbulence. Herring [3] derived a statistical dynamical closure for 2D anisotropic turbulence that generalized the EDQNM to situations without transient waves. The EDQNM is computationally efficient because it is Markovian and consists of just the equations for the single-time two-point cumulant in spectral space with an analytical expression for the triad relaxation function. However, as we discuss and analyze in detail in this article, the EDQNM is not guaranteed to be realizable for anisotropic turbulence interacting with transient waves. In contrast, the Eddy Damped Markovian Anisotropic Closure (EDMAC), for which we formulate the theoretical basis, is realizable.

One can arrive at the EDQNM in various ways, as discussed, for example, by Lesieur [4], with our preference being from a reduction of the more fundamentally based Eulerian Direct Interaction Approximation (DIA) of Kraichnan [5]. The DIA is a more complex non-Markovian closure that is founded on renormalized perturbation theory (reviews of closures are presented in Refs. [4,6,7,8,9,10,11,12,13]). These reviews present analyses of the relative performance of the various closure models. The form of the eddy damping in the EDQNM is chosen to be consistent with the forward $k^{-\frac{5}{3}}$ energy cascading inertial range of 3D isotropic turbulence and the forward $k^{-3}$ enstrophy cascading inertial range of 2D turbulence.

Starting from the DIA, the EDQNM may be obtained by using a Markovian form for the response function and using the current-time FDT (Fluctuation Dissipation Theorem)

$$C_{\mathbf{k}}(t,t^{\prime })=R_{\mathbf{k}}(t,t^{\prime })C_{\mathbf{k}}(t,t)$$

(1)

for $t\ge t^{\prime }$ and $C_k(t,t^{\prime })=C_{-\mathbf{k}}(t^{\prime },t)$ for $t^{\prime }>t$. In Equation (1), $R_{\mathbf{k}}(t,t^{\prime })$ is the response function, $C_{\mathbf{k}}(t,t^{\prime })$ is the two-time spectral cumulant, and $C_{\mathbf{k}}(t,t)$ the current-time single-time cumulant at wavenumber (or strictly wavevector) $k$. In the EDQNM, the response function is also specified by an analytical form consistent with either the 2D or 3D inertial ranges. With these modifications, the time-history integrals of the DIA can be performed analytically to determine a triad relaxation function that enters the cumulant equation (see Section 7). In the EDQNM, the strength of eddy damping is specified by an empirical constant, as distinct from the DIA, which has no empirical parameters.

The EDQNM closure has the advantage of being very computationally efficient since it scales like $O(T)$, where $T$ is the length of the time integration. This contrasts with the DIA non-Markovian closure that scales as $O(T^3)$, because the two-point cumulants and response functions satisfy integro-differential equations. Although, we should note that the efficiency of the DIA can be improved to scaling like $O(T^2)$. This can be performed by using the three-point cumulant in periodic restarts of the DIA closure numerical model [14,15,16,17].

The DIA is closely related to the Self Consistent Field Theory closure (SCFT) of Herring [18,19] and the Local Energy-Transfer Theory closure (LET) of McComb [6,20,21]. The SCFT and LET closures are also non-Markovian and, despite being originally derived by different methods, can in fact also be obtained by modifications of the DIA equations. These three Eulerian non-Markovian closures all have the same single-time two-point closure equation. However, the SCFT and LET use the prior-time FDT [12] (Equation (1)) defined by

$$C_{\mathbf{k}}(t,t^{\prime })=R_{\mathbf{k}}(t,t^{\prime })C_{\mathbf{k}}(t^{\prime },t^{\prime })$$

(2)

for $t\ge t^{\prime }$ where $C_{\mathbf{k}}(t^{\prime },t^{\prime })$ is the prior-time single-time cumulant. In the SCFT, the response function, identical to that for the DIA, is taken as fundamental and the two-time cumulant is obtained from the prior-time FDT. For the LET the two-time cumulant, identical to that for the DIA, is fundamental, and the response function is derived from the FDT in Equation (2). In this study, we focus on the derivation of Markovian closures from the DIA closure, but we could equally consider the SCFT or LET as the starting points.

The EDQNM closure has been widely used in many studies and applications, as reviewed in Refs. [4,9,10,11,12,22]. As emphasized in the detailed study of Bowman et al. [23], the EDQNM may not be realizable in the presence of transient waves such as drift waves in plasmas or, equivalently, Rossby waves in geophysical flows. In some studies, the waves have been assumed to be slowly varying, and the quasi-steady state form of the triad relaxation function is used to avoid possible singular behavior [22,23,24,25,26,27,28,29]. In other studies, modified forms of quasi-normal Markovian closures have instead been used [30,31].

The reason that the EDQNM may not necessarily be realizable with transient waves is because the oscillations of the propagating waves mean that the real part of the triad relaxation functions in the EDQNM may, on occasions, become negative. Bowman et al. [23] found that the same problem also occurs if the prior time FDT in Equation (2) is used. However, they showed that a realizable Markovian closure could be obtained by using a FDT that is essentially a half-way house between the current-time and prior-time FDTs. This FDT, which we call the correlation FDT, is given by

$$C_{\mathbf{k}}(t,t^{\prime })={[C_{\mathbf{k}}(t,t)]}^{\frac{1}{2}}{R}_{\mathbf{k}}(t,t^{\prime }){[C_{\mathbf{k}}(t^{\prime },t^{\prime })]}^{\frac{1}{2}}$$

(3)

for $t\ge t^{\prime }$. Bowman et al. [23] called their closure the Realizable Markovian Closure (RMC). Unlike the EDQNM closure, for which the triad relaxation functions have analytical expressions, the RMC relaxation functions need to be evaluated by solving a system of auxiliary differential equations. Hence the RMC is less computationally efficient than the EDQNM. Hu et al. [32] used the RMC in studies of the Hasegawa–Wakatani two-field equation of plasma physics. Bowman and Krommes [33] further developed and applied realizable Markovian closures for homogeneous turbulence in the presence of waves. They considered the interaction of drift waves (essentially Rossby waves) and anisotropic turbulence with the single-level Charney–Hasegawa–Mima model. They also formulated a realizable test-field model (RTFM), a generalization of Kraichnan’s [34] test-field model (TFM), for turbulence in the presence of transient waves.

Traditionally, closure theory, including Markovian closures, like the EDQNM [1,2], have focused on the forward energy cascade of 3D turbulence and the forward enstrophy cascade of 2D turbulence in decaying and steady state systems. The EDQNM damping is normally specified to be consistent with these forward inertial ranges. Two-dimensional turbulence also has an inverse energy cascade, with a steady state inertial range of $k^{-\frac{5}{3}}$ [35]. With differential rotation on a $\beta$–plane, supporting Rossby waves, the turbulence becomes anisotropic, with preferential zonal motion. This was noted by Rhines [36] based on theory and direct numerical simulations (DNS) and by Holloway and Hendershott [37] in studies with a modified TFM. Their TFM did not cater for transient Rossby waves, but required an approximate steady state triad relaxation function to guarantee realizability [33]. Rhines [36] proposed that, under suitable conditions, the fluid velocities follow an inverse $k^{-5}$ power law in steady state where $k$ is the wavenumber, the magnitude of the total wavevector. Chekhlov et al. [38] performed a set of DNS experiments of 2D turbulence on a $\beta$–plane forced by small scale forcing. They found support for the scaling of Rhines, for the zonal components with $k_x=0$ which followed a $k_y^{-5}$ power law. However, the rest of the spectrum followed the $k^{-\frac{5}{3}}$ inverse energy cascade, emphasizing the anisotropy induced by the Rossby waves [22,24,25,26,27,28,29,37]. This is further discussed in the review by Galperin et al. [39]. Krommes and Parker [40] also discuss the application of homogeneous closures to circulations with zonal flow generation through inverse cascades, such as on the giant planets and their moons, and to plasma physics (see also the reviews in Refs. [41,42]). They note, however, that real flows are generally inhomogeneous and “statistically inhomogeneous solutions with nontrivial mean fields can emerge from a statistically homogeneous state by spontaneous symmetry breaking”. Indeed, Frederiksen [43] shows that just an infinitesimal topographic perturbation can cause such symmetry breaking with the mountain torque driving a mean westward zonal flow on the $\beta$–plane and solid body rotation on the sphere. This symmetry breaking may considerably modify the structures of the larger scale flows in the inverse cascades.

For the earth’s atmosphere, with a Rossby radius of deformation around 1000 km, the larger scales of turbulent 2D and quasigeostrophic flows are determined by heating, topography, and baroclinic instability that override the inverse cascades [4]. It is 2D homogeneous flows, in typical atmospheric parameter ranges, with forward enstrophy cascades, and the analogous 3D homogeneous flows with forward energy cascades, that are of primary interest in our current study. For flows in these parameter ranges, the large-scale flows can be specified to have realistic initial spectra or driven towards such states with the smaller scales evolving towards their statistically steady states. In this way, studies with homogeneous closures can throw considerable light onto the essential issues of the more complicated inhomogeneous turbulence interactions.

In the last few years, Markovian closures have also been developed for inhomogeneous turbulence interacting with mean flows, Rossby waves, and topography by Frederiksen and O’Kane [12,44,45]. Three types of Markovian Inhomogeneous Closures (MICs) were developed [44], as well as three types of abridged MICs [45]. For the abridged MICs, the mean field trajectory in the time history integrals is replaced by the current-time mean field and results in slightly more efficient closures. To date, all the MIC variants have been developed as modifications of the inhomogeneous non-Markovian QDIA closure [46,47,48].

The QDIA was developed for 2D inhomogeneous turbulent flows interacting with mean fields and topography by Frederiksen [46]. It was extended to include Rossby waves (Frederiksen and O’Kane [47]) and to multi-level and multi-field models for classical and quantum field theories (Frederiksen [48,49]). It was numerically implemented by O’Kane and Frederiksen [17] and subsequently used in further studies of turbulence dynamics, interactions with topography and Rossby waves, predictability of transitions between strong zonal flows and blocking, in data assimilation and for determining subgrid-scale parameterizations. The literature is reviewed in Refs. [45,49].

The MICs are obtained from the QDIA closure by assuming one of the three FDTs. They can be combined as:

$$C_{\mathbf{k}}(t,t^{\prime })={[C_{\mathbf{k}}(t,t)]}^{1-X}{R}_{\mathbf{k}}(t,t^{\prime }){[C_{\mathbf{k}}(t^{\prime },t^{\prime })]}^X$$

(4)

for $t\ge t^{\prime }$ and $C_k(t,t^{\prime })=C_{-\mathbf{k}}(t^{\prime },t)$ for $t^{\prime }>t$. The current-time FDT is recovered when $X=0$, the prior-time FDT has $X=1$ and the correlation FDT has $X=\frac{1}{2}$.

In the numerical studies of Frederiksen and O’Kane [44,45], it was found that both the MICs [44] and abridged MICs [45] closely reproduced the statistics of large ensembles of DNS in simulations of inhomogeneous turbulence interacting with mean flows, Rossby waves and topography at low Reynolds numbers. While the MIC and abridged MIC models are more computationally efficient than the non-Markovian QDIA closure, the relaxation functions still need to be calculated by solving auxiliary differential equations, as for the RMC of Bowman et al. [23]. A still more efficient inhomogeneous Markovian closure with analytical relaxation functions, the Eddy Damped Markovian Inhomogeneous Closure (EDMIC), was also developed by Frederiksen and O’Kane [45].

In a recent tribute to Jack Herring [12], and review of his major achievements in statistical dynamical fluid dynamics, it was noted that it is possible to generalize the EDQNM to a closure that is realizable in the presence of transient Rossby waves interacting with anisotropic turbulence. The resulting Eddy Damped Markovian Anisotropic Closure (EDMAC) involves a frequency renormalization of the eddy damping and by construction is as efficient as the EDQNM.

The major aims of this article are as follows:

• To provide a theoretical motivation for the realizable EDMAC based on renormalized perturbation theory.
• To generalize the EDMAC model for the interaction of anisotropic turbulence with transient Rossby waves in the presence of transient large-scale flows.
• To establish conditions under which the real part of the triad relaxation functions is positive semi-definite when the large-scale flow and Doppler shifted Rossby wave frequency have general time dependencies.
• To examine the extent to which the frequency-dependent contribution to the eddy damping in the EDMAC model changes the evolved energy and palinstrophy spectra and Reynolds number and skewness, compared with the EDQNM, for rapidly evolving moderate Reynolds number turbulence interacting with Rossby waves.

The methodology we apply for developing the realizable EDMAC model has parallels with the approach used by Frederiksen [46] to formulate the inhomogeneous quasi-diagonal direct interaction approximation (QDIA) closure. Frederiksen [46] developed the QDIA using a renormalized perturbation theory approach in which the isotropic problem was treated as the zero-order problem. The anisotropic and inhomogeneous terms were considered as small and multiplied by a perturbation parameter $\lambda$, prior to renormalization. Here we apply this approach to formulate the EDMAC.

The plan of this article is as follows. In Section 2, we present the equations for two-dimensional barotropic flows on a generalized $\beta$–plane and with large-scale advection by a mean flow. The dynamical flow equations are then converted to the equivalent form in Fourier space in Section 3. This is necessary for the subsequent formulation, in Section 4, of non-Markovian closures for the interaction of anisotropic turbulence with Rossby waves. These generalized DIA, SCFT, and LET closures have the same single-time cumulant equation and are equally suitable for reducing to the Markovian anisotropic closures (MACs) formulated in Section 5. There, three variants of the MACs are formulated using the three versions of the FDT that are summarized in Equation (4). As well, in Section 5, the prognostic equations for the relaxation functions that close the Markovian statistical dynamical equations are derived. In Section 6, we focus on the Markovian variant for which the current-time FDT in Equation (1) is imposed, and note the simplifications in the structure of the closure that then entail. The EDMAC model is formulated in Section 7 and in Section 8 the sufficient conditions for the realizability of EDMAC in the presence of transient Rossby waves are established. Generalization of the realizable EDMAC to other physical systems with transient waves and to higher dimension and higher order nonlinearity is considered in Section 9. In Section 10, numerical integrations of the EDMAC and EDQNM models are reported on for isotropic and anisotropic 2D turbulence of various complexity. The implications of the results and our conclusions are presented in Section 11. The derivation, based on renormalized perturbation theory, of the generalized DIA closure from which our analysis starts is located in Appendix A. The Langevin equation for the EDMAC model that guarantees realizability is presented in Appendix B.

2. Two-Dimensional Barotropic Flows on a $\mathit{\beta }$–Plane with Large-Scale Advection

We formulate our statistical dynamical equations for the case of two-dimensional or barotropic flows, although the generalization to three-dimensional flows can also be accomplished, as described in Section 9. Turbulent flows in planar geometry are considered on a generalized $\beta$–plane, and include a large-scale advection by a wind $U$. The total flow is described by the streamfunction $\mathrm{\Psi }=\psi -Uy$ where $\psi$ represents the small scales. Throughout this paper, we present theoretical and numerical results for flows on the doubly periodic plane $0\le x\le 2\pi , 0\le y\le 2\pi$ with $x=(x,y)$.

2.1. Large-Scale Flow Equation

The large-scale flow $U$ is kept the same in each realization of the smaller scales, and effectively just modifies the beta effect. We include possible relaxation to a mean large-scale flow $U_0$ with strength ${\alpha }_U$ so that $U$ may evolve with time according to

$$\frac{\partial U}{\partial t}={\alpha }_U(U_0-U).$$

(5)

In the absence of the relaxation forcing and dissipative drag, $U$ is separately conserved.

2.2. Barotropic Vorticity Equation for the Small Scales

The evolution of the small scales is described by the barotropic vorticity equation on a $\beta$–plane with the advective wind $U$ modifying the Rossby wave propagation:

$$\frac{\partial \zeta }{\partial t}=-J(\psi -Uy,\zeta +\beta y+k_0^{2}Uy)+{\nu }_0{\nabla }^2\zeta +f_0.$$

(6)

The Jacobian is defined by:

$$J(\psi ,\zeta )=\frac{\partial \psi }{\partial x}\frac{\partial \zeta }{\partial y}-\frac{\partial \psi }{\partial y}\frac{\partial \zeta }{\partial x}$$

(7)

and the relationship between the vorticity $\zeta$ and the stream function $\psi$ is

$$\zeta ={\nabla }^2\psi \equiv \left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)\psi$$

(8)

where ${\nabla }^2$ is the Laplacian.

The results that we present can also easily be extended to turbulent flow on a sphere, where $U$ is the analogue of the solid body rotation. Indeed, the structure of Equation (6) is the same as for flow on a sphere. The wavenumber $k_0$ is the analogue of that on the sphere with $k_0^{2}Uy$ corresponding to the vorticity of the solid body rotation. In Equation (6), ${\nu }_0$ is the viscosity, $\beta$ is the beta effect and $f_0$ specifies possible external forcing.

3. Dynamical Equations in Fourier Space

The statistical dynamics of turbulent flows is most conveniently analyzed in spectral space. For the doubly periodic domain, the equations are transformed into Fourier space by spectral representations of the fields in the form:

$$\zeta (\mathbf{x},t)={\displaystyle \sum _{\mathbf{k}}{\zeta }_{\mathbf{k}}(t)\exp (i \mathbf{k}\mathbf{.}\mathbf{x})},$$

(9)

where the spectral coefficients

$${\zeta }_{\mathbf{k}}(t)=\frac{1}{{(2\pi )}^2}{\displaystyle {\int }_0^{2\pi }{d}^2\mathbf{x}\zeta (\mathbf{x},t)\exp (-i \mathbf{k}\mathbf{.}\mathbf{x})}.$$

(10)

The wavenumbers in the $x$ and $y$ directions are $k_x$ and $k_y$ with the wavenumber vector $\mathbf{k}=(k_x,k_y)$ and the magnitude $k={(k_{\mathrm{x}}^2+k_{\mathrm{y}}^2)}^{1/2}$. The spectral coefficients satisfy the complex conjugate symmetry ${\zeta }_{-k}={\zeta }_k^{*}$ which guarantees that the fields in physical space are real. The domain in Equation (9) can be general, but we suppose that it contains the integer wavenumbers in a circular area (to be specified) that excludes the origin $0$. The vorticity equation in spectral space then takes the form:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t)\right){\zeta }_{\mathbf{k}}(t)\\ ={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q}){\zeta }_{-\mathbf{p}}(t){\zeta }_{-\mathbf{q}}(t}})+f_0(k,t).\end{array}$$

(11)

Here, we have generalized the form of the viscosity ${\nu }_0\to {\nu }_0(k)$ which corresponds to more general dissipation operators than the Laplacian in Equation (6). In Equation (11), the delta function $\delta (\mathbf{k},\mathbf{p},\mathbf{q})=1$ if $\mathbf{k}+\mathbf{p}+\mathbf{q}=0$ and 0 if $\mathbf{k}+\mathbf{p}+\mathbf{q}\ne 0$. The interaction coefficient $K(\mathbf{k},\mathbf{p},\mathbf{q})$ is defined by:

$$K(k,p,q)=\frac{1}{2}[p_x{q}_y-p_y{q}_x](p^2-q^2)/p^2{q}^2.$$

(12)

The Doppler shifted Rossby wave frequency is given by the dispersion relationship:

$${\omega }_k^U(t)={\mathrm{\Omega }}_k^U(t)+{\omega }_k^{\beta }=\frac{U(t)k_x(k^2-k_0^2)}{k^2}-\frac{\beta k_x}{k^2}$$

(13)

where the Rossby wave frequency is

$${\omega }_k\equiv {\omega }_k^{\beta }=-\frac{\beta k_x}{k^2},$$

(14)

and

$${\Omega }_k^U(t)=\frac{U(t)k_x(k^2-k_0^2)}{k^2}.$$

(15)

We note that in the case of inviscid flows Rossby waves of the form $\mathrm{expi}(k.x-{\omega }_k^{U}t)$ are exact solutions to Equation (6).

4. Non-Markovian Closures for Turbulence and Rossby Waves

The aim with closure theory is to develop statistical dynamical equations that describe the evolution of the low order moments or cumulants for infinite ensembles of flow fields. To start the process, we first represent a given member ${\zeta }_{\mathbf{k}}(t)$ by its mean, the one-point cumulant, $<{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t)$, and the deviation from the mean, ${\tilde{\zeta }}_{\mathbf{k}}(t)$. We consider homogeneous turbulence in this study for which the small scales have zero mean:

$$<{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t)=0;{\zeta }_{\mathbf{k}}(t)= {\tilde{\zeta }}_{\mathbf{k}}(t).$$

(16)

Since the mean is zero, the deviation ${\tilde{\zeta }}_{\mathbf{k}}$ also satisfies Equation (11) (with ${\zeta }_{\mathbf{k}}(t)\to {\tilde{\zeta }}_{\mathbf{k}}(t)$). The mean forcing also needs to be zero:

$$<f_0(k,t)>\equiv {\overline{f}}_0(k,t)=0;f_0(k,t)= {\tilde{f}}_0(k,t).$$

(17)

4.1. Generalized DIA Closure for Anisotropic Turbulence and Rossby Waves

To formulate the realizable EDMAC model equations, we start with the DIA closure that is slightly generalized as derived in Appendix A. For homogeneous anisotropic turbulence, the DIA closure describes the evolution of the two-time two-point cumulant

$$C_{\mathbf{k}}(t,t^{\prime })=<{\tilde{\zeta }}_k(t){\tilde{\zeta }}_{-k}(t^{\prime })>$$

(18)

and response function

$$R_{\mathbf{k}}(t,t^{\prime })=〈{\tilde{R}}_{\mathbf{k}}(t,t^{\prime })〉$$

(19)

which is the ensemble average of individual responses. The response function for an individual disturbance is

$${\tilde{R}}_{\mathbf{k}}(t,t^{\prime })=\frac{\delta {\tilde{\zeta }}_{\mathbf{k}}(t)}{\delta {\tilde{f}}_0(\mathbf{k},t^{\prime })}.$$

(20)

The response function ${\tilde{R}}_{\mathbf{k}}(t,t^{\prime })$ measures the change in the field $\delta {\tilde{\zeta }}_{\mathbf{k}}(t)$ due to an infinitesimal $\delta {\tilde{f}}_0(\mathbf{k},t^{\prime })$ perturbation in the forcing at an earlier time. Here, $\delta$ denotes the functional derivative.

The two-time two-point cumulant equation for the generalized DIA, as formulated in Appendix A, is given by

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t)\right)C_{\mathbf{k}}(t,t^{\prime })\\ +{\displaystyle \underset{t_0}{\overset{t}{\int }}ds ({\eta }_k(t,s)+{\pi }_k^{\omega }(t,s))C_{-k}(t^{\prime },s)}\\ ={\displaystyle \underset{t_0}{\overset{t^{\prime }}{\int }}ds (S_k(t,s)+P_k^{\omega }(t,s)+F_0(k,t,s))R_{-k}(t^{\prime },s)}.\end{array}$$

(21)

This is the case for $t>t^{\prime }$, while for $t^{\prime }>t$ we have $C_k(t,t^{\prime })=C_{-\mathbf{k}}(t^{\prime },t)$. Here, the last term in Equation (21) arises from

$$<{\tilde{f}}_0(k,t){\tilde{\zeta }}_{-\mathbf{k}}(t^{\prime })>={\displaystyle \underset{t_0}{\overset{t^{\prime }}{\int }}ds F_0(k,t,s)R_{-\mathbf{k}}(t^{\prime },s)}$$

(22)

where $t_0$ is the initial time and the bare noise

$$F_0(k,t,s)= <{\tilde{f}}_0(k,t){\tilde{f}}_0^{\ast }(k,s)>.$$

(23)

Also, the nonlinear damping and nonlinear noise are derived in Appendix A and given by:

$${\eta }_{\mathbf{k}}(t,s)=-4{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})R_{-\mathbf{p}}(t,s)C_{-\mathbf{q}}(t,s),}}$$

(24)

and

$$S_{\mathbf{k}}(t,s)=2{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{k},-\mathbf{p},-\mathbf{q})C_{-\mathbf{p}}(t,s)C_{-\mathbf{q}}(t,s)}}.$$

(25)

The expressions for ${\pi }_{\mathbf{k}}^{\omega }(t,s)$ and $P_{\mathbf{k}}^{\omega }(t,s)$ are also derived in Appendix A, and are given by:

$${\pi }_{\mathbf{k}}^{\omega }(t,s)=R_{\mathbf{k}}(t,s){\omega }_{\mathbf{k}}^U(t){\omega }_{\mathbf{k}}^U(s)$$

(26)

and

$$P_{\mathbf{k}}^{\omega }(t,s)=C_{\mathbf{k}}(t,s){\omega }_{\mathbf{k}}^U(t){\omega }_{\mathbf{k}}^U(s).$$

(27)

In a similar way, the response function equation can be derived as

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t)\right)R_{\mathbf{k}}(t,t^{\prime })\\ +{\displaystyle \underset{t^{\prime }}{\overset{t}{\int }}ds({\eta }_k(t,s)+{\pi }_k^{\omega }(t,s))R_k(s,t^{\prime })}=\delta (t-t^{\prime })\end{array}$$

(28)

for $t\ge t^{\prime }$, and the Dirac delta function means that $R_k(t,t)=1$.

The final equation needed for the DIA closure is that for the single-time two-point cumulant:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+2\mathrm{Re}({\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t))\right)C_{\mathbf{k}}(t,t)\\ +2\mathrm{Re}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds ({\eta }_k(t,s)+{\pi }_k^{\omega }(t,s))C_{-k}(t,s)}\\ =2\mathrm{Re}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds (S_k(t,s)+P_k^{\omega }(t,s)+F_0(k,t,s))R_{-k}(t,s)}.\end{array}$$

(29)

The system of DIA equations is started from the initial conditions $C_k(t_0,t_0)$ and $R_k(t,t)=1$. Equation (29) can also be simplified to the form:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+2{\nu }_0(k)k^2\right)C_k(t,t)+2\mathrm{Re}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds {\eta }_k(t,s)C_{-k}(t,s)}\\ =2\mathrm{Re}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds (S_k(t,s)+F_0(k,t,s))R_{-k}(t,s)}\end{array}$$

(30)

since the ${\pi }_k^{\omega }$ and $P_k^{\omega }$ terms cancel because

$$\mathrm{Re}\{{\pi }_k^{\omega }(t,s)C_{-k}(t,s)-P_k^{\omega }(t,s)R_{-k}(t,s)\}=0.$$

(31)

4.2. The Abridged DIA Closure for Anisotropic Turbulence and Rossby Waves

As a first step towards developing Markovian Anisotropic Closures (MACs) for the interaction of turbulence and Rossby waves, we present, in this subsection, an abridged generalized DIA closure. We consider the situation where the large-scale field $U(t)$ is slowly varying in the time history integrals:

$$U(s)\to U(t);{\omega }_k^U(s)\to {\omega }_k^U(t).$$

(32)

Thus, the abridged DIA closure equations are again given by Equations (21) to (25) and (28) to (30), but with Equations (26) and (27) replaced by

$${\pi }_k^{\omega }(t,s)\to {\pi }_k^{\underset{\_}{\omega }}(t,s)=R_k(t,s){[{\omega }_k^U(t)]}^2,$$

(33)

and

$$P_k^{\omega }(t,s)\to P_k^{\underset{\_}{\omega }}(t,s)=C_k(t,s){[{\omega }_k^U(t)]}^2,$$

(34)

where the Markov approximation in Equation (32) is denoted by the superscript $\underset{\_}{\omega }$.

4.3. Generalized SCFT and LET Closures for Ansotropic Turbulence and Rossby Waves

We can also arrive at generalized SCFT and LET closures from the generalized DIA closure in the usual way by invoking the prior-time FDT in Equation (2). For the generalized SCFT this FDT determines the two-time cumulant instead of Equation (21). For the generalized LET closure, the FDT determines the response function instead of Equation (28). For the corresponding abridged SCFT and LET closures, Equations (33) and (34) again replace Equations (26) and (27).

5. Statistical Dynamical Equations for Markovian Anisotropic Closures

General Formulation of Markovian Anisotropic Closures

As the next step towards formulating the realizable EDMAC model with analytical form for the relaxation function we outline the theory of associated Markovian Anisotropic Closures (MACs). The MACs are described by the single-time cumulant equation and auxiliary integral, or equivalent differential, equations for two relaxation functions. The MACS result from employing any of the three FDTs in Equation (4) and Markovianizing the generalized DIA in Section 4. As noted there the version with $X=0$ is the current-time FDT, $X=\frac{1}{2}$ the correlation FDT, and $X=1$ the prior-time FDT. The corresponding MACs are denoted by $MA{C}^X$. The MACs can be formulated using the general expressions for ${\pi }_k^{\omega }(t,s)$ and $P_k^{\omega }(t,s)$ in Equations (26) and (27). However, since our final aim is to simplify the MAC with $X=0$ to the realizable EDMAC we make the derivations from the abridged DIA with ${\pi }_k^{\underset{\_}{\omega }}(t,s)$ and $P_k^{\underset{\_}{\omega }}(t,s)$ given in Equations (33) and (34).

The single-time two-point cumulant equation for the abridged generalized DIA, (and as well for the corresponding SCFT and LET closure) that is also used for each of the MAC models can be written as:

$$\begin{array}{l}\frac{\partial }{\partial t}{C}_{\mathbf{k}}(t,t)+2 \mathrm{Re}[{\mathcal{N}}_{\eta }(k,t)+{\mathcal{N}}_{\underset{\_}{\omega }}(k,t)+{\mathcal{N}}_0(k)]\\ =2 \mathrm{Re}[{\mathcal{F}}_S(k,t)+{\mathcal{F}}_{\underset{\_}{\omega }}(k,t)+{\mathcal{F}}_0(k,t)].\end{array}$$

(35)

We note that $C_{\mathbf{k}}(t,t)=C_{-\mathbf{k}}(t,t)$ is real. The $\mathcal{F}(\mathbf{k},t)$ and $\mathcal{N}(\mathbf{k},t)$ terms can be written in the following convenient forms for subsequent Markovianization. Firstly,

$${\mathcal{F}}_S(k,t)=2{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{k},-\mathbf{p},-\mathbf{q})\Delta (}}-\mathbf{k},-\mathbf{p},-\mathbf{q})(t),$$

(36)

where

$$\Delta (-\mathbf{k},-\mathbf{p},-\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_{-k}(t,s)C_{-p}(t,s)C_{-q}(t,s).$$

(37)

As well,

$${\mathcal{F}}_{\underset{\_}{\omega }}(k,t)={[{\omega }_k^U(t)]}^2\Lambda (-k,k)(t)$$

(38)

with

$$\Lambda (-\mathbf{k},\mathbf{k})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_{-k}(t,s)C_k(t,s).$$

(39)

The bare forcing spectrum

$${\mathcal{F}}_0(k,t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{F}_0(k,t,s)R_{-k}(t,s).$$

(40)

In a similar way,

$${\mathcal{N}}_{\eta }(k,t)=-4{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})\Delta (}}-\mathbf{p},-\mathbf{q},-\mathbf{k})(t),$$

(41)

with

$${\mathcal{N}}_{\underset{\_}{\omega }}(k,t)={[{\omega }_k^U(t)]}^2\Lambda (k,-k)(t),$$

(42)

and

$${\mathcal{N}}_0(k,t)=({\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U)C_k(t,t)={\mathcal{D}}_0(k)C_k(t,t).$$

(43)

Here, the dissipation and Rossby wave dispersion term is

$${\mathcal{D}}_0(k,t)={\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t).$$

(44)

Next, the three versions of the FDT in Equation (4) are applied and this simplifies the nonlinear noise and nonlinear damping terms in Equations (36) and (41). Then, the time history integrals can then be expressed by relaxation functions ${\Theta }^X$ and ${\Psi }^X$ for $X=0,\frac{1}{2},1$. The integral representations for ${\Theta }^X$ and ${\Psi }^X$ can in turn be replaced by differential equations that augment Equation (35) for the single-time cumulant. This system is therefore Markovian with the variant with $X=\frac{1}{2}$ also guaranteed to be realizable since it is a generalization of the RMC model of Bowman et al. [23].

From Equations (36), (37) and (4) we have

$$\begin{array}{ll}{\mathcal{F}}_S(k,t)& =2{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{k},-\mathbf{p},-\mathbf{q})}}\\ & \times C_{-p}^{1-X}(t,t)C_{-q}^{1-X}(t,t){\Theta }^X(-\mathbf{k},-\mathbf{p},-\mathbf{q})(t),\end{array}$$

(45)

with the triad relaxation function

$${\Theta }^X(-\mathbf{k},-\mathbf{p},-\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_{-k}(t,s)R_{-p}(t,s)R_{-q}(t,s)C_{-p}^X(s,s)C_{-q}^X(s,s).$$

(46)

As well,

$${\mathcal{F}}_{\underset{\_}{\omega }}(k,t)={[{\omega }_k^U(t)]}^2{C}_k^{1-X}(t,t){\Psi }^X(-k,k)(t)$$

(47)

where the relaxation function

$${\Psi }^X(-\mathbf{k},\mathbf{k})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_{-k}(t,s)R_k(t,s)C_k^X(s,s).$$

(48)

In a similar way, ${\mathcal{N}}_{\eta }(k,t)$ and ${\mathcal{N}}_{\underset{\_}{\omega }}(k,t)$ simplify to

$${\mathcal{N}}_{\eta }(k,t)={\mathcal{D}}_{\eta }(k,t)C_k(t,t),$$

(49)

where

$$\begin{array}{l}{\mathcal{D}}_{\eta }(k,t)=-4{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})}}\\ \times C_{-q}^{1-X}(t,t)C_{-k}^{-X}(t,t){\Theta }^X(-\mathbf{p},-\mathbf{q},-\mathbf{k})(t),\end{array}$$

(50)

and

$${\mathcal{N}}_{\underset{\_}{\omega }}(k,t)={\mathcal{D}}_{\underset{\_}{\omega }}(k,t)C_k(t,t),$$

(51)

with

$${\mathcal{D}}_{\underset{\_}{\omega }}(k,t)={[{\omega }_k^U(t)]}^2{C}_k^{-X}(t,t){\Psi }^X(k,-k)(t).$$

(52)

The equation for the single-time cumulant can also be written in the form:

$$\left(\frac{\partial }{\partial t}+2 \mathrm{Re}({\mathcal{D}}_r(k,t))\right)C_{\mathbf{k}}(t,t)=2 \mathrm{Re}({\mathcal{F}}_r(k,t))$$

(53)

and for consistency the response function equation becomes:

$$\frac{\partial }{\partial t}{R}_{\mathbf{k}}(t,t^{\prime })+{\mathcal{D}}_r(k,t)R_{\mathbf{k}}(t,t^{\prime })=\delta (t-t^{\prime }).$$

(54)

The renormalized dissipation and noise functions in Equations (53) and (54) are defined by

$$\begin{array}{l}{\mathcal{D}}_r(k,t)={\mathcal{D}}_0(k)+{\mathcal{D}}_{\eta }(k,t)+{\mathcal{D}}_{\underset{\_}{\omega }}(k,t);\\ {\mathcal{F}}_r(k,t)={\mathcal{F}}_0(k,t)+{\mathcal{F}}_S(k,t)+{\mathcal{F}}_{\underset{\_}{\omega }}(k,t)\end{array}$$

(55)

with the contributing terms given above.

The modified form of the response function can also be written as:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t)\right)R_{\mathbf{k}}(t,t^{\prime })-\delta (t-t^{\prime })\\ =-\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}ds({\eta }_k^X(t,s)+{\pi }_k^{\underset{\_}{\omega }}(t,s))R_{-\mathbf{k}}(t,s){[C_{\mathbf{k}}(s,s)]}^X}{[C_k(t,t)]}^{-X}\right]R_k(t,t^{\prime })\end{array}$$

(56)

for $t\ge t^{\prime }$ with $R_k(t,t)=1$ and $C_k(t,t)=C_{-k}(t,t)$ is real. Here we have used the general FDT in Equation (4) and ${\eta }_{\mathbf{k}}^X$ is then also obtained from ${\eta }_{\mathbf{k}}$ in Equation (24) and has the expression:

$$\begin{array}{ll}{\eta }_{\mathbf{k}}^X(t,s)=& -4{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})}}\\ & \hspace{1em}\times R_{-\mathbf{p}}(t,s)R_{-\mathbf{q}}(t,s)C_{\mathbf{q}}^{1-X}(t,t)C_{\mathbf{q}}^X(s,s).\end{array}$$

(57)

The simplified form of the response function in Equation (54) together with the FDT in Equation (4) means that the relaxation functions ${\Theta }^X$ in Equation (46) and ${\Psi }^X$ in Equation (48) can be replaced by differential equations. The differential equation for ${\Theta }^X$ is:

$$\begin{array}{l}\frac{\partial }{\partial t}{\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(t)+({\mathcal{D}}_r(k,t)+{\mathcal{D}}_r(p,t)+{\mathcal{D}}_r(q,t)){\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(t)\\ =C_p^X(t,t)C_q^X(t,t)\end{array}$$

(58)

where ${\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(0)=0$ and ${\mathcal{D}}_k^r$ is given in Equation (55). The corresponding differential equation for ${\Psi }^X$ is:

$$\frac{\partial }{\partial t}{\Psi }^X(\mathbf{k},-\mathbf{k})(t)+({\mathcal{D}}_r(\mathbf{k},t)+{\mathcal{D}}_r(-\mathbf{k},t)){\Psi }^X(\mathbf{k},-\mathbf{k})(t)=C_{\mathbf{k}}^X(t,t)$$

(59)

with ${\Psi }^X(\mathbf{k},-\mathbf{k})(0)=0$.

Now, the closure for any of the three MACs with $X=0,\frac{1}{2},1$ is Equation (53) for the single-time cumulant $C_{\mathbf{k}}(t,t)$, augmented by the differential equations for ${\Theta }^X(\mathrm{k},p,q)(t)$ in Equation (58) and ${\Psi }^X(\mathbf{k},-\mathbf{k})(t)$ in Equation (59). The replacement of the integral forms for ${\Theta }^X$ and ${\Psi }^X$ by the differential equations has made the system Markovian. The auxiliary differential equations for ${\Theta }^X$ and ${\Psi }^X$ of course generate the same information as the time history integral forms but for long time simulations are more efficient. Nevertheless, analytical forms for ${\Theta }^X$ and ${\Psi }^X$ result in even more efficient closures, such as the EDQNM and EDMAC of Section 7. Each of the three $MA{C}^X$ formulated in this section is realizable for pure turbulence without transient wave phenomena. The MAC model with $X=\frac{1}{2}$ is also realizable when transient Rossby waves are present [23]. We note that, in fact,

$${\mathcal{D}}_{\underset{\_}{\omega }}(k,t)C_k(t,t)={\mathcal{F}}_{\underset{\_}{\omega }}(k,t)$$

(60)

so these terms can be dropped from Equation (53).

6. Markovian Anisotropic Closure with Current-Time FDT

In this Section, we consider the further simplifications of the $MA{C}^X$ models when $X=0$ so that the current-time FDT is used. This is a step towards the formulation of the EDMAC model, in the following Section, with analytical forms for the relaxation functions. When $X=0$, Equation (53), for the single-time cumulant, simplifies to:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+2{\nu }_0(k)k^2\right)C_{\mathrm{k}}(t,t)=2N_k(t)+2\mathrm{Re}{\mathcal{F}}_0(k,t)\\ =8{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{p},\mathbf{q},\mathbf{k})}}\mathrm{Re}{\Theta }^{X=0}(\mathbf{k},\mathbf{p},\mathbf{q})(t)\\ \hspace{1em}\times C_q(t,t)\left[C_k(t,t)-C_{\mathbf{p}}(t,t)\right]+2\mathrm{Re}{\mathcal{F}}_0(k,t),\end{array}$$

(61)

where, for white noise, $2\mathrm{Re}{\mathcal{F}}_0(k,t)=F_0(k,t,t)$. Here, we have used Equation (60) and the fact that the interaction coefficients have the properties that $K(\mathbf{k},\mathbf{p},\mathbf{q})=K(\mathbf{k},\mathbf{q},\mathbf{p})$, and $K(\mathbf{k},\mathbf{p},\mathbf{q})+K(\mathbf{p},\mathbf{q},\mathbf{k})+K(\mathbf{q},\mathbf{k},\mathbf{p})=0$. Also note that single-time cumulants are real. As well, to establish Equation (61), we have used the fact that the triad relaxation function ${\Theta }^{X=0}$ is symmetric in the three wavenumber indices:

$${\Theta }^{X=0}(\mathbf{k},\mathbf{p},\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_k(t,s)R_p(t,s)R_q(t,s)$$

(62)

as seen from Equation (46), when $X=0$.

Similarly, the expression for ${\Psi }^{X=0}$ simplifies to

$${\Psi }^{X=0}(\mathbf{k},-\mathbf{k})(t)={\Psi }^{X=0}(-\mathbf{k},\mathbf{k})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_{-k}(t,s)R_k(t,s).$$

(63)

The corresponding differential equation forms of ${\Theta }^{X=0}$ and ${\Psi }^{X=0}$ are given in Equations (58) and (59) with the simplification that $C^{X=0}\to 1$ on the right-hand sides. Also, when $X=0$, Equation (50) reduces to

$${\mathcal{D}}_{\eta }(k,t)=-4{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{p},\mathbf{q},\mathbf{k})}}{C}_q(t,t){\Theta }^{X=0}(\mathbf{k},\mathbf{p},\mathbf{q})(t)$$

(64)

where we have used the properties of the interaction coefficients noted above and the symmetry properties of the triad relaxation function. Also, Equation (52) becomes

$${\mathcal{D}}_{\underset{\_}{\omega }}(k,t)={[{\omega }_k^U(t)]}^2{\Psi }^{X=0}(k,-k)(t)={[{\omega }_k^U(t)]}^2{\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{R}_k(t,s)R_{-k}(t,s).$$

(65)

We also note that if the current-time FDT in Equation (1) is used (with $X=0$ in Equation (4)), then the expression for the response function simplifies to

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+\mathrm{i}{\omega }_k^U(t)\right)R_{\mathbf{k}}(t,t^{\prime })\\ +\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}ds({\eta }_k^{X=0}(t,s)+{[{\omega }_k^U(t)]}^2{R}_k(t,s))R_{-\mathbf{k}}(t,s)}\right]R_k(t,t^{\prime })\\ =\delta (t-t^{\prime })\end{array}$$

(66)

for $t\ge t^{\prime }$ with $R_k(t,t)=1$ and ${\eta }_k^X$ is defined in Equation (57). Here we note that the expression for ${\mathcal{D}}_{\eta }(k,t)$ in Equation (64) where $X=0$ can also be written as

$${\mathcal{D}}_{\eta }(k,t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds{\eta }_k^{X=0}(t,s)R_{-\mathbf{k}}(t,s)}$$

(67)

and ${\mathcal{D}}_{\underset{\_}{\omega }}(k,t)$ is given in Equation (52) with $C^{X=0}\to 1$. Thus, Equation (66) is of course just Equation (54) written out in detail when $X=0$.

For general $X=0,\frac{1}{2},1$, corresponding to the three FDTs, the solution to the response function differential Equation (54) (and Equation (66)) is

$$R_{\mathbf{k}}(t,t^{\prime })=\exp \left(-{\displaystyle {\int }_{t^{\prime }}^{t}ds{\mathcal{D}}_r(k,s)}\right)$$

(68)

where ${\mathcal{D}}_r$ is defined in Equation (55). Of course, when $X=0$, ${\mathcal{D}}_{\eta }$ and ${\mathcal{D}}_{\underset{\_}{\omega }}$ are given in Equations (64) and (65). We also note that, for $X=0$ the triad relaxation function in Equation (46) reduces to

$${\Theta }^{X=0}(\mathbf{k},\mathbf{p},\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}d{t}^{\prime }\exp \left(-{\displaystyle {\int }_{t^{\prime }}^{t}ds\left[{\mathcal{D}}_r(k,s)+{\mathcal{D}}_r(p,s)+{\mathcal{D}}_r(q,s)\right]}\right)}.$$

(69)

7. Realizable Eddy-Damped Markovian Anisotropic Closure

In this Section, we consider further simplifications that lead to the realizable EDMAC model. The approach involves the Markovian approximation for the damping terms, as in the EDQNM, and leads to analytical expressions for the relaxation functions.

7.1. Markovian Approximation and Analytical Eddy Damping

Thus, with the Markov approximation for ${\mathcal{D}}_r$, namely ${\mathcal{D}}_r(k,s)\to {\mathcal{D}}_r(k,t)$, the response function takes the simple form:

$$R_{\mathbf{k}}(t,t^{\prime })=\exp \left(-{\mathcal{D}}_r(k,t)(t-t^{\prime })\right),$$

(70)

since integral in Equation (68) can then be evaluated. Also, the integrals in the expression for the triad relaxation function ${\Theta }^{X=0}$ in Equation (69) can also be evaluated to give

$$\begin{array}{l}{\Theta }^{X=0}(\mathbf{k},\mathbf{p},\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}d{t}^{\prime }\exp \left(-\left[{\mathcal{D}}_r(k,t)+{\mathcal{D}}_r(p,t)+{\mathcal{D}}_r(q,t)\right](t-t^{\prime })\right)}\\ =\frac{1-\exp \left(-\left[{\mathcal{D}}_r(k,t)+{\mathcal{D}}_r(p,t)+{\mathcal{D}}_r(q,t)\right](t-t_0)\right)}{\left[{\mathcal{D}}_r(k,t)+{\mathcal{D}}_r(p,t)+{\mathcal{D}}_r(q,t)\right]}.\end{array}$$

(71)

For the EDMAC model, like the EDQNM closure, we use a parameterized form for the eddy damping ${\mathcal{D}}_{\eta }$, which appears in Equation (71) through ${\mathcal{D}}_r={\mathcal{D}}_0+{\mathcal{D}}_{\eta }+{\mathcal{D}}_{\underset{\_}{\omega }}$, and that is consistent with the $k^{-3}$ forward enstrophy cascading inertial range:

$${\mathcal{D}}_{\eta }(k,t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds{\eta }_k^{X=0}(t,s)R_{-\mathbf{k}}(t,s)}\to {\mu }_k^{eddy}(t)=\gamma {\left[k^2{C}_k(t,t)\right]}^{\frac{1}{2}}$$

(72)

where

$${\mu }_k(t)={\nu }_0(k)k^2+{\mu }_{\mathbf{k}}^{eddy}(t)={\nu }_0(k)k^2+\gamma {\left[k^2{C}_k(t,t)\right]}^{\frac{1}{2}}.$$

(73)

Here $\gamma$ is a positive dimensionless parameter. This local form for ${\mu }_k^{eddy}(t)$ is consistent with that of Leith [2] for 2D turbulence, and Orszag [1] used a corresponding local form for 3D turbulence applicable to the $k^{-\frac{5}{3}}$ forward energy cascading inertial range. Note, however, that our approach applies equally if an integral form over wavenumbers (Pouquet et al. [50]; Herring et al. [51]) is used instead for the eddy damping.

7.2. Frequency-Dependent Damping from Renormalized Perturbation Theory

Here, we consider the problem where the damping ${\mu }_k(t)$ is taken as the zero-order term and the Doppler-shifted frequency ${\omega }_k^U(t)\to \lambda {\omega }_k^U(t)$ is supposed to be a small perturbation of order $\lambda$. Then, we carry out the perturbation expansion as in Appendix A with the result that:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\mu }_k(t)\right)R_k(t,t^{\prime })-\delta (t-t^{\prime })\\ =-\mathrm{i}\lambda {\omega }_k^U(t)R_k^{(0)}(t,t^{\prime })-{\lambda }^2{\omega }_k^U(t)\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}ds{\omega }_k^U(s)R_k^{(0)}(t,s)R_{-k}^{(0)}(t,s)}\right]R_k^{(0)}(t,t^{\prime })\\ \approx -\mathrm{i}\lambda {\omega }_k^U(t)R_k^{(0)}(t,t^{\prime })-{\lambda }^2{[{\omega }_k^U(t)]}^2\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}ds{R}_k^{(0)}(t,s)R_{-k}^{(0)}(t,s)}\right]R_k^{(0)}(t,t^{\prime }).\end{array}$$

(74)

Here, $R_k^{(0)}(t,t^{\prime })=\exp -{\mu }_k(t)(t-t^{\prime })$ is the zero-order response function, as in Appendix A. We have assumed that $\mu$ and $\omega$ are slowly varying, and then we can evaluate the integral to obtain:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\mu }_k(t)\right)R_k(t,t^{\prime })-\delta (t-t^{\prime })\\ \approx -\mathrm{i}\lambda {\omega }_k^U(t)R_k^{(0)}(t,t^{\prime })-{\lambda }^2\frac{{[{\omega }_k^U(t)]}^2}{2{\mu }_k(t)}\left[1-\exp -(2{\mu }_k(t)t)\right]R_k^{(0)}(t,t^{\prime }).\end{array}$$

(75)

In formulating the EDQNM model, ${\mu }_k^{eddy}(t)$ is specified to have the form in Equation (72) (or the integral form of Pouquet et al. [50]) from the initial conditions rather than developing from zero. In the same way, we take the quasi-steady state expression for the contribution from ${[{\omega }_k^U(t)]}^2$. The result is:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\mu }_k(t)\right)R_k(t,t^{\prime })\\ \approx -\mathrm{i}\lambda {\omega }_k^U(t)R_k^{(0)}(t,t^{\prime })-{\lambda }^2\frac{{[{\omega }_k^U(t)]}^2}{2{\mu }_k(t)}{R}_k^{(0)}(t,t^{\prime })+\delta (t-t^{\prime }).\end{array}$$

(76)

Renormalizing, we obtain our required expression for establishing the realizable EDMAC model:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\mu }_k(t)\right)R_k(t,t^{\prime })\\ \approx -\mathrm{i}{\omega }_k^U(t)R_k(t,t^{\prime })-\frac{{[{\omega }_k^U(t)]}^2}{2{\mu }_k(t)}{R}_k(t,t^{\prime })+\delta (t-t^{\prime }).\end{array}$$

(77)

Thus, comparing the frequency-dependent term ${\mathcal{D}}_{\underset{\_}{\omega }}$, for the MAC model, in Equation (65) with Equation (77), we see that

$${\mathcal{D}}_{\underset{\_}{\omega }}(k,t)\to c\frac{{[{\omega }_k^U(t)]}^2}{{\mu }_k(t)}$$

(78)

where the above analysis indicates that $c=\frac{1}{2}$. We leave the constant $c$ in Equation (78) and in the following analysis, since it could be used as an empirical parameter in situations when the wave frequency is not small compared with the eddy damping. From Equation (52), with $C^{X=0}\to 1$, and Equation (78), we see that

$${\Psi }^{EDMAC}(\mathbf{k},-\mathbf{k})(t)=\frac{c}{{\mu }_k(t)}$$

(79)

where ${\mu }_k(t)$ is given in Equation (73).

Combining the above results with those of Section 7.1, we find

$${\mathcal{D}}_r(k,t)\to {\rho }_k(t)+\mathrm{i}{\omega }_k^U(t)$$

(80)

where

$${\rho }_k(t)={\mu }_k(t)+c\frac{{[{\omega }_k^U(t)]}^2}{{\mu }_k(t)}>0$$

(81)

is the frequency renormalized damping. Thus, the EDMAC response function becomes:

$$R_{\mathbf{k}}(t,t^{\prime })\approx R_k^{EDMAC}(t,t^{\prime })=\exp \left(-[{\rho }_k(t)+\mathrm{i}{\omega }_k^U(t)](t-t^{\prime })\right),$$

(82)

and the EDMAC triad relaxation function:

$$\begin{array}{l}{\Theta }^{EDMAC}(\mathbf{k},\mathbf{p},\mathbf{q})(t)\\ =\frac{1-\exp \left(-\left[{\rho }_k(t)+{\rho }_p(t)+{\rho }_q(t)+\mathrm{i}({\omega }_k^U(t)+{\omega }_p^U(t)+{\omega }_q^U(t))\right](t-t_0)\right)}{{\rho }_k(t)+{\rho }_p(t)+{\rho }_q(t)+\mathrm{i}({\omega }_k^U(t)+{\omega }_p^U(t)+{\omega }_q^U(t))}.\end{array}$$

(83)

The EDMAC model is then defined by Equation (61), but with ${\Theta }^{EDMAC}$ in Equation (83) replacing ${\Theta }^{X=0}$. In Section 8, to follow, and Appendix B, we determine sufficient conditions for the EDMAC to be realizable. The EDMAC has the same structure as the EDQNM, but with just ${\rho }_k$ replacing ${\mu }_k$ in the triad relaxation function. It is essentially as computationally efficient as the EDQNM, as discussed further in Section 10, but with the advantage of guaranteed realizability.

8. Conditions for Realizability of EDMAC with Variable Rossby Wave Frequency

In this Section, we determine sufficient conditions on the damping ${\rho }_k(t)>0$, in Equation (81), so that $\mathrm{Re}{\Theta }^{EDMAC}\ge 0$; that is, so that the real part of the triad relaxation in Equation (83) is a positive semi-definite. We consider the general case where the Doppler shifted frequencies ${\omega }_k^U(t)$ are time-dependent and are not necessarily monotonic.

The results of Section 7 suggest that $c=\frac{1}{2}$ is an appropriate value for small ${\omega }_k^U(t)$. More generally, for larger ${\omega }_k^U(t)$, the constant $c$ may be regarded as an empirical factor that is specified. Our aim is therefore to determine sufficient conditions on $c$ for realizability of the EDMAC model.

The triad relaxation function, ${\Theta }^{EDMAC}(t)={\Theta }^{EDMAC}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)(t)$, in Equation (83) for the EDMAC model, can be written in simplified form as:

$${\Theta }^{EDMAC}(t)=\frac{1-\exp -[\rho +\mathrm{i}\omega ](t-t_0)}{\rho +\mathrm{i}\omega }$$

(84)

where we take $t_0=0$ (without loss of generality). In Equation (84),

$$\begin{array}{l}\mu (t)={\displaystyle \sum _{j=1}^N{\mu }_{k_j}(t)}>0,\\ \omega (t)={\displaystyle \sum _{j=1}^N{\omega }_{k_j}^U(t)},\\ \rho (t)={\displaystyle \sum _{j=1}^N{\rho }_{k_j}(t)}={\displaystyle \sum _{j=1}^N\left({\mu }_{k_j}(t)+c\frac{{[{\omega }_{k_j}^U(t)]}^2}{{\mu }_{k_j}(t)}\right)}>0,\end{array}$$

(85)

with $N=3$ for triad interactions. The real part of the relaxation function that appears in the EDMAC model for the two-point cumulant (Equations (53) and (61) with superscript $X=0$ replaced by $EDMAC$) is then given by:

$$\begin{array}{l}\mathrm{Re}{\Theta }^{EDMAC}(t)\\ =\frac{1}{{\rho }^2+{\omega }^2}\left[\rho \{1-[\exp -\rho t]\cos \omega t\}+\omega [\exp -\rho t]\sin \omega t\right].\end{array}$$

(86)

Then, sufficient conditions for $\mathrm{Re}{\Theta }^{EDMAC}(t)\ge 0$ are determined as follows. In Equation (86), we can replace $\omega$ by $\left|\omega \right|$ for $t\ge 0$, since both $\cos \omega t$ and $\omega \sin \omega t$ are even functions of $\omega$. On the basis of the analysis in Section 7.2, we suppose that

$$\rho \ge \left|\omega \right|.$$

(87)

Firstly, we see that $\mathrm{Re}{\Theta }^{EDMAC}(t)\ge 0$, for $0<|\omega (t)|t\le \pi$, since then

$$1-[\exp -\rho t]\cos \left|\omega \right|t>0;\left|\omega \right|[\exp -\rho t]\sin \left|\omega \right|t\ge 0.$$

(88)

Secondly, from Equation (87), we see that for $|\omega (t)|t>\pi$,

$$\exp -\rho t\le \exp -(|\omega |t)<\exp -(\pi )<0.05.$$

(89)

Further, in general,

$$-1\le \cos \left|\omega \right|t\le 1;-1\le \sin \left|\omega \right|t\le 1$$

(90)

and thus, for $|\omega (t)|t>\pi$,

$$\left[{\rho }^2+{\omega }^2\right]\mathrm{Re}{\Theta }^{EDMAC}(t)>0.9\rho >0.$$

(91)

This then means that $\mathrm{Re}{\Theta }^{EDMAC}(t)\ge 0$, for $0\le |\omega (t)|t<\infty$. Sufficient conditions for the validity of Equation (87) is that each wavevector component satisfies

$${\rho }_{k_j}(t)=\left({\mu }_{k_j}(t)+c\frac{{[{\omega }_{k_j}^U(t)]}^2}{{\mu }_{k_j}(t)}\right)\ge |{\omega }_{k_j}^U(t)|.$$

(92)

In turn, the inequalities in Equation (92) are valid provided $c\ge \frac{1}{4}$, as seen by solving quadratic equations.

9. Generalizations of the EDMAC Model

The analysis of Section 8, determining sufficient conditions for $\mathrm{Re}{\Theta }^{EDMAC}\ge 0$, can in fact be generalized to all dimensions $d\ge 2$ and with $N\ge 3$ components interacting instead of just the three. The condition for this is again that $c\ge \frac{1}{4}$ since none of the argument in Section 8 depends on the dimension or the number of interacting components. It just requires that ${\mu }_k(t)>0$ and ${\omega }_k^U(t)$, with a general dispersion relationship, is finite.

This again has important implications for general EDQNM-type closures for which realizability is dependent on the relaxation function having positive semi-definite real part. Suppose that the relaxation function ${\Theta }^{EDQNM}(t)$ is given by the right hand side of Equation (84), but with $\mu$ replacing $\rho$. Then, including the effects of frequency renormalized damping, $\mu \to \rho$, as specified in Equation (85), results in the associated realizable EDMAC model.

It is straightforward to extend the EDMAC model to 2D turbulent flows and Rossby waves on a differentially rotating sphere [9,27]. It is also a simple matter to extend the EDMAC to 3D quasigeostrophic turbulent flows and waves in the model in Appendix B of Frederiksen [48] generalized to rotating flows on a $\beta$–plane and on a sphere. The EDQNM closure of Carnevale and Frederiksen [25] for two-dimensional internal gravity waves used the quasi-steady state form for the triad relaxation function to ensure realizability. This restriction can be removed, and transient internal gravity waves considered, by going to the EDMAC model form for the relaxation functions. Realizability can be assured by using the frequency renormalized form of the damping in Equation (81) where ${\omega }_k^U(t)$ is replaced by the internal wave frequency ${\omega }_k^s=s{k}_x/k$ with $s=\pm 1$ (Equation (2.7c) of Ref. [25]).

One might also expect to be able to develop suitable generalizations of the EDMAC to 3D Navier–Stokes turbulence [1,10,22], including with rotation and waves. In particular, the study by Cambon and Jacqui [52] considered 3D anisotropic turbulence subject to rotation with different Rossby numbers and waves in EDQNM type models. They also used the quasi-steady state form for the triad relaxation function that ensures realizability. This might be extended to the transient regime form by using a frequency renormalized damping as for the EDMAC model. Rose and Sulem [53] and Clark et al. [54] formulate the EDQNM closure for isotropic turbulence in general $d\ge 2$ dimensions, and it would be of interest to see if this could be generalized to anisotropic turbulence, with preferred zonal motion [27,36,37], subject to differential rotation and waves with the realizable EDMAC.

There are, of course, also important physical systems for which the nonlinearity is of a higher order than quadratic. Two examples, with cubic nonlinearity in the field equations, are the nonlinear Schrodinger equation (Nazarenko [55] reviews the literature) and the Klein Gordon equation with $\lambda {\varphi }^4$ Lagrangian (Frederiksen [49] reviews the literature). The nonlinear Schrodinger equation is of major importance in the study of optical turbulence, plasma physics and continuum mechanics [55,56] and the Klein Gordon equation in the classical and quantum statistical dynamics of Bose–Einstein condensation and scalar field theory [49,57,58]. In spectral space, $N=4$ in Equation (85), corresponding to quartic interactions, would then be needed in closures including EDMAC, EDQNM, and resonant interaction or wave turbulence models [25,49,54,55,56,57,58] of such phenomena.

10. Comparison of Closure Integrations for Turbulence and Rossby Wave Dynamics

In this Section, we compare closure calculations with the EDQNM and EDMAC models for isotropic turbulence and anisotropic turbulence with Rossby waves and possibly a large-scale mean flow. The integrations start from the isotropic spectrum B that was studied by Frederiksen and Davies [15,16] for closures with discrete spectra and which is very similar to the continuous spectrum II of Herring et al. [59]. These spectra were used by Herring et al. [59] and Frederiksen and Davies [15] in studies of the comparison of closures with DNS for 2D isotropic turbulence. In particular, the continuous wavenumber formulation of the DIA closure was compared with discrete wavenumber DNS by Herring et al. [59], while Frederiksen and Davies [15] used the discrete formulation for both. In both studies, the DIA closure, and its associated SCFT and LET closures in Ref. [15], were found to underestimate the small-scale amplitudes (Figures 22 and 23 of Ref. [59]; Figures 3 and 4 of Ref. [15]). However, the discrete closures [15] were found to be in better agreement with the statistics of DNS than the continuous closures [59]. O’Kane and Frederiksen [17] also used spectrum B [15,16] for the initial transient spectrum in studies of inhomogeneous turbulence interacting with mean flows and topography with the QDIA closure.

A regularized version of the DIA, the RDIA, in which the wavenumber ranges of interaction are restricted in the response function and two-time cumulant [16], was found to give quite close agreement with the statistics of DNS at the expense of an empirical parameter $\alpha$. This is shown in Figures 1 and 2 of Ref. [16], although there does seem to be a tendency of the smallest scales of the DNS to kick up somewhat in all these studies [15,16,59].

The EDQNM and EDMAC models also depend on the empirical parameter $\gamma$ in Equation (72). For the EDMAC model, we have estimated $c=\frac{1}{2}$ to be the strength of the frequency-dependent damping in Section 7.2, and we use that value. More generally, if the waves do not satisfy the formal requirement of being of small amplitude, then $c$ could also be regarded as an empirical parameter. Here, our aim is first to compare the EDQNM and EDMAC models for a commonly used value of $\gamma =0.6$ [2,9] for 2D turbulence and then to discuss the generality of our findings for a range of $\gamma$.

10.1. Diagnostics

The evolved closure integrations for the EDQNM and EDMAC models are compared in terms of the similarity of their evolved kinetic energy and palinstrophy spectra and Reynolds number and skewness. These diagnostics are defined as follows. Firstly, we specify the energy,

$$\mathcal{E}=\frac{1}{2}{\displaystyle \sum _k{C}_k(t,t)}{k}^{-2},$$

(93)

enstrophy,

$$\mathcal{F}=\frac{1}{2}{\displaystyle \sum _k{C}_k(t,t)},$$

(94)

palinstrophy,

$$\mathcal{P}=\frac{1}{2}{\displaystyle \sum _k{k}^2{C}_k(t,t)},$$

(95)

enstrophy dissipation,

$$\eta ={\displaystyle \sum _k{\nu }_0{k}^2{C}_k(t,t)=2}{\nu }_0\mathcal{P},$$

(96)

and palintrophy production,

$$\mathcal{K}={\displaystyle \sum _k{k}^2{N}_k(t)}.$$

(97)

Here, $N_k(t)$ is defined in Equation (61). The large-scale Reynolds number and skewness are then specified by

$$R_L(t)=\mathcal{E}/({\nu }_0{\eta }^{\frac{1}{3}})$$

(98)

and

$$S(t)=\mathcal{2}\mathcal{K}/(\mathcal{P}{\mathcal{F}}^{\frac{1}{2}}).$$

(99)

The transient kinetic energy spectra, and palinstrophy spectra, averaged over circular bands, are defined by:

$$E(k_i,t)=\frac{1}{2}{\displaystyle \sum _{k\in \mathcal{S}}{C}_k(t,t)k^{-2},}$$

(100)

and

$$P(k_i,t)=\frac{1}{2}{\displaystyle \sum _{k\in \mathcal{S}}{k}^2{C}_k(t,t),}$$

(101)

with the set $\mathcal{S}$ is defined by

$$\mathcal{S}=[k|k_i=\mathrm{Int}.(\left|k\right|+\frac{1}{2})].$$

(102)

The band-average is over all $k$ within a band of unit width at a radius $k_i$.

10.2. Initial Spectra and Parameter Specifications

In this subsection, we specify the parameters and initial spectra used in our closure calculations. They are based on the related studies in Refs. [9,15,16,47,59]. The length scale used is of one half the earth’s radius, $a_e/2$, and time scale is the inverse of the earth’s rotation rate, ${\Omega }^{-1}$. The the $\beta$–effect is zero in the isotropic integrations and $1.15\times {10}^{-11}{\mathrm{m}}^{-1}{\mathrm{s}}^{-1}$ (corresponding to $\beta =\frac{1}{2}$ in non-dimensional units) in the anisotropic runs. This is characteristic of the earth’s differential rotation at ${60}^{\mathrm{o}}$ latitude, and the parameter $k_0^2=\frac{1}{2}$ in Equation (6). We use a (bare) viscosity coefficient ${\nu }_0$ of $1.85\times {10}^6{\mathrm{m}}^2{\mathrm{s}}^{-1}$(non-dimensional ${\nu }_0=2.5\times {10}^{-3}$). In the closure integrations, the mean eastward flow $U$ is either zero or has a speed of $15{\mathrm{m}\mathrm{s}}^{-1}$ (non-dimensional $U=0.065$). The closure runs are all unforced, with $f_0=0$ in Equation (6) and the drag on the large-scale flow ${\alpha }_U=0$ in Equation (5). The parameter $\gamma$ in Equation (72) that specifies the strength of the eddy damping is specified at $\gamma =0.6$ as used in [9]. We have, however, checked that the broad conclusions regarding the similarity of the closure calculations described in Section 10.3 to follow hold for a wide range of $\gamma$ down to 0.01. The closure calculations are performed at a resolution of circular truncation C64 where $\left|k\right|=k\le 64$. All integrations proceed to $t_{\max }=90\mathrm{s}$ (non-dimensional $t_{\max }=0.4$) with a time step of $\Delta t=0.9\mathrm{s}$ (non-dimensional $\Delta t=0.004$). The time stepping is performed with a predictor-corrector scheme.

The closure calculations start from Gaussian initial conditions with $C_k(t_0,t_0)=C_k(0,0)$ for spectrum B of Frederiksen and Davies [15,16]:

$$C_k(0,0)=0.18k^2\exp \left(-\frac{2}{3}k\right).$$

(103)

From this spectrum, the closure integrations proceed with one run being for isotropic turbulence, denoted run $\mathcal{I}$, and three runs for anisotropic turbulence with Rossby waves, denoted runs $\mathcal{A}\mathcal{1}$ to $\mathcal{A}\mathcal{3}$. The particular values of nondimensional $\beta$, $c$, of Equation (78), and $U$ that are used in runs $\mathcal{I}$ and $\mathcal{A}\mathcal{1}$ to $\mathcal{A}\mathcal{3}$ are given in

Table 1. Non-dimensional parameters specifying the isotropic and three anisotropic closure runs.

Closure Runs	$\mathit{\beta }$	$\mathit{c}$	$\mathit{U}$
Isotropic Run  $\mathcal{I}$	0	0	0
Anisotropic Run  $\mathcal{A}\mathcal{1}$	0.5	0	0
Anisotropic Run  $\mathcal{A}\mathcal{2}$	0.5	0.5	0
Anisotropic Run  $\mathcal{A}\mathcal{3}$	0.5	0.5	0.065

.

10.3. Moderate Reynolds Number Closure Integrations

Next, we consider the evolution of the statistics of turbulence, possibly interacting with Rossby waves and a large-scale mean flow $U$ within the EDQNM and EDMAC models described in Section 7. For the isotropic run $\mathcal{I}$ and the anisotropic run $\mathcal{A}\mathcal{1}$, $c=0$, $U=0$, and the EDQNM closure has been used for these. The EDQNM is not guaranteed to be realizable for the anisotropic run with $\beta =0.5$. However, the EDQNM remained stable throughout the integration. Indeed, we found little difference in the results for the anisotropic run $\mathcal{A}\mathcal{1}$ with the EDQNM and the anisotropic run $\mathcal{A}\mathcal{2}$ with $c=0.5$ performed with the EDMAC model, which is guaranteed to be realizable in the presence of Rossby waves. This is shown in

Table 2. Initial and evolved Reynolds number and skewness for the isotropic and three anisotropic closure runs.

Closure Runs	${\mathit{R}}_{\mathit{L}}(0)$	${\mathit{R}}_{\mathit{L}}(0.4)$	$\mathit{S}(0)$	$\mathit{S}(0.4)$
Isotropic Run  $\mathcal{I}$	304.8	263.7	0	0.735
Anisotropic Run  $\mathcal{A}\mathcal{1}$	304.8	263.8	0	0.734
Anisotropic Run  $\mathcal{A}\mathcal{2}$	304.8	263.8	0	0.735
Anisotropic Run  $\mathcal{A}\mathcal{3}$	304.8	265.7	0	0.690

for the large-scale Reynolds number, and for the skewness, which are identical or essentially the same for runs $\mathcal{I}$, $\mathcal{A}\mathcal{1}$, and $\mathcal{A}\mathcal{2}$. The close similarity between these runs for the evolved transient kinetic energy and palinstrophy is also seen from Figure 1 and Figure 2. In fact, the results are so close that we have had to shift down the plots by increasing factors of 10 in order separate the results.

The inclusion of the large-scale mean flow $U=0.065$, corresponding to a dimensional value of $U=15{\mathrm{ms}}^{-1}$, in the anisotropic run $\mathcal{A}\mathcal{3}$ with the EDMAC model increases the evolved Reynolds number slightly and decreases the evolved skewness more. Thus, Rossby waves, in the presence of the large-scale flow, make the turbulence more Gaussian.

Thus, in these closure integrations, we conclude that the addition of the frequency-dependent damping in the EDMAC model does not greatly change the results compared with the EDQNM. This is the case for the strength of the eddy damping $\gamma =0.6$ used in the displayed results. We have also checked that applies equally for strengths of $\gamma$ down to 0.01. As $\gamma$ is increased to larger values, the frequency-dependent contribution to the eddy damping becomes smaller and the damping experienced by the EDMAC is closer to that of the EDQNM. While the EDQNM closure remained stable for the anisotropic run $\mathcal{A}\mathcal{1}$, this cannot always be guaranteed for the EDQNM, but it can for the EDMAC model under the conditions determined in Section 8 and Appendix B.

11. Discussion and Conclusions

The theoretical development of the Eddy Damped Markovian Anisotropic Closure (EDMAC) has been presented for anisotropic turbulence interacting with Rossby waves in the presence of advection by a large-scale mean flow. The EDMAC generalizes the Eddy Damped Quasi Normal Markovian (EDQNM) to a form that is realizable in the presence of transient waves. We have documented the equations for two-dimensional turbulence interacting with Rossby waves and a large-scale flow in physical space and in discrete Fourier spectral space.

11.1. Theoretical Results

The development of the EDMAC model has then proceeded in a number of steps requiring the formulation of generalized non-Markovian and Markovian closures.

11.1.1. Generalized Non-Markovian Closures

Firstly, renormalized perturbation theory has been employed to construct a generalization of the non-Markovian Direct Interaction Approximation (DIA) closure with additional eddy damping and nonlinear noise terms that depend on the Doppler shifted Rossby wave frequencies. Associated generalized Self Consistent Field Theory (SCFT) and Local Energy-Transfer Theory (LET) closures have also been introduced. The next step has been to slightly simplify these non-Markovian closures to abridged forms in which the large-scale mean flow is regarded as slowly varying in the time history integrals.

11.1.2. Markovian Anisotropic Closures

Markovian Anisotropic Closures (MACs) have then been developed by employing any of three forms of the Fluctuation Dissipation Theorem in Equation (4). This simplifies the time history integrals to the extent that their effects can be encapsulated in two relaxation functions that can also be determined by ordinary Markovian differential equations. While the MACs are more efficient for long time integrations than the non-Markovian closures, the calculation of the relaxation functions through differential equations is a considerable overhead compared with the analytical forms that are used for the EDQNM and established here for the EDMAC.

11.1.3. Realizable Eddy Damped Markovian Anisotropic Closure

The EDMAC has been derived from the MAC employing the current-time Fluctuation Dissipation Theorem (FDT) in Equation (1). Further simplification of the relaxation functions has been achieved by making the Markov approximation that the generalized dissipation functionals are slowly varying in the time history integrals, as for the EDQNM [1]. The consequent EDMAC model has the same eddy damping term as in the EDQNM that parameterizes turbulent mixing [1,2,50,51], but also has a frequency-dependent damping, which ensures the realizability of the EDMAC in the presence of transient waves. The strength of the frequency dependent damping has been determined under the assumption of small amplitude waves, again based on renormalized perturbation theory. More generally, the strength may be specified empirically, and sufficient conditions on the frequency-dependent damping that ensure realizability of the EDMAC have been determined.

11.1.4. Generalizations of the Eddy Damped Markovian Anisotropic Closure

We have considered the relationships between EDQNM and EDMAC models for anisotropic turbulence interacting with transient waves for systems involving higher order nonlinearity than quadratic and in higher dimension than two. The conditions that ensure that the real part of the EDMAC relaxation function is again positive semi-definite have been determined.

11.2. Numerical Closure Calculations

In this largely theoretical paper, we have also reported in detail four integrations of the closures for rapidly evolving turbulence. They have been performed for isotropic 2D turbulence, for anisotropic turbulence with Rossby waves with and without the frequency-dependent damping, and for anisotropic turbulence with Rossby waves and large-scale mean flow and frequency-dependent damping. In these numerical experiments, with empirical parameter $\gamma =0.6$ in Equation (72), the results of the evolved simulations with EDQNM and EDMAC turn out very similar to each other with little effect on one-dimensional spectra. This is also so for smaller $\gamma$ down to $\gamma =0.01$, while for larger $\gamma$, the frequency-dependent drain becomes smaller and the total drain in the EDMAC approaches that of the EDQNM.

11.3. Conclusions and Future Prospects

The realizability of the EDMAC in the presence of waves at little or no extra computational cost over the EDQNM should make it attractive provided the additional frequency-dependent damping does not change the broad properties of the numerical integrations with the closures. We have noted, in Section 9, the prospects for generalizing the realizable EDMAC model to higher dimensions and different physical systems, including with higher-order nonlinearity. In future studies, we plan to explore the performance of the EDMAC model in various settings. Also, we aim to formulate and study the performance of a similar inhomogeneous closure, the Eddy Damped Markovian Inhomogeneous Closure (EDMIC) [45], generalized to be realizable for turbulence interacting with transient Rossby waves.