## Abstract

Recently developed closure-based and stochastic model approaches to subgrid-scale modelling of eddy interactions are reviewed. It is shown how statistical dynamical closure models can be used to self-consistently calculate the eddy damping and stochastic backscatter parameters, required in large eddy simulations (LESs), from higher resolution simulations. A closely related direct stochastic modelling scheme that is more generally applicable to complex models is then described and applied to LESs of quasi-geostrophic turbulence of the atmosphere and oceans. The fundamental differences between atmospheric and oceanic LESs, which are related to the difference in the deformation scales in the two classes of flows, are discussed. It is noted that a stochastic approach may be crucial when baroclinic instability is inadequately resolved. Finally, inhomogeneous closure theory is applied to the complex problem of flow over topography; it is shown that it can be used to understand the successes and limitations of currently used heuristic schemes and to provide a basis for further developments in the future.

## 1. Introduction

We review recent progress in the development and applications of dynamical subgrid-scale parametrizations for geophysical flows based on statistical dynamical and stochastic methods. The focus is on approaches where the subgrid model is determined self-consistently from the statistics of higher resolution closure calculations or direct numerical simulations (DNSs). These methods differ from more traditional methods in that no tuning parameters are used in the lower resolution large eddy simulation (LES). We also briefly summarize the historical development and application of dynamical subgrid modelling for geophysical flows.

### (a) Deterministic parametrizations for atmospheric flows

Since the first climate simulations with atmospheric general circulation models (GCMs), it has been recognized that the specifications of subgrid processes play major roles in determining the accuracy of the large-scale flows and energy spectra [1]. The Smagorinsky [1] empirical subgrid model for eddy viscosity has been one of the most celebrated methods for parametrizing eddy viscosity. It is more appropriate for three-dimensional turbulence than for quasi-geostrophic turbulence where various Laplacian or higher-order dissipation schemes are generally used. The sensitivity of the large-scale kinetic energy to resolution and to the strength and form of the diffusion parametrizations, as found for example by Manabe *et al*. [2], has been a continuing issue with atmospheric GCMs [3,4,5,6]. The same qualitative changes in energy spectra with varying horizontal resolution were found in early studies with the barotropic vorticity equation and were analysed on the sphere using canonical equilibrium theory by Frederiksen & Sawford [7]. Frederiksen *et al*. [8] proposed an explanation for this sensitivity, on the basis of ideas of energy and enstrophy conservation and statistical mechanics.

### (b) Deterministic parametrizations for oceanic flows

In oceanic GCMs, the need to parametrize the effects of the subgrid-scale eddies is perhaps even more important because baroclinic instability for the formation of mesoscale eddies occurs at much smaller scales that are generally not resolvable in climate simulations. The early approach was also to use empirically specified diffusion operators in the horizontal and vertical [9]. However, Redi [10] proposed a diffusion operator that mixes tracers along isopycnals rather than horizontally and vertically. Gent & McWilliams [11] proposed a parametrization scheme for the skew-symmetric, or advective, eddy flux tensor for passive tracers, temperature and salinity. McDougall & McIntosh [12] developed a formalism applicable to transport of passive tracers, temperature and salinity in three-dimensional time-averaged oceanic flows and obtained a similar form to the Gent–McWilliams parametrization. As yet, there is no analogous approach for the momentum equations. In these approaches, the retained scale fields are taken to be the mean fields, and the subgrid-scale eddies are the transient fluctuations about the mean [13]. A review of this literature is given by Griffies *et al*. [14].

### (c) Statistical dynamical subgrid modelling

Kraichnan [15] used his direct interaction approximation (DIA) to show that for isotropic turbulence the inertial transfer of energy could be represented as a combination of both an eddy viscous and stochastic backscatter terms. These subgrid processes are associated with the respective eddy-damping and nonlinear noise terms that constitute the right-hand side of the DIA tendency equation [16]. Leith [17] used the eddy-damped quasi-normal Markovian (EDQNM) closure to calculate an eddy dissipation function that would preserve a stationary *k*^{−3} kinetic energy spectrum for two-dimensional turbulence. Kraichnan [18] developed the theory of eddy viscosity in two and three dimensions and identified the existence of a strong cusp in the spectral eddy viscosity near the cut-off wavenumber representing local interactions between modes below and near the cusp. Further studies pertaining mainly to three-dimensional homogeneous turbulence have been carried out by Leslie & Quarini [19], Chollet & Lesieur [20], Carnevale & Frederiksen [21], Leith [22], Chasnov [23], Domaradzki *et al.* [24], Herring [25] and McComb and coworkers [26,27,28].

Frederiksen & Davies [29], hereafter FD1997, developed representations of eddy viscosity and stochastic backscatter based on EDQNM and DIA closure models for barotropic turbulent flows on the sphere. They found that their parametrizations cured the typical resolution dependence of atmospheric energy spectra with LES being in close agreement with higher resolution barotropic DNS. Subsequent implementation of their eddy viscosity parametrization into an atmospheric GCM resulted in significantly improved circulation and energy spectra [6].

