## Abstract

Methods motivated by non-equilibrium statistical mechanics of turbulence are applied to solve an important practical problem in geophysical fluid dynamics, namely the parametrization of subgrid-scale eddies needed in large-eddy simulations (LESs). A direct stochastic modelling scheme that is closely related to techniques based on statistical closure theories, but which is more generally applicable to complex models, is employed. Here, we parametrize the effects of baroclinically unstable subgrid-scale eddies in idealized flows with broad similarities to the Antarctic Circumpolar Current of the Southern Ocean. The subgrid model represents the effects of the unresolved eddies through a generalized Langevin equation. The subgrid dissipation and stochastic forcing covariance matrices as well as the mean subgrid forcing required by the LES model are obtained from the statistics of a high resolution direct numerical simulation (DNS). We show that employing these parametrizations leads to LES in close agreement with DNS.

## 1. Introduction

Since the earliest simulations of the climate of non-equilibrium geophysical flows (Smagorinsky 1963) it has been recognized that the effects of unresolved subgrid-scale eddies must be parametrized in order to obtain realistic large scale circulations. Most parametrizations have represented the unresolved processes through heuristic specifications of eddy diffusion (Frederiksen *et al.* 1996), or through arguments based on baroclinic instability (Gent *et al.* 1995) or ideas based on equilibrium statistical mechanics or entropy production hypotheses (Frederiksen & O’Kane 2008).

Leith (1971) and Kraichnan (1976) developed a theoretical foundation for subgrid-scale modelling based on closure theory for homogeneous turbulence. Statistical closure theories of turbulence express the statistics of the flow, at least up to second order, in a closed system of equations. For high Reynolds number flows, such as those encountered in the atmosphere and oceans, developing statistical closures that are accurate at all scales is an extremely difficult problem. The most successful statistical closure theories have been based on renormalized perturbation theory, starting with the direct interaction approximation (DIA) closure of Kraichnan (1959). The subsequent development of statistical closures for homogeneous turbulence has been reviewed by McComb (1990) and Frederiksen & Davies (2004). Recently, closure theories have been generalized to inhomogeneous turbulence interacting with general mean flows (Frederiksen 1999) and applied to problems of flow over topography, Rossby wave dispersion, ensemble prediction and data assimilation (Frederiksen & O’Kane 2008; O’Kane & Frederiksen 2008*b*).

Closure theory is a natural starting point for formulating subgrid-scale eddy parametrization schemes for geophysical flows (Leith 1971; Kraichnan 1976; Carnevale & Martin 1982; Herring 1997) since it is generally only possible to represent the statistical effects of unresolved eddies while their phase relationships with the resolved scales are lost (McComb *et al.* 2001). In this approach, the closure scheme is used to calculate the effective eddy viscosity that is needed to maintain the flow statistics at the lower resolution. In the case of quasi two-dimensional turbulence, energy injected at some scale, for example at the atmospheric synoptic scale due to baroclinic instability, is transferred towards larger scales while enstrophy is transferred towards smaller scales. If the truncation scale is within the region of enstrophy transfer, as it commonly is in atmospheric simulations, the eddy viscosity then parametrizes the drain of enstrophy to the subgrid scales. The eddy viscosity typically has a cusp near the truncation scale due to the predominance of local tranfers. This is also seen in studies of three-dimensional turbulence (Young & McComb 2000), but in quasi two-dimensional turbulence there is also typically a significant negative component to the eddy viscosity at the large scales (Leith 1971; Kraichnan 1976) due to non-local transfers.

Although, statistically, the net transfer is either positive (injection) or negative (dissipation) at a given scale, turbulent eddies dissipate and inject energy and enstrophy at all scales (Rose 1977; Leslie & Quarini 1979). This injection of energy is known as backscatter. Using statistical closure theories, one may not only calculate the (drain) eddy dissipation but also calculate the variance of the stochastic forcing representing the backscatter. This was done by Frederiksen & Davies (1997) for atmospheric barotropic homogeneous flows using eddy damped quasi-normal Markovian (EDQNM) and DIA closure models. Recently, O’Kane & Frederiksen (2008*a*) have calculated eddy dissipation, stochastic backscatter and eddy-topographic force parametrizations for atmospheric inhomogeneous turbulent flows over global topography using the quasi-diagonal direct interaction approximation (QDIA) closure of Frederiksen (1999). Stochastic backscatter leads to the eventual decorrelation between large-eddy simulation (LES) and direct numerical simulation (DNS; Herring 1979, 1997); hence, its inclusion leads to models with more realistic error estimates. This feature of stochastic backscatter has been shown to improve weather and climate simulations of the atmosphere (Palmer *et al.* 2005; Shutts 2005; Berner *et al.* 2008; Kleeman 2008; Seiffert & von Storch 2008).

