Statistical Dynamics of Mean Flows Interacting with Rossby Waves, Turbulence, and Topography

Abstract

Abridged statistical dynamical closures, for the interaction of two-dimensional inhomogeneous turbulent flows with topography and Rossby waves on a beta–plane, are formulated from the Quasi-diagonal Direct Interaction Approximation (QDIA) theory, at various levels of simplification. An abridged QDIA is obtained by replacing the mean field trajectory, from initial-time to current-time, in the time history integrals of the non-Markovian closure by the current-time mean field. Three variants of Markovian Inhomogeneous Closures (MICs) are formulated from the abridged QDIA by using the current-time, prior-time, and correlation fluctuation dissipation theorems. The abridged MICs have auxiliary prognostic equations for relaxation functions that approximate the information in the time history integrals of the QDIA. The abridged MICs are more efficient than the QDIA for long integrations with just two relaxation functions required. The efficacy of the closures is studied in 10-day simulations with an easterly large-scale flow impinging on a conical mountain to generate rapidly growing Rossby waves in a turbulent environment. The abridged closures closely agree with the statistics of large ensembles of direct numerical simulations for the mean and transients. An Eddy Damped Markovian Inhomogeneous Closure (EDMIC), with analytical relaxation functions, which generalizes the Eddy Dampened Quasi Normal Markovian (EDQNM) to inhomogeneous flows, is formulated and shown to be realizable under the same circumstances as the homogeneous EDQNM.

1. Introduction

The development of statistical dynamical closure theory for fluid dynamical equations describing chaotic motion depends on a suitable truncation of the infinite hierarchy of coupled moment or cumulant equations. For systems described by quadratic nonlinearity, the prognostic equation for a given cumulant is coupled to the next higher cumulant. Examples of particular interest here are the Navier–Stokes equations, and the quasi-geostrophic and primitive equations of atmospheric and oceanic dynamics. If closed at second order, the mean, the one-point function, is coupled to the second order cumulant, the two-point function, which in turn is coupled to the three-point function and so on. The order at which this infinite hierarchy of moments is truncated defines the dynamics in terms of interacting triads in wavenumber space. In particular, there is an important class of closure theories for which the three-point function can be represented in terms of time history integrals over a product of three (renormalized) propagators–two-point cumulants and response functions–multiplied by vertex functions [1]. As noted in the review in the introduction of Frederiksen [2], essentially, the same closure problem occurs for corresponding systems in quantum field theory, such as quantum electrodynamics (QED) and the scalar Klein Gordon equation with $g{\varphi }^3$ Lagrangian.

The definitive pioneering advance in statistical fluid dynamical closure theory was made by Kraichnan [3] with his development of the non-Markovian Eulerian Direct Interaction Approximation (DIA) for three-dimensional homogeneous isotropic turbulence (HIT) in which the propagators are closed at second order and the renormalized vertex functions are replaced by the bare vertices, the interaction coefficients, of the Navier–Stokes equations.

Kraichnan’s [3] DIA was followed by independent approaches by Herring [4,5], resulting in the closely similar self-consistent field theory (SCFT) closure and by McComb [6,7,8,9] with his local energy transfer (LET) closure. With hindsight, the SCFT and LET closures can be derived formally from the DIA, as noted by Frederiksen et al. [10] and Kiyani and McComb [11], by using a fluctuation–dissipation theorem (FDT) [12,13,14,15,16], as an approximation. These related closures use the prior-time FDT [14] (Equation (3.5)) defined by

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

(1)

for $t\ge t^{\prime }$. Here, $C_{\mathbf{k}}(t,t^{\prime })$ is the two-time spectral covariance at wavenumber $k$, $R_{\mathbf{k}}(t,t^{\prime })$ is the response function and $C_{\mathbf{k}}(t^{\prime },t^{\prime })$ the prior time single-time covariance. Both the SCFT and LET have the same single-time second order cumulant equation as the DIA but with somewhat different response function or two-time cumulant equations. The SCFT has the same response function as the DIA but determines the two-time cumulant from Equation (1). The LET uses the DIA two-time cumulant and determines the response function from Equation (1). Carnevale and Martin [17] and Carnevale and Frederiksen [18] also developed two space scale versions of the Eulerian DIA closure for homogeneous anisotropic turbulence interacting with Rossby waves and internal gravity waves.

The subsequent further development of closures for homogeneous turbulence has been reviewed by Frederiksen and O’Kane [1], McComb [7,8,9], Lesieur [19], Zhou et al. [20], Cambon et al. [21], Sagaut and Cambon [22], and Zhou [23].

Kraichnan’s [3] theoretical approach to the statistical dynamics of HIT, and even more so the diagrammatic approaches of Wyld [24] and Lee [25], are based on application of renormalized perturbation theory to classical systems. Renormalized perturbation theory, through diagrammatic and functional approaches, were employed most famously by Tomonaga, Schwinger, and Feynman to develop the remarkably successful quantum electrodynamics (QED) theory in the mid-1900s (Frederiksen [2] reviews the literature). In QED the strength of the interaction, measured by the fine structure constant of ~$\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$137$}\right.$, is quite weak while in high Reynolds number turbulence the interactions are much stronger.

Martin et al. [26] and Phythian [27] generalized the Schwinger–Dyson functional operator approach for QED to classical systems while the Feynman path integral approach [2,28,29] was employed by Phythian [30] and Jensen [31] to extend the functional operator approach to more complex systems.

The performance of the DIA, as well the SCFT and LET closures, at the large energy containing scales, is quite accurate but the high Reynolds number power law behavior of the DIA and SCFT differ slightly from the $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ inertial range for three-dimensional turbulence and the $k^{-3}$ enstrophy cascading inertial range for two-dimensional turbulence [32]. Kraichnan [33] attributed these deficiencies to “sweeping effects”, spurious non-local interactions between the small and large scales of the turbulent eddies. McComb [6,7] showed that, at statistical steady state, the cause of the inconsistency of the DIA, SCFT, and Edward’s [34] steady state closure, with the Kolmogorov [35] classical $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ inertial range for high Reynolds number three-dimensional turbulence was associated with an infrared (zero wavenumber) divergence for the response function equation. Remarkably, the singularity disappears with the LET response function showing that the Eulerian LET is consistent with the classical $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ inertial range for high Reynolds number three-dimensional turbulence. However, at moderate Reynolds numbers or for any finite resolution, and notably for two-dimensional turbulence [36,37], all two-point two-time Eulerian closures such as the DIA, SCFT, and LET closures all have similar deficiencies with too little small-scale energy.

Kraichnan [38] noted that the correct inertial ranges could be obtained by localizing the interactions between wavenumbers of triads through a cut-off ratio $\alpha$. In fact, the studies of Frederiksen and Davies [37] and O’Kane and Frederiksen [39] suggest that $\alpha \approx 4$ is essentially universal for both two-dimensional homogeneous and inhomogeneous turbulence. This is also very close to the estimate of $\alpha \approx 3.5$ of Sudan and Pfirsch [40] for three-dimensional HIT.

Subsequently, Kraichnan [41,42], Kraichnan and Herring [43], Kaneda [44], and Gotoh et al. [45] used transformations of the fluid dynamical equations into quasi-Lagrangian coordinates to develop alternatives to the Eulerian closures. Unfortunately, unlike the Eulerian DIA, the equations and results of the quasi-Lagrangian closures depend on the choice of field variables used and whether the formulation is in terms of labelling time derivatives [41] or measuring time derivatives [44] as reviewed by Frederiksen and Davies [37]. The quasi-Lagrangian closures attempt to avoid the spurious convection effects by non-unique transformations and choices of variables but are still second order in perturbation theory and do not fundamentally address the problem of vertex renormalization. Martin et al. [26] asserted that “the whole problem of strong turbulence is contained in a proper treatment of vertex renormalization”. A fundamental resolution of the vertex renormalization problem would certainly be a phenomenal achievement for the theory of HIT. However, in the meantime a one parameter non-Markovian theory with $\alpha \approx 3.5 \mathrm{to} 4$, as discussed above, or a one parameter Markovian closure such as the EDQNM and EDMIC, discussed in Section 7, may suffice.

There is an equally important issue, stressed by McComb [9], and that is “Ultimately, if closures are to be useful, they must be capable of application to real-life situations”. This requires the development of computationally tractable closure theories for inhomogeneous systems. Second-order inhomogeneous closures have traditionally required the computation of the full covariance and response function matrices to close the mean field equation. This is the case for closures based on the renormalized perturbation theory approach of Kraichnan [46,47], on the functional approach of Martin et al. [26], on the path integral approaches [30,31], and on the corresponding Schwinger–Dyson and Schwinger–Keldysh methods for quantum field theory [2].

The computational costs of computing the full covariance and response function matrices at every timestep are currently too large, for high-dimensional systems, and are prohibitive for the non-Markovian closures with potentially long time-history integrals. Thus, to date, closures that require full covariance matrices including time-history information, i.e., correlations in time, such as Kraichnan’s [46,47] inhomogeneous DIA (IDIA) have not been applied to studies of real-life inhomogeneous turbulence problems. Some progress has been made in tackling these two important issues of (1) the size of the covariance and response function matrices and (2) the potentially long time-history integrals.

The quasi-diagonal direct interaction approximation (QDIA) closure, formulated by Frederiksen [48] for inhomogeneous turbulence interacting with mean flows and topography, represents the full covariance and response function matrices, and the three-point function, in terms of the diagonal elements and the mean flow and topography in spectral space. It thus reduces the computational problem of second order inhomogeneous closure from the computation of $N\times N$ covariance and response function elements at each time step to $N$ where $N$ is the total number of components of the dynamical fields. In addition, the coupled mean field and covariance and response function equations are much simpler than Kraichnan’s [46,47] IDIA closure and the related inhomogeneous closure formulations discussed above [2]. The QDIA closure has subsequently been generalized to include non-Gaussian initial conditions [39], to include Rossby waves on a $\beta$-plane [49], to general classical field theories with first order time derivatives [50,51] and to classical and quantum field theories with first or second order time derivatives and non-Gaussian noise and non-Gaussian initial conditions [2].

The QDIA has been numerically implemented and extensively tested against large ensembles of direct numerical simulations (DNS) and shown to be only slightly slower than the DIA for homogeneous turbulence. It has been implemented, in the bare vertex approximation and, at higher resolution and Reynolds numbers, in regularized form with $\alpha \approx 4$; it has been used to study the statistical dynamics of the interaction of inhomogeneous turbulent flows with topography on $f$-planes [39] and $\beta$-planes [49]. It has been tested in predictability studies for atmospheric blocking transitions [49,52], in data assimilation studies with square root and Kalman filters [53,54], and for developing subgrid scale parameterizations [55,56]. The structure of the QDIA closure equations has also been used as a framework for developing subgrid scale parameterizations for atmospheric and oceanic turbulent flows [57].

