## Abstract

Parametrizations of the subgrid eddy–eddy and eddy–meanfield interactions are developed for the simulation of baroclinic ocean circulations representative of an idealized Antarctic Circumpolar Current. Benchmark simulations are generated using a spectral spherical harmonic quasi-geostrophic model with maximum truncation wavenumber of *T*=504, which is equivalent to a resolution of 0.24^{°} globally. A stochastic parametrization is used for the eddy–eddy interactions, and a linear deterministic parametrization for the eddy–meanfield interactions. The parametrization coefficients are determined from the statistics of benchmark simulations truncated back to the large eddy simulation (LES) truncation wavenumber, *T*_{R}<*T*. A stochastic technique is used to determine the eddy–eddy coefficients, and a new least-squares regression method for the eddy–meanfield terms. Truncations are repeated for various *T*_{R}, and the resolution dependence of the subgrid coefficients is identified. The mean jet structure and the kinetic and potential energy spectra resulting from the LESs closely agree with those from the benchmark simulations.

## 1. Introduction

It is not possible to simulate the global ocean circulation by explicitly resolving all of the scales of motion. Instead one resorts to large eddy simulation (LES), where the large-scale eddies are resolved on a computational grid and the interactions between the resolved scales and the unresolved subgrid scales are parametrized. There are four major classes of subgrid interactions [1]: eddy–eddy, between subgrid and resolved eddies; eddy–meanfield, between subgrid eddies and the resolved meanfield; meanfield–meanfield, between the subgrid and resolved meanfield; and eddy–topographic, between the subgrid eddies and resolved topography. The first three classes are the focus of this study.

Various approaches have been taken to parametrize these subgrid interactions. The empirical Smagorinsky model [2], applicable for three-dimensional turbulence, is one of the most widely adopted. In quasi-two-dimensional simulations, the subgrid scales are typically represented either by Laplacian or higher order dissipation schemes, or by a modified version of the Smagorinsky scheme [3]. Regardless of the approach, if the subgrid interactions are not properly parametrized, then an increase in resolution will not necessarily increase the accuracy of the explicitly resolved scales. The dependence of the resolved planetary and synoptic scales on resolution has been a problem since the earliest simulations of weather and climate [2,4,5], and persists in even the most sophisticated general circulation models (GCMs). Simulations of the atmosphere [6] and ocean [7,8] have illustrated that the large scales continue to depend sensitively on resolution.

The subgrid parametrization in oceanic simulations is especially important as, in contrast to the atmosphere, at typical grid resolutions baroclinic instability is not explicitly resolved. In practice subgrid parametrizations of the ocean are oriented in isopycnal directions, with the diagonal elements of the dissipation operator represented by the scheme in [9]. The Gent & McWilliams-type [10] parametrizations are typically used to represent the off-diagonal elements, for temperature, salinity and tracers. At present, there is no associated approach for the momentum equations. The strengths of these parametrizations are tuned to yield desirable physical properties and numerical stability [11]. In the derivation of the scheme in [10], the resolved scales are defined to be the meanfield, and the subgrid scales defined to be the fluctuations—as is characteristic of eddy–meanfield parametrizations. One could then consider that this class of parametrizations represent the eddy–meanfield interactions; however, in practice the dissipation operator is applied to both the meanfield and the transients.

As it is only possible to parametrize the statistical effects of the subgrid eddies, statistical dynamical closure theory is the natural formulation for developing self-consistent subgrid models. Next, we provide a brief overview of statistical dynamical closure theory. In the pioneering study of Kraichnan [12], the parameter-free Eulerian direct interaction approximation (DIA) was developed for isotropic homogeneous turbulence. Related closures were developed by Herring [13] and McComb [14]. Self-consistent subgrid models were developed by Frederiksen & Davies [15], based on the DIA and the eddy damped quasi-normal Markovian closures in spherical geometry for barotropic homogeneous turbulence, which achieved excellent statistical agreement between the LES and higher resolution benchmark simulations. Following this, self-consistent statistical dynamical subgrid parametrizations were developed by Frederiksen [16] for inhomogeneous turbulent flows based on a quasi-diagonal DIA (QDIA) closure. The coefficients relating the subgrid tendency to the associated resolved quantities (mean field, fluctuating field or topography) are referred to by Frederiksen [1] by the term self-energy (SE). Generally speaking, an SE is a term that represents the effect on a given mode of the nonlinear interactions of the other modes in the system [17,18]. In this case, these are the unresolved nonlinear processes.