Frederiksen [30], hereafter F1999, derived general expressions for the eddy–topographic force, eddy viscosity and stochastic backscatter based on the quasi-diagonal direct interaction approximation (QDIA). O'Kane & Frederiksen [31] calculated QDIA-based subgrid-scale parametrizations (SSPs) considering observed atmospheric flows over global topography and quantifying the relative importance of the subgrid-scale eddy–topographic, eddy–mean field, quadratic mean and mean field–topographic terms. O'Kane & Frederiksen [31] and Frederiksen & O'Kane [32] also compared the QDIA-based SSPs with heuristic approaches, based on maximum entropy, that have been used to improve systematic deficiencies in ocean climate models [33].

A crucial question for more general flows is how to systematically parametrize each of the respective subgrid-scale drain and injection terms identified using closure methods in atmospheric and oceanic GCMs. Frederiksen [34] has recently formulated a multi-field generalization of the computationally tractable barotropic QDIA inhomogeneous closure. The multi-field QDIA was also formulated for LESs and included subgrid-scale parametrizations that ensure the same large-scale statistical behaviour as higher resolution closures. These closures were also shown to have underpinning non-Markovian stochastic models that guarantee realizability for the elements of the covariance matrices that are diagonal in spectral space. The multi-field QDIA was formulated for studying inhomogeneous turbulence and subgrid-scale parametrizations for geophysical flows described by the baroclinic quasi-geostrophic and primitive equations [35] and for three-dimensional flows [36]. However, the framework is generally applicable to classical field theories with quadratic nonlinearity.

### (d) Stochastic subgrid modelling

It has also been realized that subgrid-scale eddies not only drain energy and enstrophy from the retained scales through an effective eddy viscosity but also inject energy and enstrophy to the retained scales through a stochastic backscatter term. Leith [22] showed that empirical representations of stochastic backscatter, in combination with the Smagorinsky [1] eddy viscosity, could improve LES in a plane shear mixing layer.

While closure models may be the natural starting place for developing subgrid-scale parametrizations, their complexity makes them difficult to formulate and apply to multi-field models such as GCMs. Instead, a direct stochastic modelling approach to subgrid processes based on the statistics of DNS, such as used by Frederiksen & Kepert [37], may be a more tractable methodology. Stochastic modelling has had a long and successful history of application to problems in atmospheric science [38,39]. Recently, there has been increasing interest in exploring how parametrizations of stochastic backscatter may improve simulations and predictions of weather and climate [40,41,42]. Shutts [43], O'Kane & Frederiksen [44] and Berner *et al.* [45] have considered the role of stochastic backscatter in weather prediction studies where it increases ensemble spread and improves forecasts in a manner similar to ensemble methods [46]. Berner *et al*. [47] also found that stochastic backscatter represented through a cellular automaton (as in Shutts [43]) resulted in reductions in systematic errors and improvements in seasonal forecasts. Seiffert & von Storch [48] found that model climate sensitivity to increased carbon dioxide concentrations depends on whether a stochastic backscatter parametrization is used.

Zidikheri & Frederiksen [49] generalized the methodology of Frederiksen & Kepert [37] to the two-level quasi-geostrophic model. In this formulation, vertical inhomogeneities are represented by the off-diagonal elements of generalized damping and random forcing covariance matrices. They applied the formulation to atmospheric flow with sheared zonal jets, Rossby waves and baroclinic instability. They found that LESs using these parametrizations successfully reproduced the statistics of higher resolution DNSs. Furthermore, at sufficiently high resolutions of the LESs, the off-diagonal elements of the aforementioned matrices became less important in magnitude compared with the diagonal elements. Recently, Kitsios *et al*. [50,51] showed that universal scaling properties can be obtained for these parametrizations, making them applicable to a wide range of flows.

The two-level quasi-geostrophic model can also be applied to ocean circulations. In this case, however, the situation is quite different from that in the atmosphere because the deformation radius in ocean flows is typically less than 100 km. Therefore, the deformation scale is not fully resolved in typical ocean climate simulations. This implies that such models do not correctly capture the conversion of large-scale potential energy to smaller-scale kinetic energy. Another effect that is not captured is the injection of energy from the subgrid scales to the resolved scales as predicted by the phenomenology of quasi-geostrophic turbulence [52,53]. This contributes to the variability of the large scales, and indeed also to the organization of large-scale mean currents. This is not captured by models that are not eddy-resolving, or only use deterministic dissipation parametrizations of the subgrid scales. Zidikheri & Frederiksen [54,55] applied the stochastic modelling approach to oceanic LES and demonstrated that it can be successfully used to maintain the correct (high-resolution) spectra for simple flows and even for relatively complex simulations of the Antarctic circumpolar current (ACC).

### (e) Plan of paper