Progress has also been made on second major issue noted above, viz., reducing the cost of computing the potentially long time-history integrals of the non-Markovian closures that scale like $O(T^3)$ where $T$ is the total integration time. The cumulant update restart procedure [10,58] for homogeneous turbulence that uses non-Gaussian terms in the periodic restarts has been generalized for the QDIA closure to also improve its computational efficiency [39,49]. Recently, Frederiksen and O’Kane [1] formulated and implemented Markovian inhomogeneous closures (MICs) and showed that they performed very well compared with large ensemble of DNS and with the non-Markovian QDIA.

Orszag [59] formulated the eddy damped quasi-normal Markovian (EDQNM) closure for HIT for which the computational cost scales with integration time $T$ like $O(T)$. It is a one-parameter realizable closure that avoids the potential negative energy problems [60,61] of quasi-normal closures [62]. It can also be formally reduced from the Eulerian DIA by replacing the response function equation by an analytical form incorporating the eddy damping, with an empirical damping timescale parameter, and replacing the two-time cumulant equation by that determined from the current-time FDT

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

(2)

for $t\ge t^{\prime }$. At about the same time Kraichnan [63] introduced a slightly more complex Markovian closure for homogeneous turbulence, the test-field model (TFM), which, like the EDQNM, is consistent with the inertial range power laws for two-dimensional and three-dimensional turbulence. These Markovian closures have been extensively employed in the study of both isotropic and homogeneous anisotropic turbulence including in the presence of waves [17,18,21,22,64,65,66,67,68,69,70,71,72].

Bowman et al. [68] made an important point about the EDQNM for the interaction of homogeneous turbulence with waves demonstrating that the closure may not be realizable if a time-dependent analytical response function is employed rather than the oft-used asymptotic form. They noted that this is the case irrespective of whether a Markovian form is derived using the current-time FDT in Equation (2) or the prior-time FDT in Equation (1). Instead, they established a realizable Markovian closure (RMC) using an FDT that contains both current- and prior-time cumulants and that we call the correlation FDT

$$C_{\mathbf{k}}(t,t^{\prime })\equiv {[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 }$. It is clear from Equation (3) that the response function and correlation function are equal. The RMC statistical equations consist of a Markovian equation for the single-time cumulant, in which appears a triad relaxation function, coupled to a Markovian equation for the triad relaxation function. Bowman et al. [68] also formulated the theory for multi-field versions of the homogeneous RMC and Hu et al. [73] applied it to the two equation Hasegawa–Wakatani model of plasma physics. A related realizable test-field model (RTFM) closure was developed by Bowman and Krommes [74] who employed it and the RMC to studies of the interaction of homogeneous turbulence and plasma drift waves (essentially Rossby waves) in the Charney–Hasagawa–Mima equation (essentially the barotropic vorticity equation with long wave stabilization).

Frederiksen and O’Kane [1] developed and examined the performance of three new versions of Markovian Inhomogeneous Closures (MICs) using the three FDTs in Equations (1)–(3). These FDTs were combined to the form

$$C_{\mathbf{k}}(t,t^{\prime })\equiv {[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 corresponds to $X=0$, correlation FDT to $X=1/2$ and the prior-time FDT to $X=1$. In principle, realizability is only guaranteed for the variant employing the correlation FDT in Equation (3). However, for the numerical experiments carried out, the performance of all three variants was very similar with remarkable agreement with results of the non-Markovian QDIA and with large ensembles of DNS.

The broad aim of this article is to make further advances in the development of efficient closures for inhomogeneous turbulent flows. We focus on an issue that complicates the formulation of MICs derived from the non-Markovian QDIA: the fact that the closure equations for the QDIA require the complete trajectory of the mean field in the time history integrals of both the single-time cumulant equation and mean field equation. This dependence on the mean field trajectory is carried through to one of the three auxiliary relaxation functions in each of the three MIC variants of Frederiksen and O’Kane [1]. As we show in this study, if the decay of the two-time cumulants and response functions in the time history integrals can formally be assumed to be faster than the change in the mean field then the MIC equations simplify and become more efficient. The abridged MIC equations involve just two relaxation functions for each closure. Importantly, the new relaxation functions do not involve the trajectory of the mean field and, equally importantly, for the MIC employing the current-time FDT in Equation (2), the relaxation functions only involve the response functions as is the case for the EDQNM for homogeneous turbulence. In the case of HIT, the EDQNM triad relaxation function can be calculated analytically rather than through time integration and thus the EDQNM is computationally even more efficient. If the problem of possible non-realizability of the EDQNM for homogeneous anisotropic turbulence (HAT) with response functions involving the bare wave frequency [1,68] could be overcome by a suitable renormalized form then, again, the MIC with analytical relaxation functions would be vastly more computationally efficient.

In this study we focus on the case of two-dimensional inhomogeneous turbulent flows interacting with Rossby waves and topography on a generalized $\beta$-plane. Our specific goals are:

• Formulate the abridged QDIA closure equations in the case of formally slowly varying mean field components in the time history integrals;
• Examine the realizability of the abridged QDIA variant and its consistency with canonical equilibrium;
• Formulate more efficient abridged MICs for each of the FDT in Equations (1)–(3) based on the abridged QDIA;
• Evaluate the performance of the abridged QDIA and the three abridged MICs compared with large ensembles of DNS;
• Formulate an Eddy Damped Markovian Inhomogeneous Closure that has analytical representations of the relaxation functions similar to the EDQNM for homogeneous turbulence.

Although our abridged more efficient formulations assume that the mean field in the time history integrals varies more slowly than the decay of the two-time cumulant and response functions, we test these new non-Markovian and Markovian closures in situations where the mean field is spun up rapidly from very small amplitude. These simulations thus are severe tests of the closures as energy is drained from the large-scale flow and the turbulent eddy field to generate Rossby wave trains via topographic interactions.

The article is organized as follows. The equations for two-dimensional flows, consisting of a large-scale eastward wind and smaller scale circulations, interacting with topography on a generalized $\beta$-plane, are presented in Section 2. The large-scale wind evolves according to the form-drag equation and the smaller scales satisfy the barotropic vorticity equation on the doubly periodic domain. In Section 3 the flow fields are represented by Fourier series and equations for the spectral coefficients are detailed. In Section 4, the statistical dynamical equations for the non-Markovian QDIA closure are documented and the abridged QDIA formulated in which the current-time mean field replaces its complete trajectory from initial-time to current-time in the time history integrals. Three variants of abridged MIC models are derived in Section 5 from the abridged QDIA by employing the current-time, prior-time, and correlation FDTs for the two-time cumulants as well as a Markovian form of the response function. The auxiliary Markovian prognostic equations for the relaxation functions needed for the MICs are formulated in Section 5. In Section 6, numerical simulations with the abridged closures are described and the results compared with very large ensembles (with 1800 members) of DNS and with results from the original QDIA and MICs [1]. A new Markovian closure, the Eddy Damped Markovian Inhomogeneous Closure (EDMIC), is formulated in Section 7 from the abridged MIC using the current-time FDT. It has analytical forms for the relaxation functions, rather than prognostic equations, like the EDQNM for homogeneous turbulence, and is a generalization of the EDQNM to inhomogeneous turbulent flows with statistical equations for both the mean flow and the second order cumulant. In Section 8 we discuss the import of our findings and conclusions and ideas for a sequel to this work. The interaction coefficients for the spectral equations are presented in Appendix A. Relationships between the diagonal and off-diagonal spectral elements of the two-point and three-point cumulants and response functions, which are needed for deriving the QDIA closure, are listed in Appendix B. Appendix C presents the Langevin equation for the abridged QDIA for slowly varying mean field and Appendix D presents the Langevin equation for the EDMIC model.

2. Barotropic Flow over Topography on a β-Plane

The basic dynamical equations from which our closures are developed are those describing barotropic two-dimensional flows over topography. We do our analysis for flows in planar geometry on a generalized $\beta$-plane with the flows consisting of a large-scale eastward wind $U$ together with smaller scales described by the streamfunction $\psi$ and the total flow by $\Psi =\psi -Uy$.

2.1. Large-Scale Flow Equation

The equation for the large-scale flow $U$ is the so-called form-drag equation where the scaled topography $h(x)$ interacts with the smaller scales to change the large-scale flow. We include possible relaxation to a mean large-scale flow $\overline{U}$ with strength ${\alpha }_U$ so that

$$\frac{\partial U}{\partial t}=\frac{1}{{(2\pi )}^2}{\displaystyle {\int }_0^{2\pi }{d}^{2}x} h(x)\frac{\partial \psi (x)}{\partial x}+{\alpha }_U(\overline{U}-U)$$

(5)

determines the evolution of $U$. 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.2. Barotropic Vorticity Equation for the Small Scales

The barotropic vorticity equation for the smaller scales interacting with large-scale flow of Equation (5) and topography is given by

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

(6a)

Here, the Jacobian is

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

(6b)

and the vorticity $\zeta$ is the Laplacian of the streamfunction

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

(6c)

The above equations, on the generalized $\beta$-plane, are structurally the same as the barotropic equations on a sphere with the wavenumber $k_0$ the analogue of that determining the strength of the vorticity of the solid-body rotation term [49]. The viscosity is represented by $\widehat{\nu }$, the strength of the beta effect denoted by $\beta$ and $f^0$ allows for possible external forcing although our simulations are for $f^0=0$.

In the absence of forcing, topography, and viscosity, Rossby waves, and superpositions of Rossby waves, proportional to $\exp \mathrm{i}(k.x-{\omega }_k^{U}t)$ are solutions to Equation (6a–c). Here, the Doppler shifted frequency must satisfy the dispersion relationship

$${\omega }_k^U(t)={\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}.$$

(7a)

where the Rossby wave frequency is

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

(7b)

and

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

(7c)

3. Dynamical Equations in Fourier Space

The dynamical equations for Fourier components of the flow fields are established by Fourier transforms. For the ‘small-scales’ the spectral representations are determined by

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

(8a)

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

(8b)

with $x=(x,y), k=(k_x,k_y), k={(k_{\mathrm{x}}^2+k_{\mathrm{y}}^2)}^{1/2}$. The spectral coefficients satisfy the complex conjugation property ${\zeta }_{-k}={\zeta }_k^{*}$ that guarantees the reality of the physical space field. In Equation (8a), $R$ is a circular domain in wavenumber space excluding the origin $0$. In fact, through a judicious definition of the interaction coefficients we can combine the spectral equations for the ‘small-scales’ with that for $U$, taken as the zero-wavenumber component. This is achieved [49] by defining the associated vorticity components by

$${\zeta }_{-0}=\mathrm{i}{k}_{0}U, {\zeta }_0={\zeta }_{-0}^{*}.$$

(9)

The complete set of interaction coefficients $A(k,p,q)$ and $K(k,p,q)$ needed to combine the equations are documented in Appendix A. Finally, the combined equations take the form

$$\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right){\zeta }_{\mathbf{k}}(t)={\displaystyle \sum _{\mathbf{p}\in \mathbf{T}}{\displaystyle \sum _{\mathbf{q}\in \mathbf{T}}\delta (k,p,q)}}\left[K(k,p,q){\zeta }_{-\mathbf{p}}(t){\zeta }_{-\mathbf{q}}(t)+A(k,p,q){\zeta }_{-\mathbf{p}}(t)h_{-\mathbf{q}}\right]+f_{\mathbf{k}}^0(t)$$

(10)

where $T=R\cup 0$ [49]. In Equation (10), $\delta (k,p,q)=1$ if $k+p+q=0$ and 0 if $k+p+q\ne 0$, and the complex ${\nu }_0(k)k^2$ combines the viscosity and the Rossby wave frequency ${\omega }_{\mathbf{k}}$

$${\nu }_0(k)k^2=\widehat{\nu }{k}^2+\mathrm{i}{\omega }_{\mathbf{k}}$$

(11)

with the Rossby wave frequency given by Equation (7b). The drag and relaxation force in Equation (5) have also been included as the $k=0$ components of $f_{\mathbf{k}}^0$ and ${\nu }_0(k)$ through

$$f_0^0={\alpha }_U{\overline{\zeta }}_0,$$

(12)

$${\nu }_0(0)k_0^2={\alpha }_U.$$

(13)

4. Statistical Dynamical Equations for Inhomogeneous QDIA Closure

The theoretical formulation of closures is based on the statistical dynamics of an infinite ensemble of individual flows, which in our case satisfy the DNS in Equations (5) and (6a–c) in physical space and Equation (10) in Fourier space. Thus, closures have the important advantage of propagating the population statistics rather than the sampling statistics of a finite number of DNS ensemble members. Large errors can occur in estimating means, when the covariances are sizeable, unless considerable care is taken in setting up even very large ensembles of DNS [39,49].

In any given member of the infinite ensemble the vorticity spectral field at wavenumber $k$ can be written in terms of its mean $<{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t)$ and deviation from the ensemble mean ${\tilde{\zeta }}_{\mathbf{k}}(t)$:

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

(14)

From Equation (10) the coupled equations for $<{\zeta }_{\mathbf{k}}>$ and ${\tilde{\zeta }}_{\mathbf{k}}$ then follow as:

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)<{\zeta }_{\mathbf{k}}(t)>={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (k,p,q)[K(k,p,q)\left\{<{\zeta }_{-\mathbf{p}}(t)><{\zeta }_{-\mathbf{q}}(t)>+C_{-\mathbf{p},-\mathbf{q}}^{\oplus }(t,t)\right\}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+A(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{p}}(t)>h_{-\mathbf{q}}]+{\overline{f}}_{\mathbf{k}}^0(t),\end{array}$$

(15a)

and

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right){\tilde{\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)>{\tilde{\zeta }}_{-\mathbf{q}}(t)+{\tilde{\zeta }}_{-\mathbf{p}}(t)<{\zeta }_{-\mathbf{q}}(t)>}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\tilde{\zeta }}_{-\mathbf{p}}(t){\tilde{\zeta }}_{-\mathbf{q}}(t)-C_{-\mathbf{p},-\mathbf{q}}^{\oplus }(t,t)\}+A(\mathbf{k},\mathbf{p},\mathbf{q}){\tilde{\zeta }}_{-\mathbf{p}}(t)h_{-q}]+{\tilde{f}}_{\mathbf{k}}^0(t).\end{array}$$

