Theoretical comparison of subgrid turbulence in atmospheric and oceanic quasi-geostrophic models

Due to the massive disparity between the largest and smallest eddies in the atmosphere and ocean, it is not possible to simulate these flows by explicitly resolving all scales on a computational grid. Instead the large scales are explicitly resolved, and the interactions between the unresolved subgrid turbulence and large resolved scales are parameterised. If these interactions are not properly represented then an increase in resolution will not necessarily improve the accuracy of the large scales. This has been a significant and long-standing problem since the earliest climate simulations. Historically subgrid models for the atmosphere and ocean have been developed in isolation, with the structure of each motivated by different physical phenomena. Here we solve the turbulence closure problem by determining the parameterisation coefficients (eddy viscosities) from the subgrid statistics of high-resolution quasi-geostrophic atmospheric and oceanic simulations. These subgrid coefficients are characterised into a set of simple unifying scaling laws, for truncations made within the enstrophy-cascading inertial range. The ocean additionally has an inverse energy cascading range, within which the subgrid model coefficients have different scaling properties. Simulations adopting these scaling laws are shown to reproduce the statistics of the reference benchmark simulations across resolved scales, with orders of magnitude improvement in computational efficiency. This reduction in both resolution dependence and computational effort will improve the efficiency and accuracy of geophysical research and operational activities that require data generated by general circulation models, including weather, seasonal, and climate prediction; transport studies; and understanding natural variability and extreme events.