For complex general circulation models (GCMs) used for weather and climate prediction and climate change prognosis, it is not feasible to formulate and solve the corresponding complicated closure equations. Instead, it is desirable to calculate the subgrid-scale parametrizations from the statistics of DNS with the models as was done by Young & McComb (2000) who calculated the deterministic eddy viscosity in an LES of three-dimensional turbulence from a high-resolution DNS. Frederiksen & Kepert (2006) formulated a stochastic method motivated by closure theory, and showed that it leads to very similar forms for the eddy dissipation and backscatter as obtained using the closure calculations for atmospheric barotropic flows. In this article we apply their DNS-based subgrid modelling scheme to more complex oceanic flows with differential rotation and circumpolar currents. We study the performance of the subgrid-scale parametrizations for idealized currents with some broad similarities to the Antarctic Circumpolar Current (ACC). Zonal currents create small-scale baroclinic instability, a transient energy injection, and differential rotation causes the transfer of energy from the injection scales to the large scales to be anisotropic, with zonally elongated eddies receiving more energy (Frederiksen *et al.* 1996). The flow is also very inhomogeneous because transient activity is centred in the region where the currents are strong. Thus, it provides a good test of the robustness of the methodology.

The plan of this paper is as follows. In §2, we introduce the quasigeostrophic two-level equations for the ocean and in §3, we summarize the DNS-based subgrid-scale parametrization methodology of Frederiksen & Kepert (2006). Results from a high-resolution, eddy-permitting, DNS of idealized oceanic flows with currents broadly resembling the ACC are presented in §4. In §5, we reduce the resolution of the model so that baroclinically unstable eddies are no longer resolved and calculate the eddy dissipation, backscatter covariance, and mean subgrid tendency. These parametrizations are then employed in lower resolution LES and shown to closely maintain the statistics of the high-resolution DNS flow. Our conclusions are summarized in §6.

## 2. The quasigeostrophic equations

We use the two-level quasigeostrophic model of Frederiksen (1998) for simulations and to develop subgrid-scale parametrizations needed for lower resolution LES. The model equations can be written as
2.1where *j* = 1 is the upper level and *j* = 2 is the lower level. The (reduced) potential vorticity at level *j* is defined as
2.2where *ψ*^{j} are level-dependent streamfunctions and *F*_{L} is the layer-coupling parameter. The latter is related to the internal Rossby radius of deformation *r*_{int} = (2*F*_{L})^{−1/2}. The model includes drag specified by *α*^{j} and prescribed viscosity ; *ρ* is a positive integer that describes the order of the Laplacian operator. The coefficient *B* represents the beta effect. The equations are non-dimensional, with the Earth’s radius as length scale and the inverse of Earth’s angular velocity as a time scale; with these scalings *B* = 2. We also allow for a simple representation of the effects of wind stress through a specification of the mean forcing, , and dissipation, *κ*. The parametrized wind stress creates a mean vertical shear in the flow that results in the formation of energetic mesoscale eddies. The mesoscale eddies also send energy back to the large scales; thus, there is a two-way interaction between the large-scale flow and mesoscale eddies.

It is the goal of this study to parametrize this two-way interaction between the large-scale flow and the mesoscale eddies in a dynamically consistent fashion when the resolution is too coarse to resolve the latter. Our aim is to formulate and calculate subgrid-scale parametrizations for baroclinic and barotropic processes unemcumbered by the complexity of boundary and topographic effects (O’Kane & Frederiksen 2008*a*). We focus on the interaction of a baroclinic circumpolar current located at 60° S with baroclinic and barotropic eddies. The simulations are performed on a global aqua-planet; there is very little activity in the regions away from the zonal current, which can be regarded as an idealized model of the ACC. It is again desirable to have a sharp distinction between the resolved and subgrid-scale eddies, and this is achieved by using a spectral representation of the flow fields (Frederiksen & Kepert 2006).