(15b)

In these and subsequent equations $k,p,q$ are in the domain $T=R\cup 0$. The bare forcing function consists of deterministic and random contributions

$$f_{\mathbf{k}}^0(t)={\overline{f}}_{\mathbf{k}}^0(t)+{\tilde{f}}_{\mathbf{k}}^0(t)\equiv <f_{\mathbf{k}}^0(t)>+{\tilde{f}}_{\mathbf{k}}^0(t),$$

(16a)

and the two-time cumulant elements are given by

$$C_{-p,-q}^{\oplus }(t,s)\equiv C_{-\mathbf{p},-\mathbf{q}}(t,s)= <{\tilde{\zeta }}_{-\mathbf{p}}(t){\tilde{\zeta }}_{-\mathbf{q}}(s)>.$$

(16b)

Here, we have placed the superscript $\oplus$ on $C_{-p,-q}^{\oplus }(t,s)$ for later convenience where it will be replaced by QDIA and the QDIA expression inserted to close the equations.

4.1. General Formulation of QDIA Closure

We see from Equation (15a) that to close the mean field equation we need an expression for the off-diagonal elements of the two-point cumulant. The QDIA theory for closing the mean field and second order cumulant equations were originally derived by Frederiksen [48], extended to the current spectral Equation (10) by Frederiksen and O’Kane [49], and to general systems of classical and quantum field theory equations by Frederiksen [2,50,51]. The QDIA closure equation for $<{\zeta }_{\mathbf{k}}(t)>$ is obtained by replacing $C_{-p,-q}^{\oplus }(t,s)\equiv C_{-\mathbf{p},-\mathbf{q}}(t,s)$ in Equation (15a) by $C_{-p,-q}^{QDIA}(t,s)$ defined in Equation (A4) of Appendix B. The expression for $C_{-p,-q}^{QDIA}(t,s)\equiv C_{-p,-q}^{QDIA}(t,s)[C_k,R_k,<{\zeta }_k>,h_k]$ is a functional of the diagonal two-point cumulant $C_{\mathbf{k}}\equiv C_{k,-k}$, the diagonal response function $R_{\mathbf{k}}\equiv R_{k,k}$, mean field $<{\zeta }_{\mathbf{k}}>$ and topography $h_{\mathbf{k}}$. Thus, we need equations for $C_{\mathbf{k}}$ and $R_{\mathbf{k}}$.

The general matrix element of the response function, defined for $t\ge t^{\prime }$, is a measure of the change in an ensemble member field ${\tilde{\zeta }}_{\mathbf{k}}(t)$ at time $t$ due to an infinitesimal change in the forcing ${\tilde{f}}_l^0(t^{\prime })$ at the earlier time $t^{\prime }$. It is defined by through the functional derivative

$${\tilde{R}}_{\mathbf{k},l}(t,t^{\prime })=\frac{\delta {\tilde{\zeta }}_k(t)}{\delta {\tilde{f}}_l^0(t^{\prime })},$$

(17a)

where $\delta$ denotes the functional derivative. We shall need the ensemble average response function

$$R_{\mathbf{k},l}(t,t^{\prime })=\langle {\tilde{R}}_{\mathbf{k},l}(t,t^{\prime })\rangle .$$

(17b)

As noted above we use a shortened notation for the diagonal elements of the response function and two-point cumulant:

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

(17c)

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

(17d)

The statistical dynamical equation for the diagonal two-time cumulant $C_{\mathbf{k}}$ can be obtained from Equation (15b), by multiplying each term by ${\tilde{\zeta }}_{-\mathbf{k}}(t^{\prime })$ and averaging over the infinite ensemble. This yields the statistical equation

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)C_{\mathbf{k}}(t,t^{\prime })={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})A(\mathbf{k},\mathbf{p},\mathbf{q})C_{-p,-k}^{\oplus }(t,t^{\prime })h_{-\mathbf{q}}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\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)>C_{-q,-k}^{\oplus }(t,t^{\prime })}}+C_{-p,-k}^{\oplus }(t,t^{\prime })<{\zeta }_{-\mathbf{q}}(t)>\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+<{\tilde{\zeta }}_{-\mathbf{p}}(t) {\tilde{\zeta }}_{-\mathbf{q}}(t) {\tilde{\zeta }}_{-\mathbf{k}}(t^{\prime }){>}^{\oplus }]+<{\tilde{f}}_k^0(t){\tilde{\zeta }}_{-\mathbf{k}}(t^{\prime })>\end{array}$$

(18)

where $t>t^{\prime }$ and $C_k(t,t^{\prime })=C_{-\mathbf{k}}(t^{\prime },t)$ for $t^{\prime }>t$. Again, the superscript $\oplus$ is used to denote the terms that need to be replaced by their QDIA expressions in Appendix B to produce the QDIA closure equations. Additionally, note that

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

(19a)

where $t_0$ is the initial time and

$$F_{\mathbf{k}}^0(t,s)= <{\tilde{f}}_{\mathbf{k}}^0(t){\tilde{f}}_{\mathbf{k}}^{0\ast }(s)>.$$

(19b)

The statistical dynamical equation for the diagonal response function $R_{\mathbf{k}}$ is obtained from Equation (15b) by taking the functional derivative with respect to ${\tilde{f}}_k^0(t^{\prime })$ and ensemble averaging to give

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)R_{\mathbf{k}}(t,t^{\prime })={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})A(\mathbf{k},\mathbf{p},\mathbf{q})R_{-p,k}^{\oplus }(t,t^{\prime })h_{-\mathbf{q}}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\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)>R_{-q,k}^{\oplus }(t,t^{\prime })}}+R_{-p,k}^{\oplus }(t,t^{\prime })<{\zeta }_{-\mathbf{q}}(t)>\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+2<{\tilde{R}}_{-\mathbf{p},k}(t,t^{\prime }) {\tilde{\zeta }}_{-\mathbf{q}}(t) {>}^{\oplus }]+\delta (t-t^{\prime })\end{array}$$

(20)

where $t\ge t^{\prime }$ and the Dirac delta function means that $R_k(t,t)=1$. The terms with superscript $\oplus$ are again to be replaced by the QDIA forms in Appendix B.

The four QDIA expressions in Appendix B, that replace the terms with superscript $\oplus$ in Equations (15a), (18) and (20), at most depend functionally on the diagonal two-point cumulant $C_{\mathbf{k}}$, the diagonal response function $R_{\mathbf{k}}$, mean field $<{\zeta }_{\mathbf{k}}>$, and topography $h_{\mathbf{k}}$. Thus, these three equations together with Equations (A4)–(A7) form the QDIA closure. The mean field equation, the diagonal two-point cumulant equation and the diagonal response function can then be written out in full.

The mean field equation is

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)<{\zeta }_{\mathbf{k}}(t)> +{\displaystyle \underset{t_o}{\overset{t}{\int }}ds {\eta }_{\mathbf{k}}(t,s)<{\zeta }_{\mathbf{k}}(s)>}={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})A(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{p}}(t)>h_{-\mathbf{q}}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\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_{\mathbf{k}}^h(t)+{\overline{f}}_{\mathbf{k}}^0(t).\end{array}$$