The plan of this paper is as follows. In §2, we summarize the baroclinic quasi-geostrophic flow equations for two-dimensional and Rossby wave turbulence on a *β*-plane and on the sphere, both in physical space and in spectral space. Section 3 presents the general form of non-Markovian Eulerian closures for inhomogeneous turbulence in spectral space. The subgrid model expressions for eddy viscosities and stochastic backscatter and renormalized drain and dissipation matrices are presented in §4. In §5, comparisons of LES using these closure-based subgrid-scale parametrizations with higher resolution DNS are presented for barotropic atmospheric flows with typical kinetic energy spectra. The fact that ‘statistical closures such as the DIA and QDIA have underpinning Langevin equations from which they can be constructed’ has motivated direct stochastic approaches to the subgrid-scale parametrization problem. Section 6 reviews a direct stochastic subgrid modelling approach that accounts for memory effects of turbulent eddies. Here, the subgrid terms and consequent LESs obtained with the stochastic method are also compared with closure-based results for barotropic atmospheric flow of §5. Applications of the stochastic approach to subgrid modelling and LESs for typical baroclinic atmospheric flows, for simple baroclinic oceanic flows and for baroclinic flows with broad similarities to the ACC are also reviewed. In §7, we discuss the structure of QDIA-based subgrid-scale parametrizations for a canonical equilibrium state and for a non-equilibrium state corresponding to the atmospheric flow for January 1979.

## 2. Baroclinic quasi-geostrophic flow over topography

Taking suitable length and time scales, the non-dimensional equation for two-level baroclinic quasi-geostrophic flow over topography on a *β*-plane and on the sphere may be written in the form
2.1
Here, *a*=1 or 2, *ψ*^{a} is the streamfunction and *q*^{a}=∇^{2}*ψ*^{a}+(−1)^{a}*F*_{L}(*ψ*^{1}−*ψ*^{2}) is the reduced potential vorticity, *h*^{2}=*h*, *h*^{1}=0, where *h* is the scaled topography, *F*_{L} is the layer coupling parameter, are bare forcing functions and is the generalized bare dissipation that may include the Rossby wave frequency.

The corresponding spectral space equations are given by
2.2
where *f*^{a}(**k**,*t*) is the bare forcing, and *D*^{ab}(**k**) is the generalized bare dissipation that may include the Rossby wave frequency [34]. Also
2.3a
2.3b
2.3c
Here, *δ*(**k**,**p**,**q**) is unity if **k**+** p**+

**q**=0 and is otherwise zero, and the interaction coefficients are defined in equation (9) of Frederiksen [34] for the planar geometry case. On the sphere, the interaction coefficients are detailed in Frederiksen & Kepert [37]; the same equations apply with the equivalences

**k**=(

*m*,

*n*), −

**k**=(−

*m*,

*n*) and

*k*=

*n*(

*m*is zonal and

*n*is total wavenumber) and selection rules as in Frederiksen & Sawford [7]. Throughout, summation over repeated superscripts is implied.

## 3. Closures and subgrid-scale terms for inhomogeneous geophysical flows

In this section, we review the QDIA closure theory and present analytical expressions for the various subgrid-scale terms entering the equations for the mean and fluctuating components of the vorticity within the quasi-geostrophic equations for inhomogeneous turbulent flow over topography.

### (a) The QDIA closure equations

The QDIA closure formulated by Frederiksen [30,34] provides a general theoretical framework for the direct calculation of the subgrid-scale parametrizations of the eddy–topographic force, eddy viscosity and stochastic backscatter of non-equilibrium flows interacting with inhomogeneous turbulence and topography. O'Kane & Frederiksen [56], hereafter OF2004, implemented the QDIA, using an efficient restart method allowing for non-Gaussian initial conditions and approximation to higher-order vertex terms [57,58]. The theory and numerics were further extended by Frederiksen & O'Kane [59], hereafter FO2005, to incorporate differential rotation and Rossby waves.

We now briefly state the closure equations required for a discussion of the subgrid calculations and refer the reader to Frederiksen [34] for more details. Let us suppose we have an ensemble of flows satisfying equation (2.2) where and 〈 〉 denotes mean and denotes fluctuation about the mean. Also, where .

The mean-field closure equation takes the form 3.1 The nonlinear damping measures the interaction of transient eddies with the mean field 3.2 while 3.3 measures the strength of the interaction of transient eddies with the topography. Here 3.4a and 3.4b

The equation for the diagonal two-time two-point cumulant is
3.5a
Here, and the nonlinear noise is
3.5b
with
3.5c
and
3.5d
being the noise and dissipation terms associated with eddy–mean field and eddy–topographic interactions. The diagonal response function satisfies
3.5e
where and *δ*^{ab} is the Kronecker delta function. The single-time cumulant may be obtained from the expression
3.6
Canonical equilibrium states [60,61] are exact solutions of the inviscid unforced QDIA closure. The QDIA reduces to the DIA [15,62] for zero mean field and topography.