Broadening the applicability of the above closure methods, a stochastic subgrid modelling approach was developed by Frederiksen & Kepert [19], in which the subgrid coefficients are determined from the statistics of higher resolution benchmark simulations. This approach is essentially a numerical means of determining the SE closure coefficients, and referred to here as the SE subgrid modelling approach. This technique has been successfully applied in quasi-geostrophic (QG) simulations of the atmosphere and ocean, comprising sheared jets, Rossby waves and baroclinic instability [20,21]. The resulting coarser grid LESs using the subgrid coefficients produce kinetic energy spectra in agreement with the associated benchmark simulations. For LESs truncated within inertial ranges, the eddy–eddy subgrid coefficients have been shown to be governed by simple functions of resolution [22,23]. The main focus of this study is to develop such scaling functions for the eddy–meanfield term, in addition to presenting the eddy–eddy coefficients.

Specifically, we use the SE subgrid modelling approach to characterize the subgrid interactions for flows broadly representative of the upper layers of the Antarctic Circumpolar Current (ACC). In §2, we summarize the numerical model, with the flows characterized in §3. The subgrid modelling process is outlined in §4, with the eddy–eddy and eddy–meanfield coefficients discussed in §§5 and 6, respectively. The scaling laws presented here are applicable to ocean simulations in an eddy-permitting regime or at higher resolution. LESs are compared to the benchmark simulations in §7, with concluding remarks made in §8.

## 2. Quasi-geostrophic equations

For the present simulations, we employ a two-level QG model [24] without topography, which captures the essential dynamics of baroclinic and barotropic instabilities. The vorticity is represented at two discrete vertical levels, with *j*=1 representing the upper level centred at a depth of 200 m, and *j*=2 the lower level centred at 600 m. The extent of the computational domain starts at the ocean surface (0 m depth) and ends at a depth of 800 m. It has been shown in real-world observations and in numerical simulations that the largest fluctuations in the Southern Ocean occur in approximately the first 800 m [7,25–27], which is why the baseline simulation was designed to have this vertical domain. All variables are assumed to be non-dimensionalized by the Earth's radius (*a*=6371 km) and angular velocity (*Ω*=7.292×10^{−5} s^{−1}) unless units are specified.

In physical space, the two-level QG equations solve for the potential vorticity, *q*^{j}=*ζ*^{j}+(−1)^{j}*F*_{L}(*ψ*^{1}−*ψ*^{2}), where *ζ*^{j} is the vorticity and *ψ*^{j} the streamfunction, all at level *j*. The layer coupling coefficient, *F*_{L}, is inversely proportional to the temperature difference between the two levels, and related to the Rossby radius by . The evolution of *q*^{j} is given by
2.1where the field variables are functions of time (*t*), longitude (*λ*) and , with *ϕ* the latitude. The coefficient *B* represents the beta effect, and the function *J*(*ψ*^{j},*q*^{j}) is the Jacobian. The bare dissipation, , by definition represents the unresolved eddy–eddy interactions in the benchmark simulation. The term *α*^{j} is a Rayleigh drag used to dampen the large scales of motion. Simulations are nudged towards a target climate by relaxation parameter *κ*.