The spectral model is obtained by expanding the fields *ψ*^{j}, *ζ*^{j}, and *q*^{j} in terms of spherical harmonics:
2.3where *χ*(*λ*,*μ*,*t*) represents the above fields on a grid defined by the longitude (*λ*) and the sine of the latitude (*μ*) at time *t*. The fields *χ*(*λ*,*μ*,*t*) can then be transformed into the fields *χ*_{mn}, which are discretized in the horizontal with the zonal wavenumber *m* and total wavenumber *n* describing horizontal degrees of freedom instead of the latitude and longitude; are Legendre functions. Using equation (2.3), the spectral form of the two-level equations becomes
2.4Here the potential vorticity spectral coefficients are given by
2.5where are the vorticity spectral coefficients. Also,
2.6is a generalized complex operator whose real part describes the bare dissipation and imaginary part describes the frequency,
2.7of Rossby waves. The interaction coefficients , with explicitly given by Frederiksen & Kepert (2006), describe the nonlinear coupling between the streamfunction and potential vorticity. In equation (2.4), the resolution of the model is defined by the truncation wavenumber *T*; the wavenumbers *m* and *n* satisfy the relations
2.8which define a triangle in the wavenumber space formed by (*m*,*n*). The total number of the components (waves) in triangular truncation is given by
2.9As in equation (2.8), the summations in equation (2.4) are over wavenumbers in the set
2.10

## 3. Subgrid-scale parametrizations

Next, we summarize the formalism for subgrid-scale parametrizations based on the statistics of higher resolution DNS. Suppose we wish to reduce the resolution of the model to *T*_{R} < *T* in equation (2.4) so that
3.1Because the summations in equation (2.4) extend to wavenumbers beyond *T*_{R}, we have to distinguish those terms that can be explicitly calculated at the reduced resolution, *T*_{R}, the resolved terms, from those that have to parametrized, the subgrid terms. We define R as the set of resolved scales,
3.2in equation (2.4). The set of subgrid scales is defined by
3.3Then on reducing the resolution from **T** (for DNS) to **R** (for LES) we can write the tendency in terms of the retained scale tendency at resolution **R** and the subgrid-scale tendency, denoted by **S**, as follows:
3.4It is also important to distinguish the mean () and transient () parts of the flow. The subgrid tendency will then have corresponding mean and transient parts:
3.5where and are the mean and transient parts of the subgrid tendency, respectively. In this study, we work with steady-state statistics, and the mean part is simply calculated from a high-resolution simulation as the time-average of the subgrid tendency. The transient part of the subgrid tendency is parametrized following the methodology of Frederiksen & Kepert (2006):
3.6Here we have adopted matrix notation so that and are column vectors of order , whose components consist of and , respectively. **D**_{d} is the drain dissipation matrix and is a random forcing vector, the stochastic backscatter vector. Frederiksen & Kepert (2006) formulated the stochastic model for the drain dissipation matrix to account for memory effects of turbulent eddies by analogy with the DIA closure approach of Frederiksen & Davies (1997). Thus, the drain dissipation matrix has the integral form
3.7Here, angular brackets denote ensemble averaging. The ensemble members are generated by successively moving both integration limits forward in time by the same amount, which we take to be one timestep. Thus, for each ensemble member, the time interval of integration, *τ* = *t*−*t*_{0}, is the same. This time interval is chosen so as to capture the decorrelation time of the subgrid contributions. Additionally, † denotes Hermitian conjugate.

The covariance matrix for the subgrid nonlinear noise
3.8where
3.9may be obtained from the Lyapunov or balance equation
3.10after computing **D**_{d}. Here, **T**_{Q} is the subgrid-scale potential enstrophy transfer matrix.

In general the nonlinear noise may be coloured but we find that once it has been calculated as described above it is sufficient to model it by white noise in the subsequent LES. Thus, the forcing takes the white noise form
3.11where *δ*(*t*−*t*′) is the Dirac delta function. The random forcing vector, , may then be constructed from the covariance matrix, **F**_{b}. In general **D**_{d} and **F**_{b} are square matrices of order , and is a vector of the same order. However, we find that successful LESs result from parametrizations that are inhomogeneous in the vertical and homogeneous but anisotropic in the horizontal; this is consistent with quasi-diagonal closure theory (Frederiksen 1999; O’Kane & Frederiksen 2008*a*). The eddy dissipation, **D**_{d}, and nonlinear noise, **F**_{b}, thus reduce to 2×2 matrices and reduces to a two-dimensional vector at each wavevector (*m*,*n*), with components , , and , where *j*,*k* = 1,2.

The LES, employing these parametrizations, may then be represented by the following equation: 3.12Here, the summation is over the (reduced) set R of wavenumbers; , , and are the elements of the drain dissipation matrix, the mean subgrid forcing, and the stochastic backscatter, respectively.