## 4. Subgrid terms

Frederiksen [30,34] derived the expressions for the subgrid scales based on the QDIA equations and the associated Langevin equations. Here, we consider the case where the resolution is reduced from T to R<T with R the resolution of the resolved scales. In planar geometry, we define *T*={**p**,**q**|**p**≤T,**q**≤T}, *R*={**p**,**q**|**p**≤R,**q**≤R} and *S*=*T*−*R* is the set of the subgrid scales. We now simply state the dynamical equations for the resolved scale vorticity including subgrid terms. We define the renormalized mean force as
4.1
where, with *S* denoting subgrid scales, the residual Jacobian is given by
4.2
and the subgrid eddy–topographic force by
4.3
Similarly, we may define the renormalized random forcing variance by
4.4
4.5
4.6
Here, is the backscatter and is the bare random forcing variance. We also define renormalized drain matrices respectively by
4.7
and
4.8
where
4.9
and
4.10
The backscatter dissipation matrix elements are then defined by
4.11
and net and renormalized net dissipation matrix elements respectively by
4.12
and
4.13
The corresponding generalized viscosity matrix elements are then defined by
4.14
and
4.15
where k^{2} is *k*^{2} on the plane and *n*(*n*+1) on the sphere, and • denotes any of 0,*b*,*d*,*n*,*r*,*rn*.

## 5. Comparisons of direct numerical simulation and large eddy simulation with closure-based subgrid models

In this section, we review the performance of LESs, for two-dimensional isotropic and Rossby wave turbulence, incorporating closure-based subgrid-scale parametrizations as described by FD1997. Comparisons are made of the spectra of DNS with LES at various resolutions and with both stochastic and deterministic subgrid-scale parametrizations. FD1997 focused on DNSs, which reproduced the observed January 1979 isotropized kinetic energy spectrum as a statistically stationary spectrum in the presence of random forcing and dissipation. The DNS runs have been performed at triangular T63 resolution, with the dissipation taken as a linear combination of surface drag and a Laplacian dissipation for which the non-dimensional bare viscosity is specified in equation (4.1) of FD1997. The DNS simulations are also forced by random forcing *F*_{0}(**k**) specified in fig. 2 of FD1997.

The bare viscosity and bare forcing when used in DNS of the barotropic vorticity equation also produce a statistically steady state that is very close to the January 1979 spectrum. This is shown in figure 1*a* (adapted from FD1997), which, at T63, compares the results of a DNS run averaged over the last 100 days of a 150-day simulation with the January 1979 spectrum. The results in figure 1 are for Rossby wave turbulence with the *β*-effect due to differential rotation *B*=2. Very similar results are obtained for isotropic turbulence with *B*=0 (fig. 1 of FD1997). Figure 1*a* shows the dimensional total wavenumber kinetic energy spectra *E*(*n*). Also shown are these spectra ± the standard deviations *σ*(*n*). The DNS spectra in figure 1*a*, or, almost equivalently, the T63 January 1979 spectra, are regarded as the benchmark or truth against which LES at lower resolution are to be compared.

The various subgrid-scale terms are shown in figure 2*a* (adapted from fig. 3*b* of FD1997) for the case when the initial T63 January 1979 spectrum is truncated back to T31. Figure 2*a* depicts the eddy viscosities *ν*_{d}(*n*), *ν*_{b}(*n*) and *ν*_{n}(*n*). These eddy viscosities have been calculated for isotropic barotropic turbulence, using the EDQNM closure. This closure may be derived from the DIA closure as described in §3 of FD1997; very similar results are also obtained from the DIA closure as shown in fig. 9 of FD1997. We note that these subgrid-scale parametrizations all have a cusp behaviour at the smallest scales. Fig. 3 of FD1997 shows the qualitative similar subgrid-scale parametrizations at T15.

Figure 1*a* also compares the spectra of LES using the renormalized viscosity *ν*_{r}(*n*)=*ν*_{d}(*n*)+*ν*_{0}(*n*) and renormalized random forcing *F*_{r}(*n*)=*F*_{b}(*n*)+*F*_{0}(*n*) at the resolution T31. Again, the kinetic energy spectra shown are averaged over the last 100 days of 150-day simulations. The diagrams also show the spectra ± the standard deviation *σ*(*n*) and the energy spectra of the T63 DNS, and January 1979 observations, truncated back to T31. The agreement between LES and DNS or January 1979 kinetic energies is good at the smaller scales; at the larger scales, there are some deviations, presumably because of the few *m* components to average over and the long eddy turn-over times.

Frederiksen & Davies [29] show that the generally good agreement between DNS and LES using the EDQNM-based renormalized viscosity and renormalized noise forcing is also found if the EDQNM-based renormalized net viscosity is used instead or if the corresponding very similar DIA-based parametrizations are used. They also find that the comparison between DNS and LES is much better than when a number of ad hoc viscosity parametrizations are used in the LES.