(21)

The diagonal two-time cumulant equation is

$$\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)C_{\mathbf{k}}(t,t^{\prime })+{\displaystyle \underset{t_0}{\overset{t}{\int }}ds ({\eta }_k(t,s)+{\pi }_k(t,s))C_{-\mathbf{k}}(t^{\prime },s)}={\displaystyle \underset{t_0}{\overset{t^{\prime }}{\int }}ds (S_{\mathbf{k}}(t,s)+P_{\mathbf{k}}(t,s)+F_{\mathbf{k}}^0(t,s))R_{-\mathbf{k}}(t^{\prime },s)}$$

(22)

where $t>t^{\prime }$ with $C_k(t,t^{\prime })=C_{-\mathbf{k}}(t^{\prime },t)$ for $t^{\prime }>t$. The diagonal response function equation is

$$\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)R_{\mathbf{k}}(t,t^{\prime })+{\displaystyle \underset{t^{\prime }}{\overset{t}{\int }}ds({\eta }_{\mathbf{k}}(t,s)+{\pi }_{\mathbf{k}}(t,s))R_{\mathbf{k}}(s,t^{\prime })}=\delta (t-t^{\prime })$$

(23)

for $t\ge t^{\prime }$ and $R_k(t,t^{\prime })=1$.

In Equations (21)–(23)

$$f_{\mathbf{k}}^h(t)\equiv f_{\mathbf{k}}^{\chi }(t)=h_{\mathbf{k}}{\displaystyle \underset{t_o}{\overset{t}{\int }}ds {\chi }_{\mathbf{k}}(t,s),}$$

(24a)

$${\chi }_{\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})A(-\mathbf{p},-\mathbf{q},-\mathbf{k})R_{-\mathbf{p}}(t,s)C_{-\mathbf{q}}(t,s)}},$$

(24b)

$${\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),}}$$

(24c)

$$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)}},$$

(24d)

$$\begin{array}{l}{\pi }_{\mathbf{k}}(t,s)=-{\displaystyle \sum _p{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})R_{-\mathbf{p}}(t,s)[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{p},-\mathbf{k},\mathbf{-}\mathbf{q})<{\zeta }_{\mathbf{q}}(s)>+A(-\mathbf{p},-\mathbf{k},-\mathbf{q})h_{\mathbf{q}}],\end{array}$$

(24e)

$$\begin{array}{l}{P}_{\mathbf{k}}(t,s)={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})C_{-\mathbf{p}}(t,s)[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{k},-\mathbf{p},-\mathbf{q})<{\zeta }_{\mathbf{q}}(s)>+A(-\mathbf{k},-\mathbf{p},-\mathbf{q})h_{\mathbf{q}}].\end{array}$$

(24f)

These terms modify or renormalize the damping or forcing in the mean field, two-point cumulant and response function equations. The eddy-topographic force $f_{\mathbf{k}}^h(t)$ is a product of the topography and the time history integral of ${\chi }_{\mathbf{k}}(t,s)$, which measures of the strength of the interaction of the transient eddies with the topography. The term ${\eta }_{\mathbf{k}}(t,s)$ is the nonlinear damping due to eddy–eddy interactions that also damps the two-point cumulant and response function. The cumulant and response functions are also damped by ${\pi }_{\mathbf{k}}(t,s)$ due to interaction of transient eddies with the mean flow and topography. The two-point cumulant equation also contains the nonlinear noise terms, $S_{\mathbf{k}}(t,s)$ and $P_{\mathbf{k}}(t,s)$ that renormalize the specified bare noise spectrum $F_{\mathbf{k}}^0(t,s)$. The term $S_{\mathbf{k}}(t,s)$ is due to eddy–eddy interactions and $P_{\mathbf{k}}(t,s)$ is arises from eddy–mean flow and eddy–topographic interactions. All these noise terms are positive semi-definite and represent energy and enstrophy injection.

The last statistical dynamical equation we need 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\right)C_{\mathbf{k}}(t,t)+2 \mathrm{Re} {\displaystyle \underset{t_0}{\overset{t}{\int }}ds({\eta }_{\mathbf{k}}(t,s)+{\pi }_{\mathbf{k}}(t,s))C_{-\mathbf{k}}(t,s)}\\ =2 \mathrm{Re} {\displaystyle \underset{t_0}{\overset{t}{\int }}ds(S_{\mathbf{k}}(t,s)+P_{\mathbf{k}}(t,s)+F_{\mathbf{k}}^0(t,s))R_{-\mathbf{k}}(t,s)}\end{array}$$

(25)

where the initial conditions $C_k(t_0,t_0)$ are to be specified.

4.2. Abridged QDIA Closure with Current-Time Mean Field

Next, we consider the case when the mean field $<{\zeta }_{\mathbf{k}}(s)>$ is replaced by the current-time mean field $<{\zeta }_{\mathbf{k}}(t)>$ in the time history integrals:

$${\overline{\zeta }}_{\mathbf{k}}(s)\equiv <{\zeta }_{\mathbf{k}}(s)>\to <{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t).$$

(26)

This corresponds to the mean field being more slowly varying than the response function and the two-time cumulant in the time history integrals. This, like the different forms of the fluctuation theorem considered in Frederiksen and O’Kane [1] and below, allows further simplification and greater efficiency of the Markovian inhomogeneous closures that we studied there. For ease of discussion we will denote this abridged version of the QDIA, with the replacement (26), by $QDIA[\overline{\zeta }(t)]$ in contrast to the original QDIA which we also denote by $QDIA[\overline{\zeta }(s)]$ since for it the time history integrals depend on the whole trajectory of ${\overline{\zeta }}_{\mathbf{k}}(s)$ for $t_0\le s\le t$.

Employing Equation (26) in Equation (21), the mean field equation for $QDIA[\overline{\zeta }(t)]$ becomes

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2+{\displaystyle \underset{t_o}{\overset{t}{\int }}ds {\eta }_{\mathbf{k}}(t,s)}\right)<{\zeta }_{\mathbf{k}}(t)> ={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})A(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{p}}(t)>h_{-\mathbf{q}}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\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_{\mathbf{k}}^h(t)+{\overline{f}}_{\mathbf{k}}^0(t).\end{array}$$

(27)

It is also clear from Equation (27) that ${\displaystyle \underset{t_o}{\overset{t}{\int }}ds {\eta }_{\mathbf{k}}(t,s)}$ renormalizes ${\nu }_0(k)k^2$, the bare dissipation and Rossby wave frequency, in Equation (11). For the abridged $QDIA[\overline{\zeta }(t)]$ closure the expressions for ${\pi }_{\mathbf{k}}(t,s)$ and $P_{\mathbf{k}}(t,s)$ in Equation (24a–f) also simplify, with the replacement in Equation (26), to

$$\begin{array}{l}{P}_{\mathbf{k}}(t,s)\doteq P_k^{\overline{\zeta }(t)}(t,s)={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{k},-\mathbf{p},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{k},-\mathbf{p},-\mathbf{q})h_{\mathbf{q}}]C_{-\mathbf{p}}(t,s)\end{array}$$

(28a)

$$\begin{array}{l}{\pi }_{\mathbf{k}}(t,s)\doteq {\pi }_k^{\overline{\zeta }(t)}(t,s)=-{\displaystyle \sum _p{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{p},-\mathbf{k},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{p},-\mathbf{k},-\mathbf{q})h_{\mathbf{q}}]R_{-\mathbf{p}}(t,s),\end{array}$$

(28b)

with ${\overline{\zeta }}_{\mathbf{k}}(t)\equiv <{\zeta }_{\mathbf{k}}(t)>.$ Thus, the two-time and single-time cumulant equations are again given by Equations (22) and (25), respectively, with ${\pi }_{\mathbf{k}}(t,s)\to {\pi }_k^{\overline{\zeta }(t)}(t,s)$ and $P_{\mathbf{k}}(t,s)\to P_k^{\overline{\zeta }(t)}(t,s)$ and the response function equation is given by Equation (23) with ${\pi }_{\mathbf{k}}(t,s)\to {\pi }_k^{\overline{\zeta }(t)}(t,s)$. Again, in the abridged $QDIA[\overline{\zeta }(t)]$ closure the mean field (and topography) terms can now be taken outside the time history integrals in Equations (22), (23) and (25). However, for the sake of brevity, we shall not show this explicitly here, but this fact results in the simplifications of the new abridged Markovian inhomogeneous equations compared with those of Frederiksen and O’Kane [1].

5. Statistical Dynamical Equations for Abridged Markovian Inhomogeneous Closure

The EDQNM for homogeneous turbulence is a Markovian equation for the single time cumulant that incorporates a triad relaxation time with an empirical eddy damping. It does not evolve the mean field. In this Section, where we formulate Markovian inhomogeneous closures (MICs) we start with the single-time cumulant equation and then follow with the mean field equation. We consider three variants of the single-time MICs that are formulated from the non-Markovian abridged QDIA in Equation (25), denoted $QDIA[\overline{\zeta }(t)]$, with the replacements ${\pi }_{\mathbf{k}}(t,s)\to {\pi }_k^{\overline{\zeta }(t)}(t,s)$ and $P_{\mathbf{k}}(t,s)\to P_k^{\overline{\zeta }(t)}(t,s)$ in Equation (28a,b). We denote these abridged MICs by the notation $MI{C}^X[\overline{\zeta }(t)]$ to distinguish them from those of Frederiksen and O’Kane [1] that we shall refer to as $MI{C}^X[\overline{\zeta }(s)]$ since one of their relaxation functions involve the whole trajectory of ${\zeta }_{\mathbf{k}}(s)$ for $t_0\le s\le t$. Here, the superscript $X$ relates to that used in the combined FDT in Equation (4) with $X=0$ being the current-time FDT, $X=1/2$ correlation FDT, and $X=1$ the prior-time FDT.

The non-Markovian single-time covariance equation can be written in the form

$$\frac{\partial }{\partial t}{C}_{\mathbf{k}}(t,t)+2 \mathrm{Re}[N_k^{\eta }(t)+N_k^{\pi }(t)+N_k^0]=2 \mathrm{Re}[F_k^S(t)+F_k^P(t)+F_k^0(t)]$$

(29)

where $C_{\mathbf{k}}(t,t)$ is real. Here, with Equation (26) employed, we no longer need to split up the $P_k^{\overline{\zeta }(t)}(t,s)$ and ${\pi }_k^{\overline{\zeta }(t)}(t,s)$ terms into their mean vorticity, topographic and cross terms for developing MICs. This results in considerable simplification and efficiency of the subsequent Markovian closures derived from the QDIA closure. Some algebra shows that the $F_k^{}(t)$ and $N_k^{}(t)$ functions have the following expressions

$$F_k^S(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),$$