In this study, we solve (2.1) spectrally by expanding the field variables in spectral spherical harmonics of zonal wavenumber *m*, total wavenumber *n*, with meridional wavenumber *n*−*m*. The spectral potential vorticity is
2.2where are the spectral coefficients of the vorticity, with being the streamfunction coefficients, and −*n*(*n*+1) the discrete form of the Laplacian operator. The evolution of is given by
2.3The summations immediately after the equals sign are over the wavenumber set
2.4with *T* the triangular wavenumber. The interaction coefficients are detailed in [19]. The Rossby wave frequency is *ω*_{mn}=−*Bm*/[*n*(*n*+1)], where *B*=2 under the chosen non-dimensionalization. The bare eddy viscosity, , where *δ*_{lj} is the Kronecker delta function, ensuring the off-diagonal elements are zero. Here, is the value of the diagonal elements at truncation and the exponent determines the steepness of . The drag at each level is given by , where erf is the error function, and *n*_{c} is the wavenumber at which . All simulations are driven towards a mean state that is purely zonal (zero unless *m*=0) and corresponds to a large-scale easterly current representative of an idealized ACC. Simulations are driven towards this target state by relaxation parameter, *κ*_{n}, which is only non-zero for *m*=0 and *n*≤15.

## 3. Characterization of the benchmark oceanic flows