## 6. Stochastic modelling approach

The stochastic modelling approach is underpinned by the stochastic model representation of the DIA and QDIA equations, as described by Frederiksen [30,34]. The advantage of the stochastic modelling approach is that it can readily be generalized to more complex models without explicitly formulating the complex closure equations. Furthermore, the damping and noise covariance coefficients may be determined by DNSs. In what follows, we consider the stochastic modelling approach as it applies to quasi-geostrophic turbulence on the sphere described by the quasi-geostrophic equations.

### (a) Quasi-geostrophic two-level model

The quasi-geostrophic two-level model for baroclinic dynamics [63] is discretized by expanding each of the fields in spherical harmonics with zonal wavenumber *m* and total wavenumber *n* to yield the spectral equations (2.2). Here, we specify the form, in equations (2.2) where **k**=(*m*,*n*). This is a generalized complex bare dissipation whose real part describes dissipation, and imaginary part describes the frequency of Rossby waves; *ρ* is a positive integer discussed below. These equations are non-dimensional; we use the Earth's radius as a length scale and the inverse of the Earth's angular velocity as a time scale. The beta effect *B*, due to differential rotation, takes the value 2, with current scalings, but we also consider the case *B*=0. The layer coupling parameter *F*_{L} is inversely proportional to the potential temperature difference between the levels. We also allow for a simple representation of the effects of differential heating through a specification of . This is the value towards which is relaxed on a time scale given by *κ*^{−1}. The barotropic model may also be obtained from the baroclinic equations by setting *F*_{L}=0.

### (b) Subgrid-scale parametrizations

Frederiksen & Kepert [37] formulated a stochastic subgrid modelling method that takes into account memory effects of turbulent eddies by analogy with the DIA closure approach of Frederiksen & Davies [29] and the QDIA closure of Frederiksen [30]. Here, we briefly summarize the approach.

Let **q** denote the column vector of spectral coefficients . Then, on reducing the resolution from *T* (for DNS) to *R* (for LES), we write the tendency as **q**_{t}=(**q**_{t})_{R}+(**q**_{t})_{S}, where **q**_{t}=(∂**q**(*t*)/∂*t*). The vector (**q**_{t})_{R} represents the resolved scale tendency and (**q**_{t})_{S} represents the subgrid-scale tendency. Now, **q** consists of a transient part and a mean part . Similarly, the subgrid tendency has both mean and transient parts , where and are the mean and transient parts, respectively. In this study, the mean part is simply calculated from a high-resolution simulation as the time average of the subgrid tendency.

The subgrid terms are represented by the stochastic equation
6.1
Here, ** D_{d}** is defined as the drain dissipation matrix, and is a random forcing vector. The memory effects of turbulent eddies are accounted for and, in particular, the drain dissipation matrix has the integral form
6.2
Here, angular brackets denote either ensemble or time averaging; in this study, we consider the latter; also † denotes Hermitian conjugate.

The covariance matrix for the subgrid nonlinear noise , where may be obtained from the Lyapunov or balance equation
6.3
after computing **D _{d}**.

In general, the nonlinear noise may be coloured, but it is found that, once it has been calculated as described earlier, it is sufficient to model it by white noise in the subsequent LES. Thus, the forcing (t) takes the white noise form where *δ*(*t*−*t*′) is the Dirac delta function.

We also consider a deterministic parametrization where is represented by a backscatter dissipation matrix defined by 6.4 Then, we define the net dissipation matrix by 6.5

### (c) Atmospheric barotropic model results

In order to compare the effectiveness of the stochastic modelling methodology for deriving dynamical subgrid-scale parametrizations from DNS, with the closure-based results of §5, we again consider the isotropized monthly averaged spectrum for January 1979. We consider the case of barotropic turbulence with no mean field or mean forcing, as in the study of Frederiksen & Kepert [37], and then go on to more complex flows later in this section. As in §5, the DNS is carried out at the resolution of triangular truncation T63, and the flow is maintained by the same random forcing and the bare dissipation, *D*_{0}(*n*)=*ν*_{0}(*n*)*n*(*n*+1).

Figure 1*b* (upper results) shows the kinetic energy *E*(*n*), averaged over *m* and averaged over the last 100 days of a 150-day DNS with the barotropic vorticity equation, the standard deviations about the means and the January 1979 spectrum. These simulations and the following LES are for Rossby wave turbulence with *B*=2. As in §5, we regard the T63 DNS spectra as the benchmarks against which we compare the LES at a lower resolution, chosen as T31, with dynamical subgrid-scale parametrizations. These subgrid-scale parametrizations are calculated as described in §6*b*. Because there is only one level and because we consider isotropized matrices for comparison with §6, **D** and **F**_{b} of §6*b* become scalars *D* and *F*_{b} that depend only on the total wavenumber *n*.