(30a)

$$\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),$$

(30b)

$$\begin{array}{l}{F}_k^P(t)={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{k},-\mathbf{p},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{k},-\mathbf{p},-\mathbf{q})h_{\mathbf{q}}]\Lambda (-\mathbf{k},-\mathbf{p})(t),\end{array}$$

(30c)

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

(30d)

$$F_k^0(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{F}_k^0(t,s)R_{-k}(t,s).$$

(30e)

Additionally,

$$N_k^{\eta }(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),$$

(31a)

$$\begin{array}{l}{N}_k^{\pi }(t)=-{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{p},-\mathbf{k},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{p},-\mathbf{k},-\mathbf{q})h_{\mathbf{q}}]\Lambda (-\mathbf{p},-\mathbf{k})(t),\end{array}$$

(31b)

$$N_k^0(t)={\nu }_0(k)k^2{C}_k(t,t).$$

(31c)

The nonlinear noise and damping terms in Equations (30a–e) and (31a–c) simplify on applying the FDTs in Equation (4) with the time history integrals expressed through the relaxation functions ${\Theta }^X$ and ${\Psi }^X$. Importantly, the relaxation functions can alternatively be determined by time dependent differential equations and effect the Markovianization. Thus,

$$F_k^S(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})C_{-p}^{1-X}(t,t)C_{-q}^X(t,t){\Theta }^X(}}-\mathbf{k},-\mathbf{p},-\mathbf{q})(t),$$

(32a)

$${\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),$$

(32b)

$$\begin{array}{l}{F}_k^P(t)={\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{k},-\mathbf{p},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{k},-\mathbf{p},-\mathbf{q})h_{\mathbf{q}}]C_{-p}^{1-X}(t,t){\Psi }^X(-\mathbf{k},-\mathbf{p})(t),\end{array}$$

(32c)

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

(32d)

$$F_k^0(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}ds}{F}_k^0(t,s)R_{-k}(t,s).$$

(32e)

Additionally,

$$N_k^{\eta }(t)=D_k^{\eta }(t)C_k(t,t),$$

(33a)

$$D_k^{\eta }(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}^{1-X}(t,t)C_{-k}^{-X}(t,t){\Theta }^X(}}-\mathbf{p},-\mathbf{q},-\mathbf{k})(t),$$

(33b)

$$N_k^{\pi }(t)=D_k^{\pi }(t)C_k(t,t),$$

(33c)

$$\begin{array}{l}{D}_k^{\pi }(t)=-{\displaystyle \sum _{\mathbf{p}}{\displaystyle \sum _{\mathbf{q}}\delta (\mathbf{k},\mathbf{p},\mathbf{q})[2K(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{q}}(t)>+A(\mathbf{k},\mathbf{p},\mathbf{q})h_{-\mathbf{q}}]}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\times [2K(-\mathbf{p},-\mathbf{k},-\mathbf{q})<{\zeta }_{\mathbf{q}}(t)>+A(-\mathbf{p},-\mathbf{k},-\mathbf{q})h_{\mathbf{q}}]C_{-k}^{-X}(t,t){\Psi }^X(-\mathbf{p},-\mathbf{k})(t),\end{array}$$

(33d)

$$N_k^0(t)=D_k^0(t)C_k(t,t),$$

(33e)

$$D_k^0={\nu }_0(k)k^2.$$

(33f)

We can simplify the single-time cumulant equation further by defining

$$\begin{array}{l}{F}_k^1(t)=F_k^S(t); F_k^2(t)=F_k^P(t); \\ D_k^1(t)=D_k^{\eta }(t); D_k^2(t)=D_k^{\pi }(t);\end{array}$$

(34)

so that

$$\left(\frac{\partial }{\partial t}+2 \mathrm{Re}{\displaystyle \sum _{j=0}^2{D}_k^j(t)}\right)C_{\mathbf{k}}(t,t)=2 \mathrm{Re}{\displaystyle \sum _{j=0}^2{F}_k^j(t)}.$$

(35)

Now, to complete the Markovianization, the response function equation (23) with ${\pi }_{\mathbf{k}}(t,s)\to {\pi }_k^{\overline{\zeta }(t)}(t,s)$ must also be replaced by

$$\frac{\partial }{\partial t}{R}_{\mathbf{k}}(t,t^{\prime })+{\displaystyle \sum _{j=0}^2{D}_k^j(t)R_{\mathbf{k}}(t,t^{\prime })}=\delta (t-t^{\prime }).$$

(36)

We note from Equation (25), with the current-time expressions in Equation (28a,b), and from Equations (29), (33a–f) and (34) that

$${\displaystyle \sum _{j=0}^2{D}_k^j(t)}={\nu }_0(k)k^2+{\displaystyle \underset{t_0}{\overset{t}{\int }}ds\{({\eta }_k(t,s)+{\pi }_k^{\overline{\zeta }(t)}(t,s))C_{-\mathbf{k}}(t,s){[C_k(t,t)]}^{-1}\}}$$

(37)

with $C_{-\mathbf{k}}(t,s)$ determined by the FDT in Equation (4). Thus, under these same conditions, the response function in Equation (36) is equivalent to

$$\frac{\partial }{\partial t}{R}_{\mathbf{k}}(t,t^{\prime })+\left[{\nu }_0(k)k^2+{\displaystyle \underset{t_0}{\overset{t}{\int }}ds\{({\eta }_{\mathbf{k}}(t,s)+{\pi }_k^{\overline{\zeta }(t)}(t,s))C_{-\mathbf{k}}(t,s){[C_k(t,t)]}^{-1}\}}\right]R_k(t,t^{\prime })=\delta (t-t^{\prime }).$$

(38)

As foreshadowed above, with the response function equation having the simpler form in Equation (36), the integral forms for the relaxation functions ${\Theta }^X$ and ${\Psi }^X$ can be replaced by differential equations. Thus, the integral expression for the relaxation time

$${\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),$$

(39a)

can equally be calculated by forward integration of the ordinary differential equation

$$\frac{\partial }{\partial t} {\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(t)+{\displaystyle \sum _{j=0}^2[D_k^j(t)+D_p^j(t)+D_q^j(t)]{\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(t)}=C_p^X(t,t)C_q^X(t,t)$$

(39b)

with ${\Theta }^X(\mathbf{k},\mathbf{p},\mathbf{q})(0)=0$. As well

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

(39c)

can be calculated from

$$\frac{\partial }{\partial t} {\Psi }^X(\mathbf{k},\mathbf{p})(t)+{\displaystyle \sum _{j=0}^2[D_k^j(t)+D_p^j(t)]{\Psi }^X(\mathbf{k},\mathbf{p})(t)}=C_p^X(t,t)$$

(39d)

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

The statistical dynamics of the single-time cumulant is now determined by the Markovian form in Equation (35) with auxiliary Markovian equations for the relaxation times in Equation (39b,d).

Next, we aim to formulate manifestly Markovian equations for the mean field and thereby have a system of coupled Markovian equations for the mean and two-point cumulant. Again, we approximate the mean field $<{\zeta }_{\mathbf{k}}(s)>$ by the current-time mean field $<{\zeta }_{\mathbf{k}}(t)>$ in the time history integrals as in Equation (26). Then, from Equation (21) we have

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+{\nu }_0(k)k^2\right)<{\zeta }_{\mathbf{k}}(t)>+N_k^M(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)>}}\\ \hspace{1em}\hspace{1em}+A(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{p}}(t)>h_{-\mathbf{q}}]+f_k^h(t)+{\overline{f}}_k^0(t).\end{array}$$

(40)

Here, the nonlinear damping acting on the mean field and eddy-topographic force are given by

$$N_k^M(t)=-4<{\zeta }_k(t)>{\displaystyle \sum _p{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})}}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds{R}_{-p}(t,s)C_{-q}(t,s)},$$

(41a)

$$f_k^h(t)=2h_k{\displaystyle \sum _p{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})A(-\mathbf{p},-\mathbf{q},-\mathbf{k})}}{\displaystyle \underset{t_0}{\overset{t}{\int }}ds{R}_{-p}(t,s)C_{-q}(t,s).}$$

(41b)

Therefore, with $\Lambda$ defined in Equation (30d),

$$N_k^M(t)=-4<\zeta (t)>{\displaystyle \sum _p}{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})\Lambda (-\mathbf{p},-\mathbf{q})(t)},$$

(42a)

and

$$f_k^h(t)=2h_k{\displaystyle \sum _p}{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})A(-\mathbf{p},-\mathbf{q},-\mathbf{k})\Lambda (-\mathbf{p},-\mathbf{q})(t)}.$$

(42b)

Now, implementing the FDT in Equation (4) these expressions become

$$N_k^M(t)=D_k^M(t)<{\zeta }_k(t)>$$

(43a)

where

$$D_k^M(t)=-4{\displaystyle \sum _p}{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})K(-\mathbf{p},-\mathbf{q},-\mathbf{k})}{C}_{-q}^{1-X}(t,t){\Psi }^X(-\mathbf{p},-\mathbf{q})(t),$$

(43b)

and

$$f_k^h(t)=2h_k{\displaystyle \sum _p}{\displaystyle \sum _q\delta (\mathbf{k},\mathbf{p},\mathbf{q})K(\mathbf{k},\mathbf{p},\mathbf{q})A(-\mathbf{p},-\mathbf{q},-\mathbf{k})}{C}_{-q}^{1-X}(t,t){\Psi }^X(-\mathbf{p},-\mathbf{q})(t).$$

(43c)

The relaxation function ${\Psi }^X$ is again calculated through the auxiliary Markovian form in Equation (39d) and the mean field prognostic in Equation (40) simplifies to the abridged Markovian equation

$$\begin{array}{l}\left(\frac{\partial }{\partial t}+D_k^0+D_k^M(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)>}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+A(\mathbf{k},\mathbf{p},\mathbf{q})<{\zeta }_{-\mathbf{p}}(t)>h_{-\mathbf{q}}]+f_k^h(t)+{\overline{f}}_k^0.\end{array}$$

(44)

Equation (44), for the mean flow $<{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t)$, and Equation (35), for the single-time cumulant $C_{\mathbf{k}}(t,t)$, together with the auxiliary Equation (39b,d) for the relaxation functions ${\Theta }^X(\mathbf{k},p,q)(t)$ and ${\Psi }^X(\mathbf{k},p)(t)$, form the prognostic equations for the three abridged MICs with $X=0, \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right., 1$. For the abridged MICs there are two relaxation functions to be calculated while for the original MICs of Frederiksen and O’Kane [1] there were three. This of course reduces the computational task for the abridged MICs correspondingly. Importantly, from the abridged MIC with $X=0$ we derive the EDMIC model with analytical representation of the relaxation functions in Section 7. It is a generalization of the EDQNM to inhomogeneous turbulent flows and like the EDQNM is still more computationally efficient because the prognostic equations for the relaxation functions do not need to be solved.