## 4. Simulation of baroclinic flows with mean currents at T252

We begin our studies by performing a long DNS of the quasigeostrophic model described in §2. We choose 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 simulation 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.

The drag, *α*^{j}, has been set to a damping time of 40 days for level 1 (typically at a depth of 1000 m) and 10 days for level 2 (typically at 3000 m); the hyperviscosity m^{4} s^{−1} for both levels; the order of the Laplacian operator *ρ* = 2. The potential vorticity forcing, , has been chosen to correspond to a vertically sheared zonal (*m* = 0) current centred at 60° S. The latitudinal profile of the corresponding zonal acceleration is shown in figure 1*a*. The term in equation (2.4) is only applied when *m* = 0 to control the evolution of the zonal current; the relaxation time, *κ*^{−1}, is 116 days. The layer coupling constant *F*_{L} = 2.4 × 10^{−10} m^{−2} in dimensional units, corresponding to a deformation radius of about 50 km. The resolution is T252 with 768×384 grid points (in longitude × latitude). The model has been integrated to statistical steady state and then further stepped forward in time for 104 days, using a timestep of 600 s, for the purpose of calculating the statistics needed in subgrid-scale parametrizations. The numerics of the model are as described for the barotropic model used by Frederiksen & Kepert (2006).

Figure 1*b* shows the zonally and time-averaged zonal currents at statistical steady state. The mean flow is dominated by the current at 60° S although smaller scale motion does occur elsewhere, especially in the upper layer; the strong current is sheared in the vertical. The presence of mean shear creates baroclinic instability at the small scales. Figure 2 shows the relative distribution of mean and transient kinetic energies as functions of *n*. The large scales are seen to be dominated by the mean field while smaller scales are dominated by transients. The resolution chosen corresponds to a grid size of about 25 km at 60° S and hence the chosen deformation scale (50 km) corresponds to a wavenumber of about 125. Baroclinic instability occurs at wavenumbers slightly smaller than the deformation scale wavenumber and nonlinear processes also transfer energy towards smaller wavenumbers.

## 5. Large-eddy simulations at T63

In this section, we study subgrid-scale parametrizations at lower resolution using the methodology outlined in §3. We choose the truncation wavenumber *n* = 63 (T63), corresponding to a grid size of approximately 2°. In order to calculate the statistics required in equations (3.7) and (3.10) with a large sample, and with reasonable computational cost, we have chosen to implement a two-stage procedure. First, we have used the T252 simulation described in §4 to calculate subgrid-scale parametrizations at the reduced resolution of T126. These parametrizations were then used in an LES at T126, which, being computationally less demanding, was run for the longer time interval of 2700 days, using a timestep of 1200 s. The statistics required for calculating subgrid-scale parametrizations at T63 were obtained using this longer simulation.

As well as calculating the subgrid dissipation and stochastic forcing covariance matrices, it is important to establish the basic nonlinear transfers that take place between the resolved and subgrid scales. Salmon (1978) showed that two-layer quasigeostrophic turbulence may be understood more clearly by decomposing the two-level equations into two modes, namely, barotropic and baroclinic modes. These are defined by and . Here and are the barotropic and baroclinic vorticities, representing the vertically averaged and shear parts of the flow, respectively; and are the vorticities for levels 1 and 2, respectively. Salmon (1978) showed that the barotropic mode behaves essentially like two-dimensional turbulence in that energy injected at some scale is transferred towards the large scales. The baroclinic mode, on the other hand, transfers energy towards the small scales. These ideas are also consistent with the results of Herring (1980, 1988) who similarly found that the large scales are barotropic while small scales are baroclinic in quasigeostrophic turbulence. To check whether the subgrid transfer matrix, (or **T**_{Q} in matrix notation), is consistent with these ideas, we first transform it into the matrix , where **M** is the matrix that transforms the vector into the vector . Hence is the subgrid transfer matrix in terms of barotropic and baroclinic modes; its diagonal elements are related to transfers of barotropic and baroclinic energies between the resolved and subgrid scales.

Figure 3 shows the elements of the matrix , divided by *n*(*n* + 1) to yield the energy transfer matrix, , and averaged over the zonal wavenumber. The (barotropic) diagonal element is positive at the large scales, consistent with the phenomenology of quasigeostrophic turbulence, which predicts transfer of barotropic energy from the injection (subgrid) scales into large (resolved) scales. Additionally, the (baroclinic) diagonal element is negative, consistent with the idea that baroclinic energy is transferred from large (resolved) scales into small (subgrid) scales. There are also significant contributions from the off-diagonal elements of due to vertical coupling. As well as transient transfers of energy, there is a corresponding change in the mean energy, which is represented by the eddy-induced mean forcing field, . In figure 4 we have shown the zonally averaged mean acceleration corresponding to this eddy-induced field in physical space. It is about 10 per cent of the value of the imposed acceleration (figure 1*a*) at 60° S.