Figure 2*b* shows the viscosities *ν*_{d}=*D*_{d}[*n*(*n*+1)]^{−1}, *ν*_{b}=*D*_{b}[*n*(*n*+1)]^{−1} and *ν*_{n}=*D*_{n}[*n*(*n*+1)]^{−1}=*ν*_{d}+*ν*_{b}, which depend only on the total wavenumber *n*. These viscosities have a cusp behaviour at the smallest scales.

Comparisons of the DNS-based quantities with the corresponding EDQNM-based results in figure 2, and DIA-based results in fig. 9 of FD1997, indicate a close correspondence in the shape of the parametrizations (apart from the jaggedness of the DNS-based results) but with some differences in magnitude. We expect that these differences between the DNS-, DIA- and EDQNM-based results are partly due to sampling with DNS (leading to jaggedness in the parametrizations).

LESs have been performed at T31 with the renormalized drain dissipation *D*_{r}=*D*_{d}+*D*_{0} and white noise renormalized forcing and also with the renormalized net dissipation *D*_{rn}=*D*_{n}+*D*_{0} and bare forcing . Figure 2*b* (middle results) shows the LES kinetic energy spectra *E*(*n*), averaged over the last 100 days of a 150-day simulation with *D*_{r} and , as well as the standard deviations about these spectra. Also shown are the energy spectra of the DNS truncated back to T31. The lower results in figure 2*b* are the corresponding findings with *D*_{rn} and .

Both sets of LES results in figure 2*b* are in very good agreement with DNS truncated back to T31. They are very comparable to the corresponding EDQNM-based results in figure 2*a* and possibly slightly better than for the DIA-based results in fig. 11 of FD1997.

### (d) Atmospheric baroclinic model results

Next, we review the performance of the stochastic subgrid modelling method for more realistic atmospheric flows, with large-scale mean jets and differential rotation (*B*=2), on the basis of the baroclinic model study of Zidikheri & Frederiksen [49]. The parameters used are described there.

The DNS was performed at T126 using the quasi-geostrophic model to reproduce typical atmospheric jets and energy spectra. The layer coupling constant *F*_{L} = 100 in non-dimensional units corresponding to a deformation radius of about 500 km.

The model was integrated to statistical steady state and then further stepped forward in time for 208 days. Figure 3 shows the level 1 kinetic energy, *E*(*n*), of the DNS truncated back to T63; at smaller scales, it follows an *n*^{−3} power law.

Again, the stochastic modelling approach was used to calculate the dynamical subgrid-scale parametrizations for LES at T63. Successful LESs result from parametrizations that are inhomogeneous in the vertical and homogeneous in the horizontal. This is consistent with the quasi-diagonal turbulence closure theory [30,34] and its ability to simulate the statistics of DNS in the presence of mean flow inhomogeneities [44]. The eddy dissipation and nonlinear noise terms then take the forms *D*^{jk}(*m*,*n*) and , where *j*,*k*=1,2; that is, they reduce to 2×2 matrices at each wavevector (*m*,*n*). As for the barotropic problem of §6*c*, the parametrizations are essentially isotropic at the smallest scales where they have a cusp (see fig. 3 of Zidikheri & Frederiksen [49]). Note that the mean forcing, , is insignificant in magnitude because the baroclinic injection is fully resolved at resolution T63.

LESs were performed at T63 for three different specifications of the dissipation. Firstly, the same bare dissipation (and forcing) as used for the DNS was used. Secondly, the renormalized net dissipation, which is the sum of the net eddy dissipation and the bare dissipation, was used. Finally, LES with renormalized drain dissipation, which is the sum of the eddy drain dissipation and the bare dissipation, and stochastic backscatter was performed. The level 1 kinetic energy spectra *E*(*n*) of the three LESs are redrawn in figure 3. The top curves show the LES with bare parameters; note that the tail of the LES kinetic energy spectrum lifts due to insufficient drain of potential enstrophy. The middle curves show the LES with the deterministic renormalized net dissipation, and the bottom curves show the LES with renormalized drain dissipation and stochastic backscatter. Both parametrizations perform very well, with LES in excellent agreement with higher resolution DNS at T126. Very similar agreement was obtained for the level 2 kinetic energy and for the potential energy (not shown).

### (e) Oceanic baroclinic model results

In oceans, the internal Rossby radius of deformation is at least an order of magnitude smaller than in the atmosphere. This implies that the injection scale corresponds to wavenumbers in the vicinity of 100. If we wish to run an LES at a lower resolution, such as T63, then the energy injection region will not be fully resolved in the LES. Clearly, then, a purely dissipative parametrization cannot be even roughly correct in such a simulation, in contrast to the case of eddy-resolving atmospheric simulations. According to the phenomenology of two-layer quasi-geostrophic turbulence, energy is injected in the vicinity of the deformation scale, and is then cascaded towards the large scales via the barotropic mode [52,53]. This implies that there is a net energy injection into the barotropic modes at the large scales, which would appear as a negative dissipation coefficient if the deterministic formulation is pursued.