6. Comparison of Non-Markovian and Markovian Closure Integrations with DNS

In this Section, we compare the performance of the non-Markovian abridged $QDIA[\overline{\zeta }(t)]$ closure, and the three abridged $MI{C}^X[\overline{\zeta }(t)]$ Markovian closures, with each other, and with an ensemble of 1800 direct numerical simulations. These results are also compared with those of Frederiksen and O’Kane [1] for the $QDIA[\overline{\zeta }(s)]$ and $MI{C}^X[\overline{\zeta }(s)]$ models. Between them, these numerical simulations provide insights into the robustness of the inhomogeneous closure calculations to Markovianization, with three versions of the FDT, and to whether the time history integrals, or equivalent relaxation functions, are significantly impacted by the full trajectory of $\overline{\zeta }(s)\equiv <\zeta (s)>$ for $t_0\le s\le t$ or not.

Our aim is to provide severe tests of these different formulations in situations where the mean flow is rapidly evolving and interacting with turbulence and topography. This is achieved using the setup of Frederiksen and O’Kane [1] where a large-scale mean eastward flow $U$ interacts with topography in a turbulent environment. Because of the differential rotation on a $\beta$-plane, Rossby waves are generated that interact with the topography to provide a form drag on the large-scale flow $U$. There is a rapid transfer of energy from the large-scale flow to the smaller scale mean field, which spins up rapidly, as well as wave-turbulence interactions and changes in the wavenumber distribution of transient energy.

The experimental setup is as follows. We use length and time scales of one half the earth’s radius, $a_e/2$, and the inverse of the earth’s rotation rate, ${\Omega }^{-1}$, respectively. At the initial time the mean eastward flow $U$ has a speed of $7.5{ \mathrm{ms}}^{-1}$ (non-dimensional $U=0.0325$), the $\beta$-effect is $1.15\times {10}^{-11}{ \mathrm{m}}^{-1}{\mathrm{s}}^{-1}$ (non-dimensional $\beta =\frac{1}{2}$) representative of the earth at 60° latitude, and the coefficient of viscosity $\widehat{\nu }$ is $2.5\times {10}^4{ \mathrm{m}}^2{\mathrm{s}}^{-1}$ (non-dimensional $\widehat{\nu }=3.378\times {10}^{-5}$). Additionally, the forcing $f^0=0$, the drag on the large-scale flow ${\alpha }_U=0$ and $k_0^2=\frac{1}{2}$. The topography is a 2.5 km high cone centered at 30° N, 180° E with and a diameter of 45° latitude [49] (Figure 1) and might be seen as an idealized representation of the Himalayas. The DNS and closure calculations are performed at a resolution of circular truncation C16 where $\left|k\right|=k\le 16$ and all integrations proceed for 10 days with a time step of $\Delta t=1/30 \mathrm{day}$ (non-dimensional $\Delta t=0.21$).

Table 1. Initial conditions used in DNS and closure calculations.

$C_k(0,0)$	$<{\zeta }_k(0)>$	$a$	$b$
$\frac{k^2\times {10}^{-2}}{a+b{k}^2}$	$-10b{h}_k{C}_k(0,0)$	$4.824\times {10}^4$	$2.511\times {10}^3$

specifies the initial transient isotropic spectrum and initial ‘small-scale’ mean field, which is localized over the topography and of small amplitude. The closure calculations start from ${\overline{\zeta }}_k(t_0)=<{\zeta }_k(t_0)>=<{\zeta }_k(0)>$ and Gaussian initial conditions with $C_k(t_0,t_0)=C_k(0,0)$ in Table 1. For the DNS, an ensemble of 1800 simulations is started from the mean field plus different Gaussian isotropic perturbations with spectrum $C_k(t_0,t_0)=C_k(0,0)$ in Table 1. Further details on the setup of the DNS perturbations are given in [49] (Section 6). The time stepping in both DNS and closures is performed with a predictor-corrector procedure and the time history integrals in the non-Markovian closures are calculated using the trapezoidal rule [10,49].

The evolved DNS and closure results are compared in terms of the similarity of the mean flow fields, their pattern correlations, and their mean and transient kinetic energy and palinstrophy spectra. The mean, transient, and total kinetic energy spectra, averaged over circular bands, are defined by

$$\overline{E}(k_i,t)=\frac{1}{2}{\displaystyle \sum _{k\in S}\langle {\zeta }_k(t)\rangle \langle {\zeta }_{-k}(t)\rangle k^{-2}},$$

(45a)

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

(45b)

$$E(k_i,t)=\overline{E}(k_i,t)+\tilde{E}(k_i,t),$$

(45c)

and the set $S$ is defined by

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

(46)

The band-average is over all $k$ within a band of unit width at a radius $k_i$; the energy of the large-scale flow is plotted at zero wavenumber. Similarly, the band-averaged mean, transient, and total palinstrophy spectra are defined by

$$\overline{P}(k_i,t)=\frac{1}{2}{\displaystyle \sum _{k\in S}{k}^2\langle {\zeta }_k(t)\rangle \langle {\zeta }_{-k}(t)\rangle },$$

(47a)

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

(47b)

$$P(k_i,t)=\overline{P}(k_i,t)+\tilde{P}(k_i,t).$$

(47c)

The initial mean non-zonal streamfunction for the simulations is shown in part (a) of Figure 1. The other panels in Figure 1 show the day 10 evolved mean streamfunctions for the non-Markovian and Markovian abridged closures and the ensemble of DNS. The substantial evolution in terms of magnitude and structure from the initial mean field is very evident in all simulations with large scale Rossby wavetrains primarily downstream of the conical mountain. There is little that can be said to distinguish between the panels apart from slight variations in the magnitude of the peaks and troughs of the wavetrains. We note that, at the first and largest peak downstream of the conical mountain, the value for the $MI{C}^0[\overline{\zeta }(t)]$ (at 0.0139) agrees best with the ensemble of DNS (at 0.0137). The values for the $QDIA[\overline{\zeta }(t)]$ (at 0.0128), $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(t)]$ (at 0.0128) and $MI{C}^1[\overline{\zeta }(t)]$ (at 0.0129) are slightly smaller. We also note from

Table 2. Correlations of day 10 non-zonal streamfunction for closures with DNS.

$QDIA[\overline{\zeta }(t)]:0.9789 MI{C}^{X=0}[\overline{\zeta }(t)]:0.9999 MI{C}^{X=\frac{1}{2}}[\overline{\zeta }(t)]:0.9994 MI{C}^{X=1}[\overline{\zeta }(t)]:0.9995$

$QDIA[\overline{\zeta }(s)]:0.9998 MI{C}^{X=0}[\overline{\zeta }(s)]:0.9999 MI{C}^{X=\frac{1}{2}}[\overline{\zeta }(s)]:0.9998 MI{C}^{X=1}[\overline{\zeta }(s)]:0.9998$

that, of the abridged closures, the pattern correlation of the DNS mean non-zonal streamfunctions is largest with the $MI{C}^0[\overline{\zeta }(t)]$ at 0.9999. We also see that the pattern correlation is the least with the $QDIA[\overline{\zeta }(t)]$ (at 0.9789) and nearly equal with the $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(t)]$ (at 0.9994) and $MI{C}^1[\overline{\zeta }(t)]$ (at 0.9995). In general, these results are only slightly less notable than those of Frederiksen and O’Kane [1] where the full trajectory of $\overline{\zeta }(s)$ was used in the time history integrals. The corresponding pattern correlations of the DNS field with $MI{C}^0[\overline{\zeta }(s)]$ is identical to that with $MI{C}^0[\overline{\zeta }(t)]$ at 0.9999, and with all of $QDIA[\overline{\zeta }(s)]$, $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(s)]$ and $MI{C}^1[\overline{\zeta }(s)]$ it is 0.9998. It is interesting that the replacement $\overline{\zeta }(s)\to \overline{\zeta }(t)$ has most effect on the non-Markovian QDIA closure that does not employ any of the FDTs and least effect on the $MI{C}^0$ which uses the current-time FDT. Of course, in all cases the pattern correlations are remarkably high particularly given the dramatic evolution the flow field has undergone.

Next, we consider the energy and palinstrophy spectra of the non-Markovian $QDIA[\overline{\zeta }(t)]$ and each of the Markovian $MI{C}^X[\overline{\zeta }(t)]$ models with $X=0,\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,1$. Figure 2 and Figure 3 show the initial and 10-day evolved mean and transient kinetic energy spectra and palinstrophy spectra, respectively, for these closures and for comparison also the results for the ensemble of DNS. The palinstrophy spectra amplify any differences at small scales as expected. The $QDIA[\overline{\zeta }(t)]$, $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(t)]$ and $MI{C}^1[\overline{\zeta }(t)]$ have slightly larger evolved transient energy near the peak at $k=4$ while for $MI{C}^0[\overline{\zeta }(t)]$ it is underestimated. However, all the closures perform remarkably well compared with the ensemble of DNS. Using the current-time mean field in the time history integrals does not significantly degrade the abridged closure simulations compared with those of Frederiksen and O’Kane [1].

For higher resolution and higher Reynolds numbers we expect that the Markovian closures, like the Eulerian homogeneous non-Markovian DIA [37] and inhomogeneous QDIA [39], will need to incorporate a regularization, or empirical vertex renormalization, in order to yield the correct small-scale spectra. This regularization of the MICs is described in Appendix C of Frederiksen and O’Kane [1]. Briefly, it is achieved by replacing the interaction coefficients $A(\mathbf{k},\mathbf{p},\mathbf{q})$ by $\theta (p-k/\alpha )\theta (q-k/\alpha )A(\mathbf{k},\mathbf{p},\mathbf{q})$ and $K(\mathbf{k},\mathbf{p},\mathbf{q})$ by $\theta (p-k/\alpha )\theta (q-k/\alpha )K(\mathbf{k},\mathbf{p},\mathbf{q})$ in the response function and two-time cumulant equations, but not in the single-time cumulant or mean field equations. Here, $\theta$ is the Heaviside step function and $\alpha$ is a wavenumber cut-off parameter which plays the same role as the $\gamma$ in the eddy damping for the EDQNM closure in Equation (54) of Section 7. As noted in the Introduction, a value of $\alpha \approx 4$ appears to be essentially universal, or only weakly flow dependent, for two-dimensional homogeneous and inhomogeneous turbulence and for three-dimensional HIT.

The performance of all the abridged Markovian inhomogeneous closures studied in this Section are very encouraging, given the very rapidly changing flow evolution in the simulations. This includes the results for the $MI{C}^0[\overline{\zeta }(t)]$ model which uses the current-time mean field and the current-time FDT, such as the EDQNM closure. From the $MI{C}^0[\overline{\zeta }(t)]$ it is possible to develop a generalized EDQNM closure for inhomogeneous turbulent flows, through suitable analytical specifications of generalized eddy damping.