Introduction
Eddies in the atmosphere and ocean range in size from thousands of kilometres down to the millimetre scale, with energy and enstrophy transferred over these scales via complex non-linear inter-eddy interactions (Kraichnan, 1976).For the numerical simulation of these flows, it is clearly not possible to capture all of these interactions by explicitly resolving the smallest eddies on a computational grid whilst spanning the entire globe.One therefore resorts to large eddy simulation (LES), where the large eddies are resolved on a computational grid, with the interactions between the resolved eddies and the unresolved subgrid eddies represented by an appropriate subgrid turbulence model.If these inter-eddy interactions are not properly represented, then an increase in grid resolution will not necessarily improve the accuracy, which leads to resolution-dependent results (Manabe et al., 1979).This has been a significant problem since the earliest simulations of weather and climate (Smagorinsky, 1963), and persists today in even the most sophisticated general circulation models (GCMs) and research codes (Koshyk and Boer, 1995;Shutts, 2005Shutts, , 2013;;Tennant et al., 2011;Morrison and Hogg, 2013).A reduction of the resolution dependence will improve the efficiency and accuracy of research and operational activities that require data generated by GCMs.
The effect that the small unresolved subgrid scales have on the large resolved scales is typically parameterised by defining a form of eddy viscosity.In most subgrid models, including the most widely celebrated and adopted ones (Smagorinsky, 1963;Gent and McWilliams, 1990), physical arguments are used to justify the form of an eddy viscosity, which is then tuned to achieve numerical stability and realistic results Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.(Griffies et al., 2005).In practice, however, there is a significant range of small scales that are excessively damped (dissipation range) due to the application of heuristic subgrid turbulence models, which in turn also affect the large scales (Frederiksen et al., 2003).Ideally one would prefer not to have an artificial dissipation range, and develop a subgrid model that renders all of the scales of motion accurate.The subgrid scales also contribute to predictability limitations by injecting noise into the system.It has been shown that weather and climate models with deterministic subgrid models have insufficient ensemble spread, a situation that is improved with the injection of stochastic backscatter (Leith, 1990;Frederiksen and Davies, 1997;O'Kane and Frederiksen, 2008;Shutts, 2005;Grooms et al., 2015;Franzke et al., 2007Franzke et al., , 2015;;Shutts, 2015).
As in general it is only possible to parameterise the statistical effects of the subgrid eddies (McComb et al., 2001), statistical dynamical closure theory is the natural formulation for developing self-consistent subgrid models.In this approach one attempts to determine the statistical effect that the unresolved scales of motion have on the resolved eddies.The foundation studies in this area were the direct interaction approximation (DIA) closure and its variants for homogeneous turbulence (Kraichnan, 1959;McComb, 1974;Herring, 1965), and the quasi-diagonal DIA (QDIA) closure (Frederiksen, 1999(Frederiksen, , 2012a;;O'Kane andFrederiksen, 2004, 2008) for inhomogeneous turbulence.The general QDIA closure theory accounts for cross-correlations between field variables (e.g.fields at different vertical levels; or velocity components) and between physical space fields, but has the remarkable property that the eddy damping and stochastic backscatter terms are diagonal in spectral space.The QDIA subgrid closure terms were calculated for typical barotropic atmospheric flows in O' Kane and Frederiksen (2008).Broadening the applicability of the QDIA closure, a stochastic subgrid modelling approach was developed to determine the eddy viscosities from the statistics of high-resolution benchmark simulations (Frederiksen and Kepert, 2006), which is the approach adopted here.
We use the method of Frederiksen and Kepert (2006) to develop stochastic subgrid models for global atmospheric and oceanic flows such that practically all of the resolved scales of motion can be trusted.In contrast to the vast majority of subgrid modelling studies, the approach adopted here makes no heuristic assumptions, with the subgrid model coefficients calculated self-consistently from the statistics of high-resolution benchmark simulations.This approach has been successfully applied to quasi-geostrophic (QG) atmospheric and oceanic simulations with horizontal and vertical shears (Zidikheri andFrederiksen, 2009, 2010a, b), threedimensional wall bounded turbulence (Kitsios et al., 2015), and global primitive equation simulations of the atmosphere (Frederiksen et al., 2015).Subgrid models developed from far simpler barotropic QG models (Frederiksen and Davies, 1997), have previously been shown to improve the simulated dynamics in GCMs (Frederiksen et al., 2003).Here we adopt more complex baroclinic QG benchmark simulations of the atmosphere and ocean, which capture the essential dynamics of barotropic (horizontal shear) and baroclinic (vertical shear) instability.
Historically subgrid models for the atmosphere and ocean have been developed in isolation, with the derivation of the functional forms of the subgrid models often motivated by very different physical phenomena.Here we provide evidence that the effects of subgrid turbulence in the atmosphere and ocean actually have much in common.When non-dimensionalised appropriately, subgrid coefficients calculated from atmospheric (Kitsios et al., 2012) and from oceanic (Kitsios et al., 2013) simulations, show remarkably good agreement within the enstrophy cascading inertial range.The justification of this approach stems from the phenomenological view of turbulence in the atmosphere and ocean.In both flows the Rossby radius (r R ) is the dominant scale at which baroclinic instability injects energy (velocity variance) and enstrophy (vorticity variance) into the system (Salmon, 1998), where the non-dimensional Rossby wavenumber is k R ≡ a/r R , with a = 6371 km the radius of the Earth.In the phenomenological view of QG turbulence, enstrophy is transferred at a constant rate from wavenumber k R to larger wavenumbers (smaller eddies), whilst energy is transferred from wavenumber k R back up to the large-scale (low wavenumber) energy-containing eddies of wavenumbers less than or equal to k E (Kraichnan, 1976;Salmon, 1998).The wavenumbers, k R and k E , divide the scales into three important wavenumber (n) regimes: the non-self-similar energy-containing range (n ≤ k E ), the selfsimilar inverse energy cascade (k E < n ≤ k R ), and the selfsimilar forward enstrophy cascade (k R < n).In the ocean k E k R with all three regimes present.In the atmosphere, however, k E ≈ k R , which means that the inverse energy cascade is either very short or non-existent, due to the largescale forcing.Both wavenumbers, k R and k E , are important for the scaling of the subgrid coefficients.
Here we present a first systematic comparison of subgrid models of QG turbulence in the atmosphere and ocean, and develop simple unifying scaling laws that represent both fluid media within their enstrophy cascading inertial ranges.A large set of simulations is analysed, which covers a broad range of flow parameters, including an order of magnitude change in the Rossby radius of deformation and the energycontaining scale.By focussing on the enstrophy cascading inertial range in both media, the large number of simulations and wide parameter range has enabled the establishment of robust scaling laws.In Sect. 2 we present the numerical details of the benchmark simulations used to generate the atmospheric and oceanic flows, with these flows characterised in Sect.3. The process by which subgrid models are calculated from the reference benchmark simulations is presented in Sect.4, with the resulting subgrid coefficients illustrated in Sect. 5.The coefficients calculated from the atmospheric and oceanic simulations are then characterised into a set of unifying scaling laws representing both fluids in Sect.6.These unifying scaling laws govern how the subgrid coefficients change with resolution and flow strength, thus removing the need to generate the coefficients from benchmark simulations in the future.The scaling laws presented here are particularly simple, and are suggestive of robust fundamental properties of QG turbulence.In Sect.7, large eddy simulations adopting these scaling laws are shown to reproduce the statistics of the benchmark simulations across all scales, with drastic improvements in computational efficiency.

Numerical details of the benchmark simulations
The atmospheric and oceanic flows are generated by solving the two-level QG potential vorticity equation (QGPVE).The numerical integration of the QGPVE is a computationally efficient means of simulating geophysical flows.It captures the essential dynamics of baroclinic and barotropic instabilities, and the interaction of coherent structures with inhomogeneous Rossby wave turbulence (Frederiksen, 1998).
In the present study the vorticity is represented on two discrete vertical levels with j = 1 representing the upper level and j = 2 the lower level.In the atmospheric simulations the upper level is at 250 hPa (≈ 10 km), and the lower level at 750 hPa (≈ 2.5 km).For the oceanic simulations the upper level is at an approximate depth of 200 m, and the lower level at 600 m.The system is non-dimensionalised by using the radius of the Earth (a = 6371 km) as a length scale, and the inverse of the Earth's angular velocity ( = 7.292 × 10 −5 s −1 ) as a timescale.By default all variables are assumed to be non-dimensional unless units are specified.
The two-level QG, equations of motion in physical space are The field variables are functions of time (t), longitude (λ), and µ = sin(φ), where φ is the latitude.The vorticity at level j is ζ j , and ψ j is the stream function.The reduced potential vorticity q j ≡ ζ j + (−1) j F L ψ 1 − ψ 2 , where F L is the layer coupling coefficient, which is inversely proportional to the temperature difference between the two levels, and is related to the Rossby radius of deformation by r R = 1/ √ 2F L .In Eq. (1), the coefficient B represents the beta effect, and J (ψ j , q j ) is the Jacobian.Using standard fluid mechanical nomenclature, D j 0 is the bare dissipation operator representing the unresolved eddy-eddy (or inter-eddy) interactions in the benchmark simulation (McComb, 1990).The constant α j parameterises the drag by dampening the large scales of motion.Simulations are nudged toward a climate q j by the constant relaxation parameter κ j .
In our study we solve Eq. ( 1) by spectrally discretising the field variables in spherical harmonics (Frederiksen, 1998).This spectral discretisation allows for a clear separation of the resolved and subgrid scales of motion for the development of the subgrid parameterisations.The system solves for the spectral coefficients of the potential vorticity defined as for zonal (longitudinal) wavenumber, m, total wavenumber, n, with latitudinal (meridional) wavenumber n−m.The spectral coefficients of the vorticity are ζ j mn = −n(n + 1)ψ j mn , where ψ j mn are the spectral stream-function coefficients.The evolution of q j mn is governed by where q j −mn is the complex conjugate of q j mn , and K mpr nqs are the interaction coefficients defined in Frederiksen and Kepert (2006).No topography is represented in the present simulations.The summations immediately after the equals sign in Eq. ( 3) are over the triangular wavenumber set with T the benchmark truncation wavenumber, which is related to the angular grid spacing in degrees ( ) by T = 120/ .The highest resolution atmospheric and oceanic simulations run for the present study have maximum truncation wavenumbers of T = 504 and T = 1008, respectively.The Rossby wave frequency is ω mn = −Bm/[n(n + 1)], where B = 2 under the chosen non-dimensionalisation.In the atmospheric simulations F L = 2.5 × 10 −12 m −2 , corresponding to a Rossby radius of deformation of r R ≡ 1/ √ 2F L = 447 km, and a non-dimensional Rossby wavenumber of k R ≡ r R /a ≈ 14.In the oceanic cases In Eq. ( 3), α j (n) is the drag applied at level j .In the atmospheric simulations α j (n) = α j max for n ≤ 15, and zero otherwise, with α 1 max = 2.3 × 10 −6 s −1 and α 2 max = 5.8 × 10 −7 s −1 .For the simulations of the ocean α , where erf is the error function, and n c = 50 is the point at which α j (n c ) = α j max /2.This functional form allows us to control the location of the energycontaining wavenumber.We undertake additional oceanic simulations with alternate values of n c to produce a series of flows with different background states and with differing wavenumber ranges of the energy-containing (non-selfsimilar) scales (k E ).

V. Kitsios et al.: Theoretical comparison of subgrid turbulence
All simulations are driven toward a mean state q j mn that is purely zonal ( q j mn are zero unless m = 0).They are driven toward this state by a relaxation parameter of κ j n = 10 −6 s −1 for m = 0 and n ≤ 15, and zero otherwise.For the simulations of the atmosphere q j mn corresponds to a large-scale westerly jets centred at 45 • S and 45 • N, representing largescale jets in the Northern and Southern hemispheres.In the oceanic simulations q j mn corresponds to a large-scale westerly current centred at 60 • S of the Southern Hemisphere, broadly representative of the mean Antarctic Circumpolar Current.
By definition the bare dissipation, D j l 0 (m, n), represents the unresolved eddy-eddy interactions in the benchmark simulation.It is written in general anisotropic matrix form (dependent on zonal, m, and total, n, wavenumbers) but in our simulations it has the isotropic form (dependent only on and δ lj is the Kronecker delta function, which ensures the off-diagonal elements of ν j l 0 (n) are zero.Here ν jj 0 (T ) is the value of the diagonal elements at truncation and the power ρ j 0 determines the steepness of ν jj 0 (n).This means that the corresponding bare viscosity and bare dissipation matrices are diagonal and isotropic.Note in Eq. ( 5), the wavenumber ratio n/T is raised to the power of ρ j 0 −2 to be consistent with the definition of the subgrid eddy viscosities throughout the document.The slope and magnitude of ν 0 is determined by the scaling laws presented in the paper.An initial study was first undertaken determine the scaling laws with an estimate of ν 0 .The study was then repeated with ν 0 defined by the scaling laws themselves.There was a negligible difference between the new subgrid coefficients and those obtained in the initial study.

Characterisation of the benchmark flows
In the benchmark atmospheric simulations, the Rossby radius of deformation r R = 447 km, with an associated wavenumber of k R = 14.This means 14 eddies of this size could fit side by side along one line of latitude.The climate state contains large-scale westerly winds in the mid-latitudes of the Northern and Southern hemispheres (Kitsios et al., 2012); see Fig. 1c.Large-scale eddies are produced in both hemispheres as illustrated by the instantaneous eddy stream function and wind field in Fig. 1a.
In the initial benchmark oceanic simulation the Rossby radius is 45 km corresponding to a wavenumber of k R = 142.The Rossby radius is an order of magnitude smaller in the ocean compared with the atmosphere.This renders oceanic simulations computationally more expensive, as a finer grid is required to explicitly resolve baroclinic instability.The climate state is illustrated in Fig. 1d, and is broadly representative of the Antarctic Circumpolar Current (Kitsios et al., 2013).Figure 1b illustrates that the oceanic flow has eddies in the mid-latitudes of the Southern Hemisphere that are smaller in size than those in the atmospheric case, and is consistent with the former having a smaller Rossby radius.
The strength of the flow field on each level is quantified by the potential enstrophy flux, and is required for scaling the magnitude of the eventual subgrid coefficients.The enstrophy flux, η j k (n), is the rate of potential enstrophy transfer from level k into level j at total wavenumber n.It is defined as is the enstrophy transfer.The latter is calculated by postmultiplying the non-linear term of the equations of motion in Eq. ( 3) by q k −mn , and then summing over zonal wavenumber m.The potential enstrophy flux for the atmospheric and oceanic simulations are illustrated in Fig. 2a and b.
The wavenumber extent of the large energy-containing scales is required for scaling the spectral slope of the subgrid coefficients.Within the inertial ranges the external forcing and dissipation are negligible, and the transfer of energy is dominated by non-linear triadic interactions (Salmon, 1998).With no additional damping or excitation within the self-similar wavenumber regimes, we find that the energy transferred into the barotropic mode is in balance with that transferred out of the baroclinic mode.We define k E to be a wavenumber indicative of the non-self-similar energycontaining scales.It is quantified by the smallest wavenumber at which the energy transferred into the barotropic mode is in balance with the energy transferred out of the baroclinic mode (Kitsios et al., 2013).The kinetic energy transfers in level space are given by T j k (n) = N j k (n)/[n(n + 1)].The barotropic/baroclinic kinetic energy transfers are given by T where and the superscript T denotes the transpose operation.The index 1 refers to the barotropic mode, and 2 the baroclinic mode.For example T 12 B (n) refers to the kinetic energy transferred from the baroclinic mode into the barotropic mode.The energy transferred into the barotropic mode is T BT (n) = T 11 B (n)+T 12 B (n), and likewise the energy transferred into the baroclinic mode is To be in balance, T BT (n) must be equal to −T BC (n).For the atmospheric flow we find k E ≈ 11, and for the oceanic flow k E = 70, as illustrated in Fig. 2c and d, respectively.

Stochastic subgrid modelling approach
Using a series of the above-discussed simulations, we study the inter-eddy interactions by removing vortices smaller than a certain cut-off size, or equivalently larger than a specified truncation wavenumber (T R ).The subgrid tendency is the component of the rate of change of the resolved large-scale vortices due to their interactions with the unresolved small scale vortices.The subgrid parameterisation problem in its most basic form is the representation of the subgrid tendency in terms of the resolved field.Here we use the stochastic subgrid modelling approach of Frederiksen and Kepert (2006) to determine such a representation for the subgrid processes.This approach is outlined below.The resolution of a LES is lower than the associated benchmark simulation, and confined to the resolved scale wavenumber set where T R is the LES truncation wavenumber such that T R < T .The subgrid wavenumber set is defined as S = T − R. We define the resolved potential vorticity field at a given wavenumber pair (m, n) by the two-element column vector q = (q 1 mn , q 2 mn ) T .In this vector notation, q t (t) = q R t (t) + q S t (t), where q t is the tendency (time derivative) of q.The tendency of the resolved scales is q R t , where all triadic interactions involve wavenumbers less than T R .The remaining subgrid tendency q S t has at least one wavenumber greater than T R , which is involved in the triadic interactions.One can further decompose q S t such that q S t (t) = f + q S t (t), where q S t is the fluctuating component representing the eddy-eddy interactions, and f ≡ q S t is the ensemble-averaged subgrid tendency representing the sum of the eddy-mean-field and mean-field-mean-field interactions.
The QDIA closure provides the theoretical justification for modelling the subgrid tendency for a particular wavenumber pair as a function of the resolved fields at only that same wavenumber pair (Frederiksen, 2012a).We can then model the fluctuating subgrid tendency at each wavenumber pair, q S t , by the stochastic equation where D d is the subgrid drain dissipation matrix, q is the fluctuating component of q, and f is a random forcing vector.As the present simulations have two vertical levels, D d is a time-independent 2 × 2 matrix, and f is a time-dependent two-element column vector.An estimate of D d is then found through the generalisation of the Gauss theorem (Frederiksen and Kepert, 2006).Both sides of Eq. ( 10) are post-multiplied by q † (t 0 ), integrated over the turbulent decorrelation period τ , ensemble averaged to minimise the contribution from f , and then rearranged to produce where † denotes the Hermitian conjugate for vectors and matrices.The angled brackets denote ensemble averaging, with each ensemble member determined by shifting t 0 forward by one time step.The decorrelation time τ , is chosen sufficiently large to capture the memory effects of the turbulence (Kitsios et al., 2012).The model for f is then determined by calculating the matrix . Postmultiplying both sides of Eq. ( 10) by q † (t), and adding the conjugate transpose of Eq. ( 10) pre-multiplied by q(t) yields the Lyapunov equation Given that D d has been determined, F b can now be calculated.There is a balancing act between the linear (D d ) and stochastic (F b ) components of the subgrid model.As D d is dependent upon τ , it is τ that defines this balance.For the implementation of parameterisation, it is sufficient to assume that f can be represented as the white noise process , with an eigenvalue decomposition of F b used to produce a stochastic model for f , as detailed in Zidikheri and Frederiksen (2009).
Backscatter is the physical process by which kinetic energy is transferred from small to large scales.The subgrid model in Eq. ( 10) represents this process in its fundamental stochastic form.One can also, however, represent the subgrid interactions using the simplified deterministic form q S t (t) = −D net q(t), where D net is the net dissipation representing the net effect of the drain and backscatter (Frederiksen and Kepert, 2006).The backscatter and net linear operators are defined by

and
(13) respectively (Frederiksen and Kepert, 2006).In the present document the subgrid coefficients are presented in eddy viscosity form, where the drain, backscatter, and net eddy viscosities are related to their respective dissipations by and where n(n + 1) is the discrete form of the Laplacian.

Structure of the eddy viscosities
For the atmosphere the subgrid model coefficients are presented at a truncation of T R = 126, capturing vortices down to a radius of 50 km in the mid-latitudes.These eddies are significantly smaller than the Rossby radius (447 km), which means the energy injected into the system via baroclinic instability is explicitly resolved.In Fig. 3a the upper diagonal element of the drain eddy viscosity is divided by the kinematic viscosity of air (10 −5 m 2 s −1 ), and represented by the height of the contour surface.The coloured surface depicts the kinetic energy of the fluctuating scales at the upper level.In this figure the eddy viscosity is 10 10 times greater than the molecular viscosity, indicating that the intereddy interactions are far more important than the intermolecular ones.The drain also increases strongly with the total wavenumber (n), has only a weak dependence upon the zonal wavenumber (m) at a given n, and is hence approximately isotropic.The kinetic energy is also largely isotropic, concentrated at the largest scales (lowest wavenumbers), and decreases rapidly as the structures get smaller (wavenumbers get larger).The form and magnitude of the lower diagonal element of the drain eddy viscosity matrix are very similar to those of the upper diagonal element, with the off-diagonal elements negligible in comparison.Since the drain eddy viscosity matrix is essentially diagonal, the positive coefficients illustrated in Fig. 3a indicate that energy (and enstrophy) is being sent from the resolved to the subgrid eddies.The backscatter has a similar form to the drain, but is negative and approximately half the magnitude.
We now consider the drain eddy viscosity in the ocean at the same resolution of T R = 126, again capturing vortices of radius 50 km.Here, the energy injection via baroclinic instability is not explicitly resolved as the Rossby radius is 45 km.The upper diagonal-drain eddy viscosity component is divided by the kinematic viscosity of sea water (10 −6 m 2 s −1 ) and plotted in Fig. 3b.It again illustrates that the influence of the inter-eddy interactions is 10 10 times greater than the inter-molecular ones.The eddy viscosity is strongly dependent upon both zonal (m) and total (n) wavenumbers, and is hence anisotropic.For certain low wavenumbers (large scales) the drain is negative, which is required to further deterministically excite the flow as the injection of energy via barotropic and baroclinic instabilities is not explicitly resolved.The coloured surface depicts the upper level kinetic energy, illustrating that it is also highly anisotropic and distributed across all scales.The lower diagonal matrix element has similar properties to the upper diagonal.The off-diagonal elements are proportionally larger in this case, indicating that the removal of the small scales modifies the interactions between the vertical levels -refer to Kitsios et al. (2013) for illustrations of the off-diagonal elements.Jansen and Held (2014) developed heuristic general purpose oceanic subgrid models for this regime that also have negative viscosity.For oceanic simulations at the higher resolution of T R = 252, in which baroclinic instability is explicitly resolved, the eddy viscosities have similar properties to the atmospheric case, with matrices diagonally dominant and largely isotropic (Kitsios et al., 2013).
The self-similarity of the eddy viscosities is most clearly illustrated by the isotropised (averaged over zonal wavenumber m) profiles.For various truncations levels (T R ), the upper diagonal element of the isotropised drain and backscatter eddy viscosities is illustrated in Fig. 3c for the atmospheric flow, and in Fig. 3d for the ocean.We also show the net eddy viscosity, given by the sum of the drain and backscatter.As the resolution increases the magnitude of all of the eddy viscosities decrease.This means that as more eddies are being explicitly resolved, the enstrophy (and energy) is being transferred to fewer subgrid eddies.For cases that resolve baroclinic instability, the subgrid parameterisation represents the energy flow to the resolved scales as being completely stochastic with only the backscatter eddy viscosity negative.The positive values of the net eddy viscosity indicate that the net effect of the drain and backscatter processes is such that energy is sent out of the system.When baroclinic instability is not resolved the energy flow to the resolved scales is modelled as having a deterministic com-ponent with the drain and net eddy viscosities negative for certain low wavenumbers.The eddy viscosity coefficients with significant magnitude are concentrated within the last 70 wavenumbers for the ocean and the last 11 wavenumbers for the atmosphere.These wavenumber ranges coincide with k E -the wavenumber to which the large non-self-similar energy-containing scales extend.

Unifying scaling laws
We have calculated the subgrid parameterisation coefficients (eddy viscosities) for the atmosphere and ocean at various resolutions (T R ).We now develop scaling laws representing how these eddy viscosities change with resolution and flow strength, for truncations made within the enstrophy cascading inertial range (k R < T R ).For the diagonal element of the drain eddy viscosity associated with level j , the maximum magnitude (ν jj d (T R )) and spectral slope (ρ j d ) are quantified by least squares fitting the isotropised eddy viscosity profiles (ν There is an analogous expression for the isotropised backscatter eddy viscosity (ν jj b (n)).The scaling laws govern how the magnitudes and slopes change with truncation wavenumber and flow strength.Oceanic benchmark simulations were also undertaken, with the Rossby wavenumber (k R ) varying from 142 to 284, and the energy-containing wavenumber (k E ) varying from 40 to 70.This coupled with the atmospheric results (k R = 14, k E = 11), means that we have results spanning almost an order of magnitude in both the Rossby and energy-containing wavenumbers.
First, we present the power exponents of the drain eddy viscosities (ρ j d ), which represent how steeply the drain of enstrophy out of the system increases with resolved wavenumber (or equivalently as the size of the resolved eddies decrease).It is the extent of the energy-containing scales (k E ) that defines how far non-linear interactions can span in wavenumber space (Kraichnan, 1976), which effectively sets the size of the largest eddy that can interact with the subgrid scales.This wavenumber distance is inversely proportional to the power exponents, and is represented by the span of wavenumbers over which the eddy viscosity profiles are non-zero in Fig. 3c and d.In Fig. 4a, we therefore plot the drain power exponent against the truncation wavenumber (T R ) non-dimensionalised by k E .A strong relationship exists for all of the atmospheric and oceanic flows, with the drain exponent increasing with T R .The spectral slope has to increase with resolution to ensure that the range of significant subgrid interactions (quantified by the eddy viscosity) is confined to the last k E wavenumbers before truncation.The scaling law for ρ j d is determined by the illustrated regres-  sion line.A similar relationship is observed for the power exponents of the backscatter eddy viscosities (ρ j b ) in Fig. 4b, with the dashed line illustrating the scaling law for the drain to serve as a direct comparison.Note the backscatter power exponents are larger and also increase with resolution more quickly than the drain exponents.To put these results into context, a power exponent of 2 represents a Laplacian dissipation, or equivalently an eddy viscosity that does not depend on wavenumber.
Scaling laws for the maximum values are again nondimensionalised using the energy-containing wavenumber, and additionally a timescale based on the potential enstrophy flux (Leith, 1971).The potential enstrophy flux is the rate at which potential enstrophy is transferred from one wavenumber to the next (Salmon, 1998).We calculate the flux and find that for both flow cases it is constant for eddies smaller than the energy-containing scale, as illustrated in Fig. 2c and d.The constant flux value at level j is denoted by η jj I .To span all cases of different k R and k E , we find that the eddy viscosities need also to be scaled by √ k R /k E .With this normalisation, the magnitude of the drain and backscatter are plotted in Fig. 4c and d, respectively.The magnitude of all eddy vis-Table 1. Equivalent powers of the Laplacian for the subgrid net eddy viscosity of atmospheric and oceanic simulations at various angular grid spacings ( ).The equivalent truncation wavenumber is T R = 120/ .The energy-containing scale for the atmospheric and oceanic simulations are k E = 11 and k E = 70, respectively.The drain profiles (ν d (n)) are calculated from T R /k E using Eq. ( 15) and the scaling laws in Fig. 4, and likewise for the backscatter (ν b (n)), with the net ν net (n) = ν d (n) + ν b (n).The spectral slope of ν net (n) is determined, divided by 2, and rounded to the nearest integer to approximate the effective power of the Laplacian operator.cosities is inversely proportional to T R , which means that if the resolution doubles the eddy viscosity halves.These scaling laws allow us to determine the drain and backscatter terms at the desired resolution (T R ), given that we have estimates of the Rossby wavenumber, energycontaining wavenumber, and enstrophy fluxes.These terms can then be used to model the subgrid interactions in simulations of the climate.Whilst the scaling laws were developed from baroclinic QG simulations, they agree with the subgrid coefficients determined from the truncation of barotropic (Frederiksen and Kepert, 2006) and more complex atmospheric multi-level primitive equation simulations (Frederiksen et al., 2015).This indicates that the scaling laws can be applied more broadly.As previously mentioned, subgrid models developed from simpler barotropic QG models (Frederiksen and Davies, 1997), have been shown to improve the simulated dynamics in GCMs (Frederiksen et al., 2003).As most GCMs run with deterministic subgrid models, in Table 1 we list the effective spectral slope of the net eddy viscosity at various resolutions (T R ) for typical atmospheric (k E = 11) and oceanic (k E = 70) flows.For a given T R /k E , the drain profile (ν d (n)) is calculated using Eq. ( 15) and the scaling laws in Fig. 4, and likewise for the backscatter (ν b (n)), with the net eddy viscosity given by ν The spectral slope of ν net (n) is then calculated, divided by 2, and rounded to the nearest integer to approximate the effective power of the Laplacian.For a given resolution, atmospheric simulations are far more scale selective than oceanic ones, because the extent of the energycontaining scales (k E ) is significantly less in the atmosphere than in the ocean.

Large eddy simulation
We now determine if LES with subgrid models defined by the eddy viscosities presented above, can replicate the statistics of the higher-resolution benchmark simulations.The equation governing the LES is equivalent to that of the benchmark simulation in Eq. ( 3), with the addition of the term added to the right-hand side, and solved over the wavenumber set R instead of T. A stochastic model for f is built from an eigenvalue decomposition of F b (Zidikheri and Frederiksen, 2009).In the deterministic form, the stochastic force f is removed and ν d is replaced with ν net .In the isotropic cases, the matrices ν d , ν net , and ν b are averaged over the zonal wavenumbers m so that they are only functions of the total wavenumbers n.
We compare the benchmark simulation results to LES comprising of both stochastic and deterministic subgrid models, with the model coefficients in their original anisotropic form (as in Fig. 3a), in their isotropised form (as in Fig. 3c), and also defined by the associated scaling laws.Comparisons are made across all scales of motion on the basis of the time-averaged zonal (m) wavenumber-summed kinetic energy spectra.The upper level spectra of the benchmark simulations (black dashed line) are compared to that of the LES (red solid line) labelled by the associated subgrid parameterisation variant in Fig. 5.The top pair of spectra represent the true energy level, with the other pairs of spectra shifted down for clarity.Findings pertaining to the upper level are consistent with those for the lower level.
The atmospheric benchmark simulation of maximum wavenumber T = 504 is compared to LES with T R = 63 in Fig. 5a.The stochastic and deterministic variants with anisotropic-, isotropic-, and scaling-law-defined coefficients all reproduce the kinetic energy of the benchmark simulation across all scales of motion.As the resolution is reduced in both horizontal directions, the number of degrees of freedom is reduced by (T 2 − T 2 R )/T 2 = (512 2 − 63 2 )/512 2 = 98 %.This reduced resolution also allows us to decrease the time step proportionally, which means the computational cost of the simulation is reduced by a factor T 3 /T 3 R = 512 3 /63 3 = 537.The oceanic benchmark simulation of T = 504 is com- pared to LES with T R = 252 in Fig. 5b.Again all LES variants replicate the statistics of the benchmark simulation.This represents a 75 % reduction in the degrees of freedom, a decrease in computational cost by a factor of 67.In summary for both the atmosphere and ocean, the idealised scaling law form of the eddy viscosities is an excellent representation of the subgrid interactions within the enstrophy cascade.We have also developed scaling laws applicable to the ocean within the inverse energy cascade (k E < n < k R ), as discussed in Kitsios et al. (2013).

Conclusions
A general stochastic modelling approach (Frederiksen and Kepert, 2006) has been used to determine eddy viscosity matrices that parameterise the interactions between fields at different vertical levels and horizontal scales in the atmosphere and ocean.Additionally when truncations are made within the enstrophy cascading inertial range the subgrid parameterisation coefficients are represented by a set of unifying scaling laws.The laws govern how the form and magnitude of both the atmospheric and oceanic eddy viscosities change with flow strength and grid resolution.We have demonstrated that simulations adopting these scaling laws produce resolution-independent statistics across all scales of motion.This means no additional resolution need be wasted in order to account for the presence of an artificial dissipation range, which drastically improves the computational efficiency of the simulations.
The scaling laws developed here can be implemented directly into spectral simulations, and are expected to improve the efficiency and accuracy of numerical weather and climate simulations (Frederiksen et al., 2003(Frederiksen et al., , 2015)).There are also two possible approaches to implement these scaling laws into grid point codes.The simplest approach is to apply the subgrid model directly in grid-point space via a Laplacian operator of the appropriate power, as outlined in Table 1.More generally it is also possible to employ grid to spectral transforms, where the subgrid model is calculated in spectral space, and then applied in physical space.
Finally, the stochastic modelling approach adopted here is not confined to fluid mechanics but can also be used to represent non-linear interactions in any classical multi-scale dynamical system (Frederiksen, 2012b).

Figure 1 .Figure 2 .
Figure 1.Instantaneous fields and climate states of the benchmark simulations.Contours of instantaneous eddy (non-zonal) stream function, and vectors of instantaneous velocity (wind/current) on the upper level of the atmosphere (Northern and Southern hemisphere), and (b) ocean (Southern Hemisphere).Climate state illustrated by the time-averaged (c) atmospheric winds and (d) oceanic currents.

Figure 3 .
Figure 3. Real component of the upper diagonal subgrid eddy viscosities.Anisotropic drain eddy viscosity at T R = 126 for the (a) atmosphere, divided by the kinematic viscosity of air (10 −5 m 2 s −1 ); and (b) ocean, divided by the kinematic viscosity of water (10 −6 m 2 s −1 ).Coloured surfaces depict kinetic energy of the fluctuations at the upper level, and the black lines are the isotropised (m-averaged) drain coefficients.Isotropic drain, backscatter, and net eddy viscosities labelled by T R for the (c) atmosphere and (d) ocean; with Rossby wavenumber (k R ) and energy-containing wavenumber (k E ) labelled on the horizontal axes.

Figure 4 .
Figure 4. Scaling of the isotropic eddy viscosities.Slope of the (a) drain eddy viscosity ρ j d , and (b) backscatter eddy viscosity ρ j b .Magnitude of the (c) drain eddy viscosity ν jj d (T R ), and (d) negative backscatter eddy viscosity −ν jj b (T R ).The dashed lines in panels (b) and (d) represent the drain scaling laws to serve as a direct comparison to the backscatter scaling laws represented by the solid lines.The symbols correspond to the various cases as follows: red diamond, k R = 142, k E = 70; blue circle, k R = 284, k E = 70; green square, k R = 142, k E ∈ (40, 50, 60); magenta upward pointing triangle, k R ∈ (201, 246), k E = 70; orange downward pointing triangle, k R = 14, k E = 11 (atmosphere).Filled symbols represent j = 1 and hollow symbols j = 2.

Figure 5 .
Figure 5. Scale by scale comparison of the benchmark simulation (dashed line) to the LES variants (red solid line).Kinetic energy spectra at the upper level of the (a) atmosphere and (b) ocean.The top pair of spectra exhibit the true energy, with subsequent pairs shifted down for clarity.Spectra are labelled with the associated subgrid parameterisation of anisotropic stochastic (AS), anisotropic deterministic (AD), isotropic stochastic (IS), isotropic deterministic (ID), scaling law stochastic (LS), or scaling law deterministic (LD).The truncation (T R ), Rossby (k R ) and energy-containing (k E ) wavenumbers are labelled on the horizontal axis.