#### (i) Simple flow

Zidikheri & Frederiksen [54] considered subgrid modelling in the case of the simple flow configuration described by Salmon [52], including when the Rossby radius of deformation is typical of oceanic flows. In this type of flow, the turbulence is driven by a purely baroclinic large-scale flow that results in the average kinetic energy and enstrophy spectra at the two levels being identical. The zonal current is relaxed towards a solid-body rotation profile , where *ϕ* is the latitude, and are the maximum currents (at the equator). These currents are specified by and . The layer coupling constant *F*_{L}=10 000 in non-dimensional units and corresponds to a deformation radius of about 50 km. As in §6*c,* we set *B*=0. The resolution of the DNS is T252 with 768×384 grid points (in longitude × latitude), while the subgrid elements were calculated for LESs with a resolution of T63.

The model was integrated to statistical steady state and then further stepped forward in time for 104 days for deriving statistics of the DNS of the turbulent flow. The isotropized elements for all subgrid terms are shown in fig. 11 of Zidikheri & Frederiksen [54].

In principle, it is also necessary to account for the conversion of mean potential energy into transient kinetic energy, through . However, with the mean field confined to the (*m*,*n*)=(0,1) mode, this term is negligible. Again, the performance of the LES using these parametrizations, when compared with higher resolution simulation is excellent, as shown in figure 4. This is also the case for the potential energy spectra (not shown).

#### (ii) Antarctic circumpolar current

The ACC of the Southern Ocean provides a further, more realistic, example of the application of subgrid-scale parametrizations to ocean flows. Zidikheri & Frederiksen [55] studied the performance of the stochastic subgrid modelling method when applied to simulated flows broadly resembling the ACC. They used forcing and dissipation parameters so that the statistically steady state has a zonal mean current centred at 60^{°}S with values broadly resembling those found in the ACC. This DNS at T252, corresponding to a grid size of approximately 0.5^{°}, will be regarded as the benchmark for comparison with the LES in subsequent sections.

In figure 4*b*, we show comparisons of the kinetic energy spectra of LES at T63, with and without the stochastic subgrid-scale parametrizations, and the DNS at T252 truncated back to T63. The subgrid parametrizations of dissipation and stochastic backscatter are shown in figs 5–8 of Zidikheri & Frederiksen [55] while is shown in their fig. 3. The LES with the subgrid-scale parametrizations is clearly in far better agreement with the DNS than the LES with bare dissipation. LESs with deterministic subgrid-scale parametrizations have also been studied for oceanic-type turbulent flows but, depending on the flow parameters, the purely deterministic net dissipation may, on occasions, lead to numerical instability when baroclinic instability is inadequately resolved.

## 7. Subgrid modelling with QDIA closure for global atmospheric flow

In this section, QDIA-closure-based subgrid-scale parametrizations are shown for calculations [31] both for a canonical equilibrium case (including forcing and dissipation balanced at each wavenumber) and for a non-equilibrium dissipative case corresponding to observed January 1979 flow over topography. A crucial question is the relative strengths of the residual Jacobian, subgrid-scale eddy–mean field and eddy–topographic interactions for flows far from equilibrium. In the following discussion, the subgrid-scale parametrizations are calculated for circularly truncated wavenumber *R*=24 LESs based on the results of higher resolution *T*=48 DNSs.

### (a) Canonical equilibrium and maximum entropy

For the canonical equilibrium case, the integral terms appearing in the single-time cumulant equation in equations (3.5e) and (3.6) are real and homogeneous, whereas the mean-field subgrid terms in equation (3.1) are complex and highly anisotropic. All closure integral terms balance with their counterparts so that canonical equilibrium is exactly maintained. Figure 5 clearly shows that, consistent with the F1999 statistical theory, the eddy–topographic force is balanced exactly by the damping due to eddy–mean field interactions and in physical space (figure 6*a*) resembles a high-pass-filtered form of the topography. Crucially, for the canonical equilibrium state, the residual Jacobian is identically zero. Energy injection from the subgrid scales to the retained scales not only is seen to occur in the form of stochastic backscatter as a result of eddy–eddy interactions but also arises due to nonlinear noise terms associated with eddy–mean field and eddy–topographic interactions (figure 5). These forcings balance the drain of energy from the retained to subgrid scales due to eddy–eddy, eddy–mean field and eddy–topographic dampings. The spectra of the QDIA-based subgrid terms, including the renormalized and net eddy viscosities and eddy–topographic force, are all cusp-shaped at the cut-off scale (figure 5). These results are consistent with corresponding subgrid terms for isotropic turbulence based on DIA and EDQNM closures and on stochastic modelling, as discussed in §1. For the current inhomogeneous case, this again results in the subgrid terms being of small scale and more specifically that the subgrid-scale eddy–topographic force will be a high-pass-filtered version of the topography. Figure 6*b* depicts the strong correlation of the zonally asymmetric component of the streamfunction and the mean vorticity with the topography. We have found that qualitatively similar results hold for the corresponding inviscid unforced canonical equilibrium case (not shown).