This then removes the need for the auxiliary prognostic equations for the relaxation functions ${\Theta }^0(\mathbf{k},\mathbf{p},\mathbf{q})(t)$ and ${\Psi }^0(\mathbf{k},\mathbf{p})(t)$ and yields a very efficient closure. By analogy, we call this closure the Eddy Damped Markovian Inhomogeneous Closure (EDMIC) that we formulate next.

7. The Eddy Damped Markovian Inhomogeneous Closure

The formulation of the abridged MIC closures, in Section 6, with the current-time mean field approximation ${\overline{\zeta }}_{\mathbf{k}}(s) \equiv <{\zeta }_{\mathbf{k}}(s)>\to <{\zeta }_{\mathbf{k}}(t)>\equiv {\overline{\zeta }}_{\mathbf{k}}(t)$ (Equation (26)) in the time history integrals becomes particularly simple when $X=0$ in Equations (32a–e) and (33a–f). Then, the FDT is the current-time FDT of Equation (2). This is of course the FDT that can be used to derive the EDQNM equations for homogeneous turbulence from the DIA closure [17,18,64,68]. In this Section, we show that from our abridged $MI{C}^{X=0}[\overline{\zeta }(t)]$ it is possible to establish the EDMIC model which is suitable generalization of the EDQNM for homogeneous turbulence to inhomogeneous turbulent flows.

The prognostic equations for the EDMIC model are again Equation (35) for the second order cumulant $C_{\mathbf{k}}(t,t)$ and Equation (44) for the mean field $<{\zeta }_{\mathbf{k}}(t)>$. However, rather than determining the relaxation functions at a given time from the auxiliary prognostic equations in Equation (39b,d) we seek analytical parameterized expressions as for the EDQNM model.

We note first that the solution to the response function differential equation in Equation (36) is

$$R_{\mathbf{k}}(t,t^{\prime })=\exp \left(-{\displaystyle {\int }_{t^{\prime }}^{t}ds{\displaystyle \sum _{j=0}^2{D}_k^j(s)}}\right)$$

(48)

where $D_k^j$ are defined in Equations (33a–f) and (34). Thus, from Equation (39a–d) the relaxation functions with $X=0$ simplify to

$${\Theta }^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{\displaystyle \sum _{j=0}^2\left[D_k^j(s)+D_p^j(s)+D_q^j(s)\right]}}\right)}$$

(49)

and

$${\Psi }^0(\mathbf{k},\mathbf{p})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}d{t}^{\prime }\exp \left(-{\displaystyle {\int }_{t^{\prime }}^{t}ds{\displaystyle \sum _{j=0}^2\left[D_k^j(s)+D_p^j(s)\right]}}\right)}.$$

(50)

Next, we make the Markovian approximation for the $D_{}^j$ such that $D_{}^j(s)\to D_{}^j(t)$; these terms can therefore be taken outside the integrals in Equations (48)–(50) giving

$$R_{\mathbf{k}}(t,t^{\prime })\doteq R_k^{\otimes }(t,t^{\prime })=\exp \left(-{\displaystyle \sum _{j=0}^2{D}_k^j(t)}(t-t^{\prime })\right),$$

(51)

$$\begin{array}{l}{\Theta }^0(\mathbf{k},\mathbf{p},\mathbf{q})(t)\doteq {\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}d{t}^{\prime }\exp \left(-{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)+D_q^j(t)\right]}(t-t^{\prime })\right)}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{1-\exp \left(-{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)+D_q^j(t)\right]}(t-t_0)\right)}{{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)+D_q^j(t)\right]}},\end{array}$$

(52)

and

$$\begin{array}{l}{\Psi }^0(\mathbf{k},\mathbf{p})(t)\doteq {\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)={\displaystyle \underset{t_0}{\overset{t}{\int }}d{t}^{\prime }\exp \left(-{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)\right]}(t-t^{\prime })\right)}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{1-\exp \left(-{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)\right]}(t-t_0)\right)}{{\displaystyle \sum _{j=0}^2\left[D_k^j(t)+D_p^j(t)\right]}}.\end{array}$$

(53)

Here, the superscript symbol $\otimes$ can denote either $EDMIC$ or $EDQNM$.

7.1. Analytical Triad Relaxation Functions for EDQNM

In the EDQNM equations for two-dimensional homogeneous turbulence [64] the mean fields and topography are zero ($<\zeta >\equiv \overline{\zeta }=0;<U>\equiv \overline{U}=0; h=0$) and the only prognostic equation is Equation (35) for the second order cumulant $C_{\mathbf{k}}(t,t)$. This also means that $D_{}^2\equiv D_{}^{\pi }=0$. Moreover, the eddy damping $D_{}^1\equiv D_{}^{\eta }$ in the triad relaxation function in Equation (51) is generally specified by an analytical form that is consistent with the $k^{-3}$ enstrophy cascading inertial range:

$$D_k^1(t)\equiv D_k^{\eta }(t)\to {\mu }_k^{eddy}(t)=\gamma {\left[k^2{C}_k(t,t)\right]}^{\frac{1}{2}}.$$

(54)

Here, $C_{\mathbf{k}}(t,t)$ is real and positive and $\gamma$ is a positive empirically determined dimensionless coefficient. Thus, the EDQNM for homogeneous turbulence has the considerable simplification and computational efficiency of having an analytical expression for the triad relaxation time ${\Theta }^{EDQNM}(\mathbf{k},\mathbf{p},\mathbf{q})(t)$, given in Equation (52), with the superscript $\otimes \to EDQNM$.

For homogeneous turbulence on an $f$-plane, and particularly for HIT, the Rossby wave frequency ${\omega }_k^{}$ vanishes, $D_k^0=\widehat{\nu }{k}^2>0$ is real as is $D_k^1(t)\equiv D_k^{\eta }(t)\to {\mu }_k^{eddy}(t)>0$ in Equation (54). Thus,

$${\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)=\frac{1-\exp \left(-\left[{\mu }_k^{}(t)+{\mu }_p^{}(t)+{\mu }_q^{}(t)\right](t-t_0)\right)}{{\mu }_k^{}(t)+{\mu }_p^{}(t)+{\mu }_q^{}(t)}$$

(55)

where

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

(56)

Now, ${\Theta }^{\otimes }={\Theta }^{EDQNM}(\mathbf{k},\mathbf{p},\mathbf{q})(t)=\mathrm{Re}{\Theta }^{EDQNM}(\mathbf{k},\mathbf{p},\mathbf{q})(t)\ge 0$ is both real and non-negative and this ensures that the cumulant $C_{\mathbf{k}}(t,t)$ is also real and non-negative from Equation (35), and thus realizable, as also shown in Appendix D with superscript $EDMIC\to EDQNM$.

For homogeneous anisotropic turbulence interacting with Rossby waves on a $\beta$-plane $D_k^0=\widehat{\nu }{k}^2+\mathrm{i}{\omega }_k$ and again taking $D_k^1(t)\equiv D_k^{\eta }(t)\to {\mu }_k^{eddy}(t)$ we have

$$R_{\mathbf{k}}(t,t^{\prime })\doteq R_k^{\otimes }(t,t^{\prime })=\exp \left(-[{\mu }_k^{}(t)+\mathrm{i}{\omega }_k](t-t^{\prime })\right),$$

(57)

and

$${\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)=\frac{1-\exp \left(-\left[{\mu }_k^{}(t)+{\mu }_p^{}(t)+{\mu }_q^{}(t)+\mathrm{i}({\omega }_k+{\omega }_p+{\omega }_q)\right](t-t_0)\right)}{{\mu }_k^{}(t)+{\mu }_p^{}(t)+{\mu }_q^{}(t)+\mathrm{i}({\omega }_k+{\omega }_p+{\omega }_q)}.$$

(58)

Unfortunately, as noted by Bowman et al. [68], in the presence of waves it cannot always be guaranteed that $\mathrm{Re}{\Theta }^{\otimes }=\mathrm{Re}{\Theta }^{EDQNM}(\mathbf{k},\mathbf{p},\mathbf{q})(t)\ge 0$ and so there may be situations where $C_{\mathbf{k}}(t,t)$ is no longer realizable.

In the presence of waves, it has been customary to employ the steady state form of the triad relaxation time