We have parametrized the subgrid transfer matrix, **T**_{Q}, in terms of drain dissipation, , and the stochastic forcing covariance matrix, . We have used an integration time interval, *τ*, of 2 days in calculating the drain dissipation in equation (3.7). Larger values of *τ* were not found to change the results significantly. The real parts of the drain dissipation are shown in figure 5, after averaging over the zonal wavenumber. The profiles are scale dependent, with maxima near the truncation scale, typical of subgrid-scale parametrizations, although the scale selectiveness is not as severe as seen in some atmospheric simulations (Frederiksen & Davies 1997; Frederiksen & Kepert 2006). The, generally smaller, corresponding imaginary parts are shown in figure 6. These do not display the typical cusp profiles of the real parts. They are strongly anisotropic (depend on both *m* and *n*) as can be seen in figure 7*b*,*d*. The real parts are also somewhat anisotropic (figure 7*a*,*c*), but to a lesser extent. The elements of (averaged over *m*) are shown in figure 8. The eigenvalues of this matrix, at each wavenumber, are positive as required for it to represent injection of energy via white noise.

We have also performed LES at T63 to check the performance of these parametrizations; the resulting kinetic energy spectra, averaged over zonal wavenumber, are depicted in figure 9. Firstly, we have run the low-resolution simulation without any form of subgrid-scale parametrization (the so-called bare parameter case). In that case the tails of the spectra are seen to rise dramatically, and there is a corresponding drop in large scale energy. These problems are largely overcome by employing the drain dissipation matrix, the stochastic forcing, and the eddy-induced mean forcing in the LES. We see that the T63 LES is then in reasonably good agreement with the higher resolution (T252) simulation. Sampling error is the likely cause of the remaining discrepancies between the LES and DNS spectra.

## 6. Conclusion

In this article we have developed and calculated subgrid-scale parametrizations for unresolved eddies in a quasigeostrophic two-level model with mean flow structure broadly resembling the ACC. We have used the methodology for parametrizing subgrid-scale transients articulated by Frederiksen & Kepert (2006). The method is stochastic, and employs a generalized Langevin equation for the subgrid-scale tendency. The model parameters are in the form of 2×2 dissipation and stochastic forcing covariance matrices at each wavenumber. Underpinning this formulation is the fact that statistical closure theories have stochastic representations. Although it is motivated by closure theories, the method has the advantage of requiring only the statistics obtained from DNS for its implementation; thus it has a wide range of applicability in geophysical flows of varying complexity. Starting with a benchmark T252 (0.5°) simulation, we have truncated the system down to T63 (2°) and calculated the associated drain dissipation and stochastic forcing covariance matrices. The matrices are in general anisotropic in wavenumber space, but when averaged over the zonal wavenumber they have the typical cusp profiles also seen in the simpler barotropic studies mentioned in §1. However, the scale selectiveness of the profiles is not as severe as for typical atmospheric simulations. Physically, the subgrid model parametrizes the net injection of barotropic energy into the resolved scales and drain of baroclinic energy into the subgrid scales. We have also calculated the effective eddy-induced forcing field that modifies the mean field pattern due to the extraction of potential energy by baroclinically unstable eddies. We have found that these parametrizations perform well when implemented in the lower resolution (T63) LES. The parametrizations work at resolutions as low as T15 although these results have not been presented here as even T63 would be considered a coarse resolution for current ocean models. In future studies we hope to implement the methodology in more complex models.

## Acknowledgements

This work was partially supported by CSIRO Complex Systems Science and contributes to the Wealth from Oceans Flagship and to the ACCESS/Climate and Atmosphere Program. A major part of this work is based on a thesis submitted in 2008 by M.J.Z. to the Australian National University for the degree of doctor of philosophy (Zidikheri 2008). M.J.Z. wishes to acknowledge invaluable support and encouragement from Rowena Ball and Bob Dewar of the Australian National University as well as material support from the Australian Research Council, through the Discovery Projects funding scheme (project no. DP0343765). Some of the figures in this article were produced with the assistance of Stacey Osbrough of CSIRO Marine and Atmospheric Research.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society