The numerical results for subgrid-scale parametrizations described in this section confirm the analytical findings and general principles of the F1999 statistical theory.

### (b) Dissipative January case

For general non-equilibrium flows, there are unfortunately no analytic solutions and numerics must be our guide. Again, we refer to the study of O'Kane & Frederiksen [31], who considered the case of a typical 500 mb January atmospheric flow over topography. They found that the subgrid contributions in the single-time cumulant equation are in general complex, prior to taking the real parts, and an imbalance between the energy injection and drain terms occurs at each wavenumber. In addition, a similar imbalance occurs between eddy–topographic, eddy–mean field and eddy–eddy interactions contributing to the mean flow. Importantly, unlike the canonical equilibrium case, the residual Jacobian not only does not vanish but also makes a contribution comparable in magnitude to both the damping due to eddy–mean field interaction and to the eddy–topographic force. Large differences exist between zonally asymmetric components of the mean streamfunctions for equilibrium and non-equilibrium flows, respectively. In figure 6*b*, we see that at equilibrium the zonally asymmetric streamfunction is strongly correlated with the contours of the topography, whereas the non-equilibrium case, in figure 6*d*, has a more complex structure with Rossby wavetrains in the mid-latitudes of the Northern Hemisphere in the band of the strong zonal jets.

Figure 6*a* shows the structure of the eddy–topographic force for the equilibrium case with strong interactions occurring over the large-scale topographic structures such as the Himalayas, the Andes, Greenland and to a lesser extent the Rocky Mountains; the non-equilibrium case is remarkably similar (not shown). For the non-equilibrium case, the residual Jacobian (figure 6*c*) has significant contributions from both the mean field–mean field interaction and the mean field–topographic interaction. In addition, the quadratic mean–field terms in the residual Jacobian consist of small-scale Rossby wavetrains stretching from the Himalayas over the North Pacific, whereas the mean-field topographic terms have features broadly similar, but of larger scale, to those seen in the eddy–topographic force.

## 8. Conclusion

We have reviewed recent progress in the development and application of dynamical subgrid models needed for LESs of geophysical flows based on statistical dynamical and stochastic methods. Our focus has been on approaches where the subgrid models are determined self-consistently from the statistics of higher resolution simulations rather than specified through tuning parameters.

We have discussed studies showing that closure models can be used to calculate the eddy dissipation and stochastic backscatter parameters needed in LESs of homogeneous barotropic turbulence so as to maintain the spectra of high-resolution simulations. Some of these works have generalized the closures to spherical geometry and applied them to improving the simulations of atmospheric GCMs. The statistical closures that we have discussed have underpinning Langevin equations from which they can be reconstructed. This suggests a direct stochastic modelling approach that does not require the explicit construction of closure models. We have discussed such a stochastic modelling approach that gives similar results for the eddy dissipation and stochastic backscatter as the closure theory for barotropic flows. It is also able to successfully maintain the spectra of high-resolution baroclinic simulations when applied to LESs.

We have distinguished subgrid-scale parametrizations in atmospheric and oceanic models. In atmospheric models, the deformation scale, which is related to the scale of turbulent energy injection, is quite large and can be easily resolved by numerical models. As a result, deterministic schemes, such as hyperviscosity models that drain enstrophy near the truncation scale, have some justification. In ocean models, however, the deformation scale is over an order of magnitude smaller and is therefore typically not resolved in long climate simulations. In that case, deterministic schemes may not always be adequate as there is a net injection of energy from the subgrid scales onto the resolved scales. We have demonstrated that this injection can be modelled by a combination of drain dissipation and stochastic backscatter.

We have discussed the application of the QDIA closure theory to subgrid modelling for flow over topography. At canonical equilibrium, the QDIA-based eddy–topographic force is a high-pass-filtered form of the topography. For non-equilibrium flows, the picture is much more complex. Again, the QDIA-based eddy–topographic force is a high-pass-filtered version of the topography. However, the residual Jacobian, which vanishes at canonical equilibrium, now makes an important contribution to the mean flow forcing of comparable magnitude to the eddy–topographic force and the eddy–mean field interactions. Moreover, the residual Jacobian is inhomogeneous, containing structures correlated to both the topography and the quadratic mean-field term. Although mathematically complex, the natural manner in which the closure subgrid-scale terms arise and the transparent physical interpretation of the component terms are remarkable.

The complexity of the closure subgrid models required for self-consistency also makes clear the enormously difficult problem presented in trying to develop subgrid-scale parametrizations for general geophysical flows using heuristic approaches. We expect that, in the future, subgrid modelling will see important advances through the development of closures for inhomogeneous baroclinic turbulence and general mean flows and corresponding stochastic modelling motivated by closure theory.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales I’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.