$$\mathrm{Re}{\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(\infty )=\frac{{\mu }_k+{\mu }_p+{\mu }_q}{{({\mu }_k+{\mu }_p+{\mu }_q)}^2+{({\omega }_k+{\omega }_p+{\omega }_q)}^2}$$

(59)

with $\otimes =EDQNM$. This has been the approach in a number of studies of turbulence interacting with Rossby waves on a $\beta$-plane [17,75] or on a sphere [76] or interacting with internal gravity waves [18,77]. Then, again $\mathrm{Re}{\Theta }^{\otimes }=\mathrm{Re}{\Theta }^{EDQNM}\ge 0$ and $C_{\mathbf{k}}(t,t)\ge 0$ is realizable.

The EDQNM closure reduces to the Boltzmann equation in the resonant wave interaction limit [78], also called wave-turbulence [22,79], as the damping vanishes. Thus,

$$\begin{array}{c}\lim \\ \mu \to 0\end{array}\mathrm{Re}{\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(\infty )=\pi \delta ({\omega }_k+{\omega }_p+{\omega }_q)$$

(60)

where, by $\mu \to 0$ is implied that all damping terms ${\mu }_k$ vanish, and $\otimes =EDQNM$. The resonant interaction limit was discussed by Holloway and Hendershott [75], Carnevale and Martin [17], and Carnevale and Frederiksen [18] and reviewed by Newell and Rumpf [79] and Sagaut and Cambon [22]. As noted by Carnevale and Frederiksen [18], the Boltzmann equation, to which the EDQNM model reduces in the limit in Equation (60), has an additional conservation law and the limit is singular.

7.2. Analytical Relaxation Functions for EDMIC

For the EDMIC model the essential difference from the EDQNM is the fact that $D_k^2(t)\equiv D_k^{\pi }(t)$ needs to be parameterized as well as $D_k^1(t)\equiv D_k^{\eta }(t)$ and there is the relaxation function ${\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)$ in addition to ${\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)$ where $\otimes =EDMIC$. If these terms are successfully represented by analytical expressions, then there are no additional obstacles to efficient closure from having the mean field prognostic equation for $<{\zeta }_{\mathbf{k}}(t)>$ in Equation (44) in addition to Equation (35) for $C_{\mathbf{k}}(t,t)$. In the presence of mean field and topography one would again expect to represent $D_k^1(t)\equiv D_k^{\eta }(t)\to {\mu }_k^{eddy}(t)$ by the enstrophy cascading inertial range form in Equation (54) for two-dimensional and quasigeostrophic (QG) turbulence. In contrast, for $D_k^2(t)\equiv D_k^{\pi }(t)$ a wider choice of parameterizations may be of interest depending on the form and strength of the mean field and topography. This is a topic that we shall study in some detail in a sequel to this paper. However, there are some general points that can be made about the parameterization of $D_k^2(t)\equiv D_k^{\pi }(t)$.

In the enstrophy cascading inertial range of typical atmospheric and oceanic turbulent flows [57,80] the mean energy spectra fall off much faster than the transients. Thus,

$$D_k^1(t)+D_k^2(t)\approx {\mu }_k^{eddy}(t)=\gamma {\left[k^2{C}_k(t,t)\right]}^{\frac{1}{2}},$$

(61)

with ${\mu }_k^{}(t)=\widehat{\nu }{k}^2+{\mu }_k^{eddy}(t)$ given in Equation (56), may be quite a good approximation in the inertial range. This may also be the case more generally if the mean field and topography are relatively weak or if forcing functions are strong outside the inertial range. In such situations we again have the same four cases to consider as for the EDQNM.

For turbulent flow on an $f$-plane ${\Theta }^{EDMIC}(\mathbf{k},\mathbf{p},\mathbf{q})(t)$ is given by the right-hand side of Equation (52) and $\mathrm{Re}{\Theta }^{EDMIC}\ge 0$ is both real and non-negative. Similarly,

$${\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)=\frac{1-\exp \left(-\left[{\mu }_k^{}(t)+{\mu }_p^{}(t)\right](t-t_0)\right)}{{\mu }_k^{}(t)+{\mu }_p^{}(t)}$$

(62)

with $\otimes =EDMIC$ and $\mathrm{Re}{\Psi }^{EDMIC}\ge 0$ is both real and non-negative. This ensures the realizability of the EDMIC model.

For turbulent flow interacting with Rossby waves and topography on a $\beta$-plane, $R_k^{\otimes }(t,t^{\prime })$ is again given by Equation (57), ${\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)$ by Equation (58) and

$${\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)=\frac{1-\exp \left(-\left[{\mu }_k^{}(t)+{\mu }_p^{}(t)+\mathrm{i}({\omega }_k+{\omega }_p)\right](t-t_0)\right)}{{\mu }_k^{}(t)+{\mu }_p^{}(t)+\mathrm{i}({\omega }_k+{\omega }_p)}$$

(63)

with $\otimes =EDMIC$. As noted above for the EDQNM, and as discussed by Bowman et al. [68], the wave terms mean that $\mathrm{Re}{\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)\ge 0$ and $\mathrm{Re}{\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)\ge 0$ cannot always be guaranteed and so there may be situations where the EDMIC is not realizable. We note that in principle the realizability of the $MI{C}^{X=0}[\overline{\zeta }(t)]$ and $MI{C}^{X=1}[\overline{\zeta }(t)]$ of Section 5 and Section 6, and the $MI{C}^{X=0}[\overline{\zeta }(s)]$ and $MI{C}^{X=1}[\overline{\zeta }(s)]$ of Frederiksen and O’Kane [1], cannot be guaranteed. In practice they performed very similarly to the associated realizable $MI{C}^{X=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(t)]$ and $MI{C}^{X=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(s)]$, and this in situations of very rapidly evolving flows from small amplitude. It should also be noted that the quasi-Lagrangian closures for HIT [41,42,43,44,45], unlike the Eulerian DIA and QDIA, do not guarantee realizability. Nevertheless, it is of course desirable for closures to be manifestly realizable in all situations and in a sequel to this paper we aim to generalize the representations of the time dependent relaxation functions in Equations (58) and (63) to ensure realizability.

Realizability of the EDMIC model is again ensured with the steady state form of the relaxation time in Equation (59) and with

$$\mathrm{Re}{\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(\infty )=\frac{{\mu }_k+{\mu }_p}{{({\mu }_k+{\mu }_p)}^2+{({\omega }_k+{\omega }_p)}^2}$$

(64)

where $\mathrm{Re}{\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(t)\ge 0$ and $\mathrm{Re}{\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(t)\ge 0$ and $\otimes =EDMIC$.

The EDMIC model equations can also be developed for the resonant interaction limit where ${\Theta }^{\otimes }(\mathbf{k},\mathbf{p},\mathbf{q})(\infty )$ is given by Equation (59) and

$$\begin{array}{c}\lim \\ \mu \to 0\end{array}\mathrm{Re}{\Psi }^{\otimes }(\mathbf{k},\mathbf{p})(\infty )=\pi \delta ({\omega }_k+{\omega }_p)$$

(65)

with $\otimes =EDMIC$.

The EDMIC closure reflects the fact that the cumulant and response function relationships in Appendix B simplify with the current-time mean field and current-time FDT in the time-history integrals and with the analytical form of the response function in Equation (51).

7.3. EDMIC for Three-Dimensional Turbulent Flows

The EDQNM closure was of course originally formulated for three-dimensional Navier–Stokes HIT by Orszag [59] and implemented for the corresponding two-dimensional HIT by Leith [64]. The EDQNM has subsequently been applied with various modifications to the eddy damping such as integral forms over wavenumbers [81] and including rotation for problems of HIT and homogeneous anisotropic turbulence (HAT) [21,22]. We expect that the EDMIC can, with effort in generalizing the parameterization of the eddy damping, also be extended to three-dimensional inhomogeneous turbulent flows.

One of the most straightforward generalizations of the EDMIC closure equations is to three-dimensional QG inhomogeneous turbulent flows. The QG dynamical equations, with continuous vertical variations, for flow over topography on an $f$-plane, are presented in Appendix B of Frederiksen [50]. The generalization to flow on a $\beta$-plane and including a large-scale flow $U(z)$ which varies in the vertical direction $z$ is straightforward as is the derivation of the consequent EDMIC model. Moreover, the structures of the EDMIC model equations for two-dimensional and three-dimensional QG inhomogeneous turbulent flows are essentially the same or isomorphic.

8. Discussion and Conclusions

We have formulated statistical dynamical closure equations for inhomogeneous turbulent flows interacting with Rossby waves and topography, at several levels of simplification, and tested their performance against large ensembles of direct numerical simulations on a $\beta$-plane. Firstly, the non-Markovian Quasi-diagonal Direct Interaction Approximation (QDIA) closure [48,49], $QDIA[\overline{\zeta }(s)]$, has been abridged to the $QDIA[\overline{\zeta }(t)]$ variant in which the current-time mean field $\overline{\zeta }(t)$ replaces the complete mean field trajectory $\overline{\zeta }(s)$ between the initial and current times, $t_0\le s\le t$, in the time history integrals. Secondly, from the abridged $QDIA[\overline{\zeta }(t)]$ closure, three variants of Markovian Inhomogeneous Closures (MICs) have been formulated based on different versions of the Fluctuation Dissipation Theorem (FDT). These are the current-time FDT, the prior-time FDT and the correlation FDT. The computational cost of the abridged MICs, like the original MICs formulated by Frederiksen and O’Kane [1], scales like $O(T)$ where $T$ is the total integration time. In contrast, the cost for the QDIA closures scales like $O(T^3)$. The MICs do however need to compute auxiliary prognostic equations for the relaxation functions, which replace the time history integrals of the QDIA, and contain similar information. The abridged MICs have the computational advantage of there only being two such relaxation functions while for the original MICs [1] there are three.

The efficacy of the abridged closures in capturing the evolved statistical dynamics of large (1800-member) ensembles of direct numerical simulations (DNS) has been tested in stringent numerical experiments with rapidly developing mean fields on a $\beta$-plane. The numerical experiments start with a large-scale eastward mean flow $U$ impinging on a mid-latitude conical mountain in the northern hemisphere within a turbulent environment and with an initial small-scale mean field of much lesser amplitude. Simulations are performed for 10 days during which the non-zonal streamfunction of the mean field rapidly develops into a large-scale Rossby wavetrain downstream of the mountain and wave-turbulence and eddy-eddy interactions change the transient kinetic energy and palinstrophy wavenumber distribution.

Pattern correlations of the 10-day evolved mean non-zonal streamfunction between the abridged closures and DNS ensemble range from 0.9789 for $QDIA[\overline{\zeta }(t)]$, through 0.9994 for $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(t)]$ and 0.9995 for $MI{C}^1[\overline{\zeta }(t)]$ to 0.9999 for $MI{C}^0[\overline{\zeta }(t)]$. Interestingly, for the original closures of Frederiksen and O’Kane [1] the pattern correlation is identical, at 0.9999, for the $MI{C}^0[\overline{\zeta }(s)]$ and 0.9998 for the $QDIA[\overline{\zeta }(s)]$, $MI{C}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}[\overline{\zeta }(s)]$, and $MI{C}^1[\overline{\zeta }(s)]$. Thus, the QDIA is most sensitive to using the current-time mean field and the MIC that also employs the current-time FDT is least sensitive. However, the mean field results are remarkably good in all cases, abridged or original. This also extends to the transient kinetic energy and palinstrophy spectra where there are just slight differences with DNS near the peak at wavenumber $k=4$.

The fact that the abridged closures perform so well even when the mean flow is rapidly evolving indicates that the perturbation fields are also rapidly decorrelating. This suggests that rapid Rossby wave growth is associated with rapid error growth and loss of deterministic predictability. Indeed, this agrees with the finding of Frederiksen [82] that instabilities tend to grow fastest when storms and blocks intensify. It also agrees with the results of ensemble weather forecasts where errors tend amplify when dynamical development of Rossby waves is fastest, as shown in Figure 9 of Frederiksen et al. [83] and further discussed by O’Kane and Frederiksen [52].

The robustness of the performance of the inhomogeneous closures suggest that it may be possible to replace the auxiliary prognostic equations for the relaxation functions by analytic expressions as in the Eddy Damped Quasi Normal Markovian (EDQNM) for homogeneous turbulence. We demonstrate that the abridged $MI{C}^{X=0}[\overline{\zeta }(t)]$ model, with both the current-time mean field and current-time FDT, can be adapted to an Eddy Damped Markovian Inhomogeneous Closure (EDMIC) consisting of a mean field equation and single-time cumulant equation with analytical expressions for the relaxation functions. We suggest that the EDMIC is the natural generalization to inhomogeneous turbulence of the EDQNM including its very high computational efficiency. We note that the EDMIC model is realizable under the same conditions as the EDQNM.

In a sequel to this work, we plan to study the performance of the EDMIC model, and its dependence on the eddy damping parameters, and we aim to include wave renormalization effects that ensure realizability in the presence of time dependent waves for both the EDMIC and EDQNM models.