Temporal integrations of (2.3) produce the statistics needed to determine the subgrid eddy viscosities. Initial simulations have *κ*_{n}=10^{−6} s^{−1}, , with *n*_{c}=50. The layer coupling parameter is specified as *F*_{L}=2.5×10^{−10} m^{−2}, which corresponds to a Rossby radius of deformation of , and a non-dimensional Rossby wavenumber of *k*_{R}=*a*/*r*_{R}=142. We also test the sensitivity of the subgrid parametrizations to the specification of the Rossby radius, climate relaxation and drag. The highest resolution benchmark simulations within have maximum zonal and total truncation wavenumbers of *T*=504, equivalent to 1536×768 grid points (in longitude × latitude), or a grid point at every *Θ*=120/*T*_{R}=0.24^{°}. Baroclinic instability is resolved in all benchmark simulations with *k*_{R}<*T*.

The dominant structures are located in the mid to high latitudes of the Southern Hemisphere. This is illustrated by the instantaneous level 1 eddy streamfunction field in figure 1*a*. The corresponding level 1 and level 2 time-averaged zonal velocity fields () are shown as functions of latitude in figure 1*b*, along with the climate states that they are nudged towards (). The maximum velocities of the time-averaged current at a depth of 200 m (level 1) and 600 m (level 2) are 0.6 and 0.2 ms^{−1}, respectively. This is consistent with the ACC measurements in [28].

The kinetic energy spectra illustrate the relative energy across all scales. The mean and eddy energy components of the kinetic energy spectra (summed over the *m* wavenumbers) for a two-level simulation are given by
3.1and
3.2where the superscript * is the complex conjugate operation. The total eddy kinetic energy is given by , and the total mean kinetic energy is given by , both of which are is illustrated in figure 1*c*. The mean kinetic energy dominates in the low wavenumber region (large scales), and the transient energy dominates in the higher wavenumber region (smaller scales). The extent of the energy containing scales is *k*_{E}=70 (labelled on the horizontal axis of figure 1*c*), past which one expects self-similar behaviour. The means of determining *k*_{E} is outlined in [23].

Enstrophy flux is the rate at which enstrophy is transferred from one wavenumber to the next. The potential enstrophy flux in level space is defined as
3.3
where
3.4is the potential enstrophy transfer from level *k* into level *j* at wavenumber *n*. In barotropic/baroclinic vorticity coordinates, the enstrophy flux is given by
3.5where is the barotropic/baroclinic vorticity enstrophy transfer. The index 1 refers to the barotropic mode, and 2 the baroclinic mode. For each wavenumber, in matrix form is given by , where the superscript ^{T} denotes the transpose operation, and the transformation matrix
3.6with *c*_{n}=1+2*F*_{L}/[*n*(*n*+1)]. In figure 1*d*, we present the level-based potential enstrophy fluxes *η*^{11}(*n*), *η*^{22}(*n*) and the baroclinic enstrophy flux . The constant values of the enstrophy fluxes within the inertial ranges are denoted by , and , respectively. These quantities are important for the scaling laws of our subgrid parametrizations.

## 4. Stochastic self-energy modelling of the subgrid scales

Here, we present the stochastic SE modelling approach for self-consistently determining the eddy–eddy and eddy–meanfield coefficients from the statistics of two-level QG benchmark simulations. LES resolution is lower than the associated benchmark simulation and confined to the resolved scale wavenumber set
4.1where *T*_{R} is the LES truncation wavenumber such that *T*_{R}<*T*. The subgrid wavenumber set is defined as **S**=**T**−**R**. To facilitate the following discussion, we let for a given wavenumber pair. In this vector notation, the time derivative of **q** can be decomposed such that **q**_{t}(*t*)=**q**^{R}_{t}(*t*)+**q**^{S}_{t}(*t*), where is the resolved tendency involving only triadic interactions between wavenumbers less than or equal to *T*_{R}, and is the subgrid tendency with at least one wavenumber involved in the triadic interactions greater than *T*_{R}. The subgrid tendency is further decomposed as , where is its time average representing the sum of the eddy–meanfield and meanfield–meanfield interactions, and is its fluctuating component representative of the eddy–eddy interactions.

### (a) Subgrid eddy–eddy parametrization

We adopt the approach developed in [19] to represent . It is modelled by the stochastic equation
4.2where **D**_{d} is the subgrid drain dissipation matrix, is the fluctuating component of **q** and is a random forcing vector. As the present simulations have two vertical levels, **D**_{d} is a 2×2 matrix, and is a 2 element column vector. We determine **D**_{d} by post-multiplying both sides of (4.2) by , integrating over the turbulent decorrelation time *τ*, and ensemble averaging. The expression for **D**_{d} is based on a generalization of Gauss's theorem of least squares and given by
4.3where † denotes the Hermitian conjugate for vectors and matrices. The angled brackets denote ensemble averaging, with each ensemble member determined by shifting the initial time *t*_{0} and the final time *t*=*t*_{0}+*τ* forward by one time step. The turbulence decorrelation time, *τ*, is chosen sufficiently large to capture the memory effects of the subgrid turbulence. The model for is determined by calculating the nonlinear noise covariance matrix , where . Again post-multiplying both sides of (4.2) by , and adding the conjugate transpose of (4.2) pre-multiplied by yields
4.4which can be rearranged for . In the implementation of the subgrid model it is sufficient to assume that can be represented as the white noise process , as is witnessed by the agreement between the benchmark simulations and stochastic LESs presented in §7. An eigenvalue decomposition of is used to produce a stochastic representation for [20]. The backscatter dissipation is then given by . One can also represent the subgrid interactions in a simplified deterministic form , where **D**_{net} is given by the sum of the drain and backscatter, representing the net effect.

### (b) Subgrid eddy–meanfield parametrization

We parametrize the eddy–meanfield interactions, using a new least-squares approach to capture the relationship between the ensemble averaged subgrid tendency and the ensemble averaged field . For each wavenumber pair,
4.5where is a 2×2 dissipation matrix representing the eddy–meanfield interactions, and **b** is a 2 element vector of constant coefficients representing the meanfield–meanfield interactions. We assume that (4.2) also holds for small perturbations of the climate centred at the ensemble averaged climate, such that
4.6where and are the time averaged meanfield and subgrid tendency calculated over the *i*th non-overlapping time window of length *τ*_{M}, and *ϵ*_{i} is the associated 2 element error vector. The ensemble averages of each of the terms are , and 〈*ϵ*_{i}〉=0. The dissipation is solved for in a least-squares sense by subtracting (4.2) from (4.3), post-multiplying by , ensemble averaging both sides, and rearranging for to produce
4.7where we assume that the error term *ϵ*_{i} is uncorrelated with . Once is known we can determine the offset by rearranging (4.2) for **b**. The dissipation can be converted to an eddy viscosity via , and likewise for the eddy–eddy operators.

## 5. Eddy viscosities representing the subgrid eddy–eddy interactions

We now present the eddy–eddy coefficients determined from a benchmark simulation with *T*=504 truncated back to various reduced resolutions. The time-step size used is *Δt*=600 s and the eddy–eddy subgrid statistics are accumulated over a period of 6 years. The first truncation is made at *T*_{R}=252 (or *Θ*=0.48^{°}), with a decorrelation time of *τ*=288*Δt*=4 days. The vertical level-based drain dissipation matrix, **D**_{d}, is calculated using (4.3). To facilitate a direct comparison with the eddy–meanfield dissipation to follow, we illustrate the eddy–eddy drain dissipation in barotropic/baroclinic space. The conversion between level and barotropic/baroclinic coordinate systems is given by *ζ*_{B}=**C****q**, where , is the barotropic component, and the baroclinic. The drain dissipation **D**_{d} is transformed into barotropic/baroclinic space via **D**_{d,B}=**CD**_{d}**C**^{−1}, and scaled into eddy viscosity units such that *ν*_{d,B}=**D**_{d,B}/[*n*(*n*+1)].

We present the properties of the lower diagonal element of the eddy–eddy operator *ν*^{22}_{d,B}, to serve as a direct comparison to , which is the dominant term of the eddy–meanfield operator characterized in the following section. The eddy–eddy element *ν*^{22}_{d,B} represents the dissipation of the resolved baroclinic field by the subgrid baroclinic eddies. Note the baroclinic streamfunction is proportional to temperature. The real component of *ν*^{22}_{d,B}(*m*,*n*) is illustrated in figure 2*a*. At this truncation level, *ν*^{22}_{d,B}(*m*,*n*) increases strongly with *n*, has only a weak dependence on *m* for *m*<150, and hence is approximately isotropic for these scales. The upper diagonal element *ν*^{11}_{d,B}(*m*,*n*) has a similar form, and the off-diagonal elements are small in comparison. These observations are also true for the backscatter and net viscosities. At lower truncation levels, the coefficients become more anisotropic, and the off-diagonal elements become proportionally more important [23].

The self-similarity of the eddy viscosities is most clearly illustrated by the isotropized (averaged over *m*) profiles. The real component of the isotropized lower diagonal element of the drain eddy viscosity, *ν*^{22}_{d,B}(*n*), is illustrated in figure 2*b* at various truncation levels for which *T*_{R}>*k*_{E}. As the resolution increases the maximum value decreases. This is consistent with more scales being resolved by the computational grid, and hence less subgrid interactions need to be parametrized. Also at all resolutions, significant values of eddy viscosity are concentrated within the last *k*_{E} wavenumbers before *T*_{R}. This is because the distance in wavenumber space that dominant triadic nonlinear interactions can span is determined by the most energetic large scales [29], which means the significant subgrid interactions involve resolved wavenumbers greater than *T*_{R}−*k*_{E}.

One can characterize these isotropized eddy viscosities by fitting *ν*^{22}_{d,B}(*n*) to
5.1where *ν*^{22}_{d,B}(*T*_{R}) is maximum value at truncation, and the power exponent. Additional simulations were also undertaken to determine the sensitivity of the results to the energy containing scales and Rossby radius. Simulations with modified values of *F*_{L} produced different Rossby wavenumbers of *k*_{R}∈{201,246,284}, with associated Rossby radii of *r*_{R}∈{32,26,22}km. Cases with modified values of *n*_{c} produced results with energy containing scales of *k*_{E}∈{40,50,60}. The power exponents determined from the truncation of all benchmark simulations are plotted against *T*_{R} non-dimensionalized by *k*_{E} in figure 2*d*. With this scaling all of the power exponents truncated outside of the energy containing scales (*T*_{R}>*k*_{E}) collapse to one scaling law increasing with resolution. The slope (or scale selectivity) increases with resolution due to the fact that the width of the eddy viscosities is proportional to *k*_{E} and independent of *T*_{R}. This means that as *T*_{R} increases, the slope increases to confine the significant eddy viscosities to the last *k*_{E} wavenumbers. The maximum value is plotted in figure 2*c*, non-dimensionalized by the length scale 1/*k*_{E}, and the time scale (*η*_{I,B})^{−1/3}, with an additional factor of required to collapse the data [23]. The eddy viscosity maximum value decreases with resolution, as more scales are explicitly resolved. Similar properties concerning both the slope and magnitude are observed for all eddy viscosities in both barotropic/baroclinic space and level coordinates. Scaling laws for the eddy–eddy drain, backscatter and net eddy viscosities in baroclinic and level coordinates are listed in table 1. Note the fundamental form of the net profiles is given by the sum of the drain and the backscatter. Using the properties of the coordinate transforms, the off-diagonal drain matrix elements in level space are approximated by and , and likewise for the backscatter and net [23].

## 6. Eddy viscosities representing the subgrid eddy–meanfield interactions

Determining the eddy–meanfield dissipation () requires very long benchmark simulation integration periods to capture the relationship between the slowly varying climate state () and the associated mean subgrid tendency (). Longer time integrations are required for determining , as opposed to only determining and , because the state space must be sampled multiple times to solve the least-squares regression problem in (4.4). We therefore adopt a lower resolution benchmark simulation with *T*=252 (or *Θ*=0.48^{°}). The time-step size is *Δt*=1200 s and the eddy–meanfield subgrid statistics are accumulated over 100 years.

Firstly, we present the eddy–meanfield coefficients at a truncation level of *T*_{R}=126 (or *Θ*=0.95^{°}). The only mechanism for symmetry breaking in (2.3) is the forcing term . As *κ*_{n} is only non-zero for *m*=0 and *n*≤15, it is only the values of and at these wavenumber components that are significant under sufficient sampling. We therefore consider only these scales of motion in the following analysis. The mean potential vorticity subgrid tendency, , is plotted in physical space in figure 3*a*. When compared to the meanfield in figure 1*b*, it is clear that the amplitude of the mean subgrid tendency is large in regions where the meanfield amplitude is large, and inactive where the meanfield is negligible. Figure 3*a* also indicates that the mean subgrid tendency at level 1 is a mirror image of that at level 2. This feature is also observable from the mean subgrid tendency plotted in spectral space in figure 3*b*, with the mirror image of . This indicates that the mean subgrid tendency modifies the mean shear. It therefore makes more sense to present in baroclinic space. The mean subgrid tendency in barotropic/baroclinic space is given by , where , with and the respective barotropic and baroclinic components. Recall that the flow is governed by the evolution of the potential vorticity which, consistent with its definition in (2.2), can also be represented by . For large *F*_{L}, the flow evolves based on the baroclinic component scaled by effectively *c*_{n}. We then present the mean baroclinic subgrid tendency also scaled by *c*_{n} in figure 3*c*, illustrating that the baroclinic component dominates.

The relationship between the mean subgrid tendency and the resolved meanfield is given by the eddy–meanfield dissipation, . We calculate using (4.4), and find the results are insensitive to *τ*_{M}, for *τ*_{M} greater than one week. The dissipation is then transformed into barotropic/baroclinic space via , with each of the components of plotted in figure 3*d* and the baroclinic components again scaled by *c*_{n}. The lowest wavenumber components (*n*<4) are not included as they have very long eddy turnover times and are not statistically converged. The dominant component, , represents the dissipation of the resolved mean baroclinic field by the subgrid baroclinic eddies.

Subgrid truncations are repeated for various *T*_{R}, with the dominant component illustrated in figure 4*a*. As the system is truncated less heavily (*T*_{R} increasing), there are fewer subgrid eddy–meanfield interactions and consequently approaches zero. For each truncation, we calculate and scale the dissipation into eddy viscosity units by . The dominant component is plotted for various resolutions as the solid lines in figure 4*b*, illustrating that as *T*_{R} increases the required eddy viscosity decreases. This is consistent with the reduction of as resolution (*T*_{R}) increases in figure 4*a*. We also test the sensitivity of the results to the relaxation parameter. In the results discussed thus far *κ*_{n}=10^{−6} s^{−1}. In figure 4*b*, for truncations of *T*_{R}=126, the dashed and dotted lines represent eddy viscosities calculated with *κ*_{n}=5×10^{−7} s^{−1} and *κ*_{n}=2×10^{−7} s^{−1}, respectively. There is an offset between each of the eddy viscosities as they have different baroclinic enstrophy fluxes.

We quantify the magnitude, , and slope, , of the baroclinic term by least-squares fitting the profiles to
6.1for *n*≤15. In addition to the previous cases with differing *κ*_{n}, simulations are also undertaken with a relaxation time of *κ*_{n}=10^{−6} s^{−1} and various energy containing scales *k*_{E}∈{40,50,60}. The magnitude of all of these cases is illustrated in figure 4*c*. This figure illustrates that when non-dimensionalized by the associated baroclinic enstrophy flux () the eddy viscosity magnitude collapses for all flow cases and is hence independent of the energy containing scales (characterized by *k*_{E}) and the climate state (modified by *κ*_{n}). The dependence of the eddy viscosity magnitude on resolution and flow strength is given by the scaling law . This means that as *T*_{R} increases (more scales resolved), the strength decreases at a rate proportional to *T*^{−5}_{R}. This is in contrast to the eddy viscosities representing the subgrid eddy–eddy interactions presented in §5, which decrease in magnitude at a rate of *T*^{−1}_{R}. As the resolution increases the eddy–meanfield interactions, therefore, become proportionally less important. The steepness, , is illustrated in figure 4*d* and maintains a relatively constant value of approximately 2 for all *T*_{R}. It is also independent of the energy containing scales and the background climate state. This is again in contrast to the eddy–eddy coefficients, which have exponents that are dependent upon *k*_{E} and increase strongly with resolution. The eddy–meanfield scaling laws are also listed in table 1.

## 7. Large eddy simulation

LESs adopting the eddy–meanfield and eddy–eddy subgrid parametrizations are compared to the benchmark simulations with *k*_{E}=70 and *k*_{R}=142. The most general anisotropic stochastic LES equations are the same as for the benchmark simulation in (2.3), with
7.1added to the right-hand side, and solved over the wavenumber set **R** instead of **T**. The first two terms on the right-hand side of (7.1) represent the anisotropic stochastic form of the eddy–eddy model. In the anisotropic deterministic variant, the drain dissipation *D*^{jl}_{d}(*m*,*n*) in (7.1) is replaced with the net dissipation *D*^{jl}_{net}(*m*,*n*) and is removed. In the eddy–eddy scaling law variants, the subgrid operators are dependent only on *n* (isotropic), with their functional forms listed in table 1. The third term in (7.1) is the eddy–meanfield model, where is represented by the scaling laws summarized in table 1. The final term in (7.1) represents the meanfield–meanfield interactions; however, **b** is found to be negligible and set to zero in all LESs presented here.

The spectral energies of the benchmark simulation are compared to that of the LES. The kinetic energy at level *j* is given by , where and are defined in (3.1) and (3.2), respectively. The potential energy spectra are given by
7.2The spectra *e*_{1} are illustrated in figure 5*a*, *e*_{2} in figure 5*b* and *e*_{P} in figure 5*c*. The structure of *e*_{P} in figure 5*c* reflects the complex form of the mean jet shear. In each figure, the pairs of spectra are shifted down multiples of one decade for clarity, and labelled by the associated subgrid eddy–eddy parametrization variant and LES truncation wavenumber. The top two pairs of spectra in each of the plots compare the benchmark simulation (dashed line) at a resolution of *T*=252, to LESs (solid line) of resolution *T*_{R}=126 using anisotropic deterministic (aD) and anisotropic stochastic (aS) eddy–eddy parametrizations in conjunction with the eddy–meanfield model. The middle two pairs of spectra compare the benchmark simulation at *T*=504, to anisotropic deterministic and stochastic LES variants at *T*_{R}=252. The bottom two pairs of spectra compare the *T*=504 benchmark simulation to LESs of resolution *T*_{R}=252 this time using scaling law representations for the deterministic (lD) and stochastic (lS) eddy–eddy parametrizations. All LESs closely reproduce the kinetic and potential energy spectra of the benchmark simulations.

The LESs are now compared to the benchmark simulations on the basis of the mean current at each level for the *T*_{R}=252 case. The results of the anisotropic deterministic and stochastic LES variants are illustrated in figure 6*a*, the isotropic variants in figure 6*b* and finally the scaling law parametrizations in figure 6*c*. As is evident from each of the figures, all of the LES variants are able to closely reproduce the mean current from the benchmark simulation. In fact, we find the pattern correlation of the mean current between the benchmark simulation and the LESs is greater than 0.9994 at both levels for all LES variants.

Comparisons are made for LESs with eddy–eddy scaling law parametrizations at other resolutions and Rossby wavenumbers in [23].

## 8. Concluding remarks

Subgrid parametrizations have been developed for oceanic simulations representative of an idealized ACC, in an eddy-permitting regime or at higher resolution. Using the SE modelling framework, subgrid parametrizations for the eddy–eddy and eddy–meanfield interactions were calculated from the statistics of high-resolution benchmark simulations. LESs using these parametrizations closely reproduce the mean jet, and energy spectra of the benchmark simulations.

Scaling laws were developed for the dominant baroclinic component of the eddy–meanfield dissipation, describing how its magnitude and form change with resolution and flow strength. These properties were shown to be independent of the energy containing scales and the background climate state. The magnitude of the eddy–meanfield dissipation decreases with resolution far more rapidly than that of the eddy–eddy parametrization, indicating that as resolution increases the eddy–meanfield terms become proportionally less important. The power exponent of the eddy–meanfield dissipation is found to be independent of wavenumber, and equivalent to Laplacian dissipation. This means the effect subgrid eddies have on the climate state is represented by the Laplacian operator multiplied by the resolved mean field. By contrast, the eddy–eddy power exponents increase with resolution and hence become more scale selective. This means the effect subgrid eddies have on the resolved eddies is represented by the Laplacian operator raised to the appropriate power multiplied by the resolved fluctuating field.

To accentuate these differences, we compare the baroclinic power exponents of the eddy–meanfield dissipation and the eddy–eddy drain, backscatter and net dissipations, evaluated at various resolutions using the scaling laws listed in table 1. For illustrative purposes, the power exponents are divided by 2, rounded down to the nearest integer to present them as effective powers of the Laplacian, and presented in table 2. As the resolution increases, the order of the Laplacian increases strongly for the eddy–eddy dissipation, while the eddy–meanfield dissipation remains Laplacian. Note the exponent of the net dissipation required for the deterministic LES is significantly lower than that of the drain and backscatter required for the stochastic LES.

We have demonstrated the structure of the subgrid parametrizations that are ideally required to undertake simulations with resolution independent statistics. At certain resolutions, this requires high-order Laplacian operators, as outlined in table 2. Such high-order operators can be readily implemented into both spectral and spectral element models. We are currently implementing even higher order dissipation operators in a grid point atmospheric GCM via grid to spectral transforms.

## Funding statement

The authors acknowledge the Commonwealth Scientific and Industrial Research Organisation and the Australian Research Council for funding the research.

## Footnotes

One contribution of 14 to a Theme Issue ‘Stochastic modelling and energy-efficient computing for weather and climate prediction’.

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