Emergent eddy saturation from an energy constrained eddy parameterisation

The large-scale features of the global ocean circulation and the sensitivity of these features with respect to forcing changes are critically dependent upon the inﬂuence of the mesoscale eddy ﬁeld. One such fea-ture, observed in numerical simulations whereby the mesoscale eddy ﬁeld is at least partially resolved, is the phenomenon of eddy saturation, where the time-mean circumpolar transport of the Antarctic Circumpolar Current displays relative insensitivity to wind forcing changes. Coarse-resolution models employing the Gent–McWilliams parameterisation with a constant Gent–McWilliams eddy transfer coe ﬃ cient seem unable to reproduce this phenomenon. In this article, an idealised model for a wind-forced, zonally sym- metric ﬂow in a channel is used to investigate the sensitivity of the circumpolar transport to changes in wind forcing under different eddy closures. It is shown that, when coupled to a simple parameterised eddy energy budget, the Gent–McWilliams eddy transfer coe ﬃ cient of the form described in Marshall et al. (2012) [ A framework for parameterizing eddy potential vorticity ﬂuxes , J. Phys. Oceanogr., vol. 42, 539–557], which includes a linear eddy energy dependence, produces eddy saturation as an emergent property.


Introduction
Studies of the response of the large-scale ocean circulation to changing forcing scenarios in numerical ocean models require long time integrations that are prohibitively expensive even at mesoscale eddy permitting resolutions. Since this is expected to remain the case for the foreseeable future, an ongoing challenge in numerical ocean modelling is the representation of the unresolved mesoscale eddy field in coarse resolution models. A particularly successful scheme that is employed is the Gent-McWilliams (GM) parameterisation ( Gent and McWilliams, 1990;Gent et al., 1995 ), which parameterises mesoscale eddies via the introduction of a non-divergent eddy transport velocity. The eddy transport velocity can be interpreted as arising from the difference between the Eulerian average of the velocity at fixed height and the thicknessweighted average of the velocity at fixed density ( McDougall and McIntosh, 2001 ), and modifies the advective transport of tracer quantities. By definition, the non-divergent eddy transport velocity conserves all moments of the advected quantities, and is thereby * Corresponding author.
adiabatic. The property of adiabatic stirring is particularly attractive, being shown to remove spurious heating and cooling in the deep ocean, such as that associated with the Deacon cell in the Southern Ocean ( Danabasoglu et al., 1994 ).
To this point, studying the modelled oceanic response to changing atmospheric forcing in conjunction with the GM parameterisation is of particular importance for emergent climatologies under different forcing scenarios. Two important large-scale Southern Ocean phenomena are of particular interest in this regard. The first is "eddy saturation", originally discussed in Straub (1993) from an argument based on critical stability, and here to be understood as the relative insensitivity of the time-mean circumpolar transport with respect to wind forcing changes. The other is "eddy compensation", here to be understood as the reduced sensitivity of the residual meridional overturning circulation with wind forcing changes (e.g., Meredith et al., 2012;Viebahn and Eden, 2012;Munday et al., 2013 ), which has consequences for the meridional transport of important tracers such as heat, salt and carbon. This article focuses on eddy saturation.
As argued by Straub (1993) , if fluid interaction with topography is the main sink for momentum input by wind stress, and consequently the zonal abyssal flow is weak, then thermal wind shear is the dominant contribution to circumpolar transport; Peña-Molino et al. (2014) suggest that thermal wind shear accounts for at least 75% of the net circumpolar transport in the Southern Ocean. Thus circumpolar transport is intimately linked to isopycnal slope, with the slope steepness limited by baroclinic instability. Eddy saturation arises through a balance between steepening of isopycnals by wind stress, and flattening of isopycnals by the presence of the mesoscale eddy field. While the question of whether the ocean is in an eddy saturated state remains unconstrained by current observations, the reduction in circumpolar transport sensitivity with varying wind stress has been observed in a variety of numerical models that at least partially resolve a mesoscale eddy field (e.g., Hallberg and Gnanadesikan, 2006 ;Hogg and Blundell, 2006;Hogg et al., 2008; . In Munday et al. (2013) , an eddy permitting one-sixth degree model of a 20 °wide ocean sector was integrated with varying wind forcings. This eddy permitting model, employing a very small value of the GM eddy transfer coefficient, showed near complete eddy saturation. By contrast, in lower resolution half degree and two degree variants of the same model, where larger values of the GM eddy transfer coefficient were utilised, the resulting time-mean circumpolar transport displayed significant sensitivity with respect to the wind forcing. Hogg and Munday (2014) found that although the value of the time-mean circumpolar transport was affected by the domain geometry, the relative insensitivity with changing wind stress at eddy permitting resolution was robust.
Thus it has been found that the GM scheme with a spatially and temporally constant GM eddy transfer coefficient is unable to represent eddy saturation (see also Farneti et al., 2015 ). With increased wind forcing, a more vigourous eddy field is to be expected. Since the GM eddy transfer coefficient in some sense specifies the intensity and efficiency of the parameterised eddy field, it is expected that a positive correlation between the strength of wind forcing and the magnitude of the GM coefficient eddy transfer is minimally required for emergent eddy saturation. Various proposals already exist with a non-constant GM eddy transfer coefficient. In Visbeck et al. (1997) , using linear stability arguments, a GM eddy transfer coefficient is proposed which depends upon the stratification profile, as well as a mixing length. In Ferreira et al. (2005) the eddy-mean-flow interaction in a global ocean model is determined via an optimisation procedure, yielding diagnosed values for the GM eddy transfer coefficient. Their optimisation is used to infer a GM eddy transfer coefficient which depends on the vertical stratification, and has subsequently been incorporated into a number of ocean general circulation models (e.g., Danabasoglu and Marshall, 2007;Gent and Danabasoglu, 2011 ). The simulations described in Gent and Danabasoglu (2011) do show some eddy compensation, as a consequence of the dependence of the GM eddy transfer coefficient on Southern Ocean stratification. However, as discussed in Munday et al. (2013) , this mechanism precludes the model from achieving full eddy saturation.
Through the consideration of the eddy kinetic energy budget, Cessi (2008) proposes a mixing length based eddy parameterisation, with a GM eddy transfer coefficient depending on the ocean state and explicitly depending on the strength of the bottom drag. An approach also based upon consideration of the eddy kinetic energy budget is discussed in Eden and Greatbatch (2008) (see also Marshall and Adcroft, 2010 ), also employing a mixing length argument but utilising a local parameterised eddy kinetic energy budget to inform the magnitude and spatial structure of the resulting GM eddy transfer coefficient. While there is no conclusive observational evidence to suggest that the ocean is in an completely eddy saturated regime, there is ample evidence from mesoscale eddypermitting model experiments that coarse resolution models with current parameterisations appear unable to replicate eddy satura-tion in a self-consistent way (e.g., the work Fyfe et al., 2007 varies the GM eddy transfer coefficient manually with changing wind stress).
In Marshall et al. (2012) a geometric interpretation of the eddymean-flow interaction for the quasi-geostrophic equations was derived. A horizontally down-gradient closure for the horizontal eddy buoyancy fluxes leads to a GM eddy transfer coefficient of the form where E is the total (kinetic plus potential) eddy energy, and N/M 2 = T is an Eady time-scale which depends on the mean stratification, with N 2 = −(g/ρ 0 ) ∂ ρ z /∂z and M 2 = (g/ρ 0 ) |∇ H ρ z | , where g is the gravitational acceleration, ρ 0 is a reference density, ρ z is the mean density averaged at fixed height, and ∇ H is its horizontal gradient operator. A crucial point is that, if the eddy energy is known, there are no undetermined dimensional parameters; the only freedom is to specify the non-dimensional geometric parameter α of magnitude less than or equal to one (see, e.g., Bachman et al., 2017 ). A form similar to (1) also appears in Jansen et al. (2015) -implied by their Eqs. (9) and (11) -but with the eddy kinetic energy in place of the full eddy energy, and motivated by the inverse energy cascade being controlled by the rate of eddy energy generation through baroclinic instability as per Larichev and Held (1995) . However the form derived in Marshall et al. (2012) provides an explicit upper bound on the relevant geometric parameter α; no other dimensional scaling is possible provided the geometric parameter α is bounded away from zero. Moreover, here the eddy energy is determined prognostically via the solution of a dynamical equation which is coupled to the equations for the mean state. This article assesses the ability of the Marshall et al. (2012) GM eddy transfer coefficient to reproduce eddy saturation, via numerical calculations in an idealised, zonally averaged, two-dimensional ocean channel model. The idealised numerical model is motivated by the physical model discussed in Marshall et al. (2017) , where eddy saturation was demonstrated through considerations of the momentum and eddy energy budget, together with the scaling for the GM eddy transfer coefficient given by Eq. (1) . The ability of the Marshall et al. (2012) scheme to produce eddy saturation is compared against a number of alternative approaches, including approaches based upon mixing length arguments, and based upon the Visbeck et al. (1997) proposal. Since the Marshall et al. (2012) variant requires information about the eddy energy, the evolution of the mean state is coupled to a prognostic equation for the parameterised domain integrated eddy energy (cf. the local budget for the eddy kinetic energy in Eden and Greatbatch, 2008 ).
The paper proceeds as follows. In Section 2 the GM scheme and the Marshall et al. (2012) parameterisation variant are revisited, focusing in particular on the energetics of the problem, and providing physical and mathematical arguments as to why the Marshall et al. (2012) variant may be expected to have skill in producing emergent eddy saturation. Section 3 contains the details of the idealised numerical model and of the other parameterisation variants considered in this work. The implementation of the parameterisation variants and their results are presented in Section 4 for a case where the GM eddy transfer coefficient is assumed to be constant over the domain, and in Section 5 for a case where the GM eddy transfer coefficient is spatially varying, focusing on the case where a spatial structure depending upon the vertical stratification is enforced. The paper concludes in Section 6 , where the results are discussed, and a recipe for implementation in a global circulation models is proposed.

The Gent-McWilliams scheme and the energetic consequences
The GM scheme parameterises the effects of baroclinic eddies via the introduction of an adiabatic stirring of the mean density, acting to decrease the available potential energy of the system (e.g., Gent and McWilliams, 1990 ). Limiting consideration to the Boussinesq case, the mean density equation, zonally averaged at constant density ( Andrews, 1983;McDougall and McIntosh, 2001;Young, 2012 ), is where ρ # is the mean density field associated with averaging at constant density (see Young, 2012 ), with the thickness-weight averaged meridional velocity at constant density given by and w # defined such that (again following Young, 2012 ) Following McDougall and McIntosh (2001) , where u z is the velocity zonally averaged at constant height, and u * is the eddy transport velocity, with The GM scheme then takes the form where κ is the GM eddy transfer coefficient, and s is the slope of the mean density surfaces The energetic consequences of the GM scheme are as follows. Consider the zonally averaged hydrostatic Boussinesq equations in the form Here contributions from Reynolds stresses are neglected, it is assumed that all significant forcing F z and dissipation D z occurs in the zonal mean momentum equation, and ρ # is used in place of ρ z in the hydrostatic relation (consistent with the discussion in appendix B of McDougall and McIntosh, 2001 ). A budget for the mean energy may be obtained by multiplying by the mean velocity, integrating over the domain, using incompressibility and the mean density equation, and assuming that the normal components of both u z and u * vanish on all boundaries. The resulting budget becomes d d t The last term is a conversion term which, via substituting w * from Eq. (7) and performing an integration by parts, results in with horizontal and vertical buoyancy frequencies M and N respectively, where The final term in Eq. (11) is the conversion of eddy energy to mean energy. It follows that the eddy energy equation takes the form (see Appendix A for a more complete derivation) where it is assumed that the eddy energy source from external forcing is negligible compared to the eddy energy generation given in the first term on the right hand side of (13) . Here, ρ 0 E is the eddy energy density, and represents the dissipation of eddy energy, for example via topographic form stress. A simple model for this dissipation term is where λ is a dissipation time scale. The eddy energy budget (13) then becomes The first right-hand-side term in Eq. (13) , which is a stratification weighted integral of the GM eddy transfer coefficient, is a consequence of the GM scheme but is independent on the precise variant of the GM eddy transfer coefficient used.

Marshall et al. (2012) geometric framework and consequences
In Marshall et al. (2012) a geometric framework for the eddy fluxes is proposed. In particular a horizontally down-gradient closure for the horizontal eddy buoyancy fluxes yields where α is a non-dimensional geometric eddy efficiency parameter that is bounded in magnitude by one. Provided αN / M 2 is bounded away from zero and infinity, this implies that the magnitude of the GM eddy transfer coefficient should scale with the eddy energy E . This is the case if the mean density has a non-trivial gradient in both the horizontal and vertical directions, and if the geometric parameter α is bounded away from zero. Note that the dependence on the eddy energy is linear, as opposed to a square root dependence that is suggested by a mixing-length based argument (e.g., Cessi, 2008;Eden and Greatbatch, 2008 ). A linear dependence of the eddy energy may be obtained as in Jansen et al. (2015) if the length scale has a dependence on square root of the eddy energy also. With this form, once information about the eddy energy is known, for example from the solution of a parameterised eddy energy budget, then the only remaining freedom is in the specification of the non-dimensional geometric parameter α bounded in magnitude by one (see, e.g., Bachman et al., 2017 ). The physical implications of this closure are described in Marshall et al. (2017) . Here we highlight the relevant properties in terms of the expected scaling of the eddy energy on α and the dissipation, and further the implications for the scaling of the emergent zonal transport, eddy energy, and GM eddy transfer coefficient.
With Eq. (16) , the eddy energy budget (15) becomes where s = −M 2 /N 2 . In particular, in steady state, the balance holds. Note that, from thermal wind shear, and hence This is an expression of an eddy energy weighted balance between the eddy energy generation rate due to the eddies, given by the first integrand term, and the eddy energy dissipation rate, given by the second. The integral balance can be achieved if the vertical shear increases with the dissipation rate λ, and decreases with the geometric parameter α. Note that, following the argument of Straub (1993) , the zonal mean transport, defined as transport = and making the assumption that u z (−L z ) = 0 , scales with the vertical shear appearing as a factor in the first integrand term. Hence Eq. (20) suggests that the zonal transport defined in Eq. (21) scales with the dissipation rate λ, and scales inversely with the geometric parameter α, but not explicitly on the external forcing F z . For an appropriately smooth eddy energy the following scaling (see Appendix B for details) is further suggested, again indicating increased transport with increasing λ, and decreased transport with increased α, but not on the external forcing F z (cf. Marshall et al., 2017 ).
These scalings may be interpreted as follows. Increasing λ means the emergent eddy generation rate needs to increase to maintain the integral balance (20) , which is achieved via changes in the emergent stratification profile, resulting in steeper isopycnals and thus a larger transport. An analogous explanation for decreasing α suggests an increase in the transport.
Similar scalings of the emergent eddy energy and GM eddy transfer coefficient may be derived. Consider the mean energy equation along with (16) . At steady state and assuming the dissipation of the mean is small, the mean energy equation is For fixed λ, and assuming saturation such that u z and N do not depend on the forcing parameter, Eq. (23) On the other hand, the functional dependence of the emergent eddy energy and GM eddy transfer coefficient on varying λ and α is not so straightforward, since it is the vertical stratification weighted transport that is more directly influenced by these two parameters. The suggested dependencies and scalings for the emergent properties are then: (i) the transport (21) to be independent of the magnitude of forcing, increasing with increased dissipation and decreasing with increased α; (ii) GM eddy transfer coefficient κ to scale linearly with the magnitude of wind forcing; (iii) eddy energy level to increase linearly with the magnitude of wind forcing. These scalings are confirmed later via diagnosing the emergent properties from the simulation data.

Numerical implementation
The Marshall et al. (2012) variant for the GM eddy transfer coefficient given by Eq. (16) , together with the parameterised eddy energy budget in Eq. (15) , is implemented in a simplified two-dimensional model, similar to that employed in Marshall (1997) and Marshall and Radko (2003) (see also Gent et al., 1995 ). The channel model is described in Section 3.1 , and the other parameterisations to be tested against the Marshall et al. (2012) parameterisation are detailed in Section 3.2 .

Channel model
pansion coefficients, and T 0 and S 0 are a reference temperature and salinity respectively. The prognostic equation for the mean density is (cf. Gent et al., 1995 ) ∂ρ # ∂t The mean density is advected by the residual velocity ˆ u = (0 , ˆ v , w # ) T , which is the sum of the Eulerian circulation u z and the eddy induced transport velocity u * . The domain is chosen to be ( y, z ) ∈ (0 , L y ) × (−L z , 0) , and the equation is solved with no-normal-flow boundary conditions ˆ u · n = 0 on boundaries. The model is integrated in time until it reaches a steady state, with the convergence criterion to be defined. As a simple model for a forced-dissipative configuration, the Eulerian circulation appearing in the prognostic Eq. (24) is taken to satisfy the f -plane steady state equation and the thermal wind Eq. (19) , where τ s is the surface wind stress, and τ b is a representation of the bottom form stress (see Marshall, 1997 ). At the surface, the wind stress is with peak wind stress τ 0 , linearly decreasing to zero within the upper grid cell of the model. The bottom stress τ b is chosen to exactly cancel the local surface wind stress, i.e., τ b = τ s , applied within the bottom cell, representing the bottom form stress across topographic barriers ( Munk and Palmén, 1951 ). A choice of τ 0 effectively specifies the Eulerian circulation associated with the Deacon cell (cf. Marshall, 1997 ). The state ρ # determines the eddy induced transport velocity (0 , v * , w * ) T through Eq. (7) , which is then used to form the residual velocity to time step the prognostic Eq.
(24) . Assuming that u z vanishes at the sea floor due to topography blocking the flow, u z may be diagnosed via thermal wind shear relation (19) .
factor for ρ # (t = 0) given in (28) c 750 m e -folding depth for ρ # (t = 0) given in (28) The prognostic equations are discretised in space using a uniform resolution Arakawa C-grid ( Arakawa and Lamb, 1977 ) with y -and z -direction grid spacings of y and z respectively. The density ρ # is defined at the cell centres, fluxes and derivatives of ρ # on cell interfaces, with appropriate interpolation of the fields where required. The boundary conditions are implemented by setting boundary fluxes to zero. The forcing and dissipation are taken to be applied over the top and bottom cells. A fourth order Runge-Kutta method is employed to time step the prognostic Eq. (24) and eddy energy Eq. (15) , with a variable t chosen at the end of each time step so as to target a desired Courant number C ( Courant et al., 1928 ) For numerical stability in integrating the parameterised eddy energy equation, the variable time step is further restricted so that t ≤ 12 hours. The calculations are initialised with an exponential density profile where a , b , c are as in Table 1 , informed by the World Ocean Circulation Experiment ( Gouretski and Kolterman, 2004;Koltermann et al., 2011 ) data.
To avoid unbounded velocities associated with weak stratification, the slope tapering of Gerdes et al. (1991) with tapering function is employed, and it is ˜ s = F gkw (y, z) s that is used in computing the parameterised eddy transport velocity (7) . The advective form of the GM scheme is employed (e.g., Griffies et al., 1998;Griffies, 1998 ). In equations in which a division by the vertical stratification N 2 appears (e.g. the first right-hand-side term of Eq. (15) ) this is replaced with with a chosen value of the minimum vertical stratification N 2 min . Tests have shown the value N 2 min determines somewhat the eddy energy growth: too large a value and the eddy energy growth is inhibited, while too small a value results in very large eddy energy growth with non-negligible contributions from the model at depth.
During time stepping a basic convection scheme is applied, with each vertical water column sorted by density within each Runge-Kutta stage. The convection scheme facilitates the development of out-cropping at the surface, which would otherwise be constrained by the initially constant surface density and the nonormal-flow boundary condition. The convection scheme is disabled when Model parameter values are provided in Table 1 .

Alternative GM eddy transfer coefficients
For comparison, a number of alternative variants based on existing parameterisation schemes are also implemented in the idealised numerical model. A scheme that employs a mixing length assumption and has dependence on the eddy energy is given by where α ML is some non-dimensional parameter (without a formal bound) and L is a mixing length scale to be specified. This scheme has a weaker dependence on the eddy energy. An approach of this form is described in Eden and Greatbatch (2008) , where the eddy energy is replaced with the eddy kinetic energy, and the length scale is taken to be the minimum of the Rhines scale and the Rossby deformation radius (their Eq. 25 ). Setting the mixing length equal to the Rhines scale increases the eddy kinetic energy exponent to 3/4, and hence this is closer to the linear energy scaling in Eq. (16) . A similar mixing length approach is taken in Cessi (2008) where a statistically steady version of (15) is utilised to derive a form of the GM eddy transfer coefficient that has explicit dependence on the bottom drag. In Cessi (2008) , the eddy kinetic energy is used in place of the eddy energy, and L is chosen to be the Rossby deformation radius. Note that the derivation of Eden and Greatbatch (2008) , in their Eq. (26) , suggests that the GM eddy transfer coefficient should have a linear dependence on the eddy kinetic energy. However in their work the chosen length scale implicitly sets the magnitude of the eddy kinetic energy. Here, instead, the eddy energy is parameterised directly. In Jansen et al. (2015) a mixing length which scales with the square root of the eddy kinetic energy is discussed, yielding a form equivalent to the scaling of (16) , with the eddy kinetic energy again used in place of the eddy energy.
Based on instability arguments, Visbeck et al. (1997) proposed where α VMHS is some non-dimensional parameter (again without a formal bound). This variant has no explicit dependence on the eddy energy, and instead depends only on the mean stratification.
In Section 3 d of Visbeck et al. (1997) the length scale L is related to the grid scale, Rossby deformation radius, and the width of the baroclinic zone. Diagnosing diffusivities from a 4 °global numerical ocean model constrained using observation data and via an adjoint based optimisation, in Ferreira et al. (2005) , it is suggested that where κ 0 is some reference GM eddy transfer coefficient value, and S imparts a spatial structure to the GM eddy transfer coefficient that is dependent on the vertical stratification. The use of such a structure function results in a GM eddy transfer coefficient that is large towards the ocean surface whilst being small in the deep ocean where the stratification is weak. The reference value κ 0 is normally taken to be constant (e.g., Ferreira et al., 2005;Danabasoglu and Marshall, 2007;Gent and Danabasoglu, 2011 ).

Summary
In summary, the four variants for the GM eddy transfer coefficient considered in this article are:  Eden and Greatbatch (2008) and Cessi (2008) , denoted ML; • a scheme similar to that described in Visbeck et al. (1997) , denoted VMHS * .
Each of these four variants are considered subject to two approximations, with implementation details given in the appropriate sections. This first is where the GM eddy transfer coefficient is assumed to be spatially constant. The second is one where the GM eddy transfer coefficient has an imposed spatial structure set by S from Eq. (34) , to be in line with more modern numerical models (e.g., Danabasoglu and Marshall, 2007;Gent and Danabasoglu, 2011 ). Where relevant all length scales are set equal to the Rossby deformation radius L D = NH/ | f | . All the implemented variants are coupled to the parameterised eddy energy Eq. (15) , although this plays a prognostic role only for the GEOM and ML variants.
The use of a prescribed spatial structure for the GM eddy transfer coefficient contradicts somewhat with the original intention of the scheme described in Visbeck et al. (1997) . Thus a variant, denoted VMHS, is additionally considered, which uses the full local dependence as specified in Eq. (33) . Note, however, the length scale is still set equal to the Rossby deformation radius, which differs from the length scale used in Visbeck et al. (1997) .
The four parameterisation variants differ in their functional dependence on the eddy energy and the mean stratification, as summarised in Table 2 .

Implementation details
In this section the case of spatially constant GM eddy transfer coefficient is considered, employing the CONST, GEOM, ML and VMHS * variants described in Section 3.3 . The CONST variant is simply employed by taking a constant value of κ. To obtain a spatially constant GM eddy transfer coefficient for the GEOM variant with κ = αE(N/M 2 ) , the terms are appropriately re-arranged, and integrating over the domain leads to The domain integrated eddy energy E d y d z is computed by solving (15) .
For the ML variant, an analogous approach yields However the domain integral of the square root of the eddy energy is not available. Use of the Cauchy-Schwarz inequality (e.g., Eq. B.1 with p = q = 2 ) leads to and so the ML variant is implemented as Here a prescribed value for the new parameter α 1 is chosen.
For the VMHS * variant the form is used. The VMHS variant with fully local dependence on the mean state is considered in Section 5 . The initial state is spun up from rest first using the GEOM variant, with τ 0 = 0 . 2 N m −2 , λ = 2 × 10 −7 s −1 and α = 0 . 1 . The associated initial and equilibrium states are shown in Fig. 1 . The equilibrium state here has a transport of around 77 Sv and a domain average parameterised eddy energy (divided by the reference density ρ 0 ) of around 0 . 01 m 2 s −2 , the latter being similar to the level given in the observations of Meredith and Hogg (2006) . From this control run and taking the mixing length L to be the Rossby deformation radius L D = NH/ | f | for the ML and VMHS * variants, the emergent κ and end state ρ # are used to calibrate κ for CONST, α 1 for the ML variant in (36) and α 2 for the VMHS * variant in (37) , which are used for subsequent calculations where τ 0 and λ are varied. The GEOM parameter values used in the parameter sweep is given in Table 3 .

Table 3
Parameter values employed in the GEOM calculations. The value for the control run displayed in Fig. 1

Results
The transport associated with the equilibrium states with varying values for τ 0 and λ are shown in Fig. 2 . It is clear that CONST and VMHS * show significant sensitivity of the mean transport with respect to the peak wind stress. By contrast, the ML variant shows reduced sensitivity. Notably, the GEOM variant shows very low sensitivity to varying wind stress, and thus exhibits emergent eddy saturation. For varying eddy energy dissipation rate λ, CONST and VMHS * are by construction independent of λ, while the ML and GEOM variants show increased transport with increased dissipation. These observed behaviours are consistent with the analysis given in Section 2.2 . Denoting the domain average by the emergent κ and ⟨ E ⟩ are shown in Fig. 3 . The ML and VMHS * variants show a sub-linear dependence of the emergent κ on the peak wind stress τ 0 , while the GEOM variant exhibits an almost linear dependence. For the ML and GEOM variants the emergent κ decreases with increasing λ. The scaling of κ and ⟨ E ⟩ with peak wind stress τ 0 , for the GEOM variant, is consistent with the arguments given in Section 2.2 . It is found here that increasing the dissipation decreases the emergent eddy energy level.
The emergent eddy saturation property of the GEOM variant is not limited to this parameter set. Fig. 4 shows contour plots of the transport in ( τ 0 , λ) and ( τ 0 , α) parameter space. As expected, there is very little dependence of the transport on τ 0 and only at extreme parameter values is a variability seen in the contour plot. This shows robustness of the insensitivity to strength of peak wind over a range of parameters.
To show how the other emergent properties of the GEOM variant depend on α and λ, the transport, GM eddy transfer coefficient, and domain averaged eddy energy over ( λ, α) parameter space are shown in Fig. 5 . Increasing α reduces the mean transport as expected, from the discussion in Section 2.2 . The GM eddy transfer coefficient κ is found to increase with increasing α. The values of the emergent κ are consistent with the emergent transport, although large values are observed where the parameterised eddies are very efficient (small λ, large α). The eddy energy has a more complex dependence on α, but for weaker dissipation increasing α leads to a decrease in the eddy energy.

Implementation details
In this section a dependence of the GM eddy transfer coefficient on the vertical stratification is introduced, again with four variants based upon the CONST, GEOM, ML, and VMHS * discussed in Section 3.2 . The simplest CONST variant is now replaced with the form proposed in Ferreira et al. (2005) (39)  The values for the parameter sweep are given in Table 3 .
This imparts a vertical as well as horizontal spatial structure to the GM eddy transfer coefficient. Unlike Ferreira et al. (2005) , however, since the simple dynamical model employed here has no restoring boundary conditions at the surface to maintain a surface stratification, the use of N 2 ref at the ocean surface as in Ferreira et al. (2005) is perhaps not a suitable choice. Instead, N 2 ref is taken to be N 2 min here. With this, no tapering of S is employed. A Gaussian ten by ten point smoother with a three point standard deviation was applied to the density field ρ # before taking derivatives to form S, for numerical stability reasons.
The GEOM variant becomes κ = κ 0 S = αE(N/M 2 ) , where again α is a prescribed constant. Re-arranging, integrating over the domain, and now assuming that κ 0 is a constant in space leads to The domain integrated eddy energy E d y d z is computed by solving Eq. (15) as before.
For the ML variant, an analogous approach yields and use of the Cauchy-Schwarz inequality leads to Here a prescribed value for the parameter α 1 is again chosen.
For the variant based on Visbeck et al. (1997) , with the GM eddy transfer coefficient given by κ = κ 0 S = α VMHS (M 2 /N) L 2 , two forms are used. Assuming κ 0 is a constant in space results in the for the parameter sweep are given in Table 3 . Alternatively the form (33) may be used directly, resulting in the VMHS variant where now α 3 is a prescribed constant. This latter form introduces an additional explicit dependence on the local value of M and the local mixing length L . As for the previous constant GM eddy transfer coefficient case, the initial state is spun up from rest first using the GEOM variant, with τ 0 = 0 . 2 N m −2 , λ = 2 × 10 −7 s −1 and α = 0 . 1 . The initial state is the same one shown in Fig. 1 , and Fig. 6 shows the equilibrium stratification profile and the associated spatially varying GM eddy transfer coefficient. This equilibrium state here has a transport of around 66 Sv and a domain average parameterised eddy energy of around 0 . 009 m 2 s −2 . From this control run and taking the mixing length L to be the Rossby deformation radius L D = NH/ | f | as before for the ML and VMHS * variants, the emergent κ and end state ρ # are used to calibrate κ 0 for CONST in (39) , α 1 for ML in (41) and α 2 for the VMHS * in (42) , which are used for subsequent calculations where τ 0 and λ are varied. For the direct VMHS variant, the functional dependence of κ on M and N differs from the functional dependence specified by S. Here an initial value of α 3 in (43) is chosen manually so that a similar level of transport is obtained, though note that the results presented here are still slightly detuned. Other simulation details are kept the same as in Table 1 except that the convergence tolerance ξ 2 is now set to 5 × 10 −13 , as there is more variability given that κ is allowed to vary over space. No minimum or maximum values of κ are imposed. Calculations are restarted from a previously converged calculation for 500 model years, and those that do not converge within the period are excluded from the diagrams. The same values displayed in Table 3 are used for the parameter sweep.

Results
The resulting transport with varying τ 0 and λ for the five parameterisation variants is shown in Fig. 7 . The GEOM variant, though possessing a slight increase in transport as τ 0 is increased, again exhibits relative insensitivity of the transport to changes in wind forcing. It is notable that the ML variant also show a reduction of the sensitivity of the mean transport with respect to  the peak wind stress compared to its respective case with spatially constant GM eddy transfer coefficient. The GEOM variant once again exhibits a strong dependence on the eddy energy dissipation rate λ. The emergent eddy saturation for the GEOM variant is again found to be robust over a range of parameters, as shown in Fig. 9 . Fig. 10 shows contour plots of the emergent properties with varying λ and α. In both figures, non-converged states have been greyed out. Although showing much more variability than the analogous spatially constant κ case in Fig. 5 , there is a pattern of increased transport at increasing λ or decreasing α, and of decreased κ 0 at increasing λ or decreasing α. Note the region with low λ and large α has small transport, large ⟨ E ⟩ and thus large κ max . The resulting parameterised eddies in the large λ and small α region are very weak, and it may be seen that the oscillations of the mean state appears to persist and have not been deemed to converge according to the imposed criterion.

Summary
In this article the problem of emergent eddy saturation in coarse resolution ocean modelling with parameterised mesoscale eddies has been considered. Specifically, an idealised zonally averaged channel configuration was used to test the sensitivity of mean zonal transports with respect to the strength of surface wind forcing, and additionally with respect to the strength of total eddy energy dissipation and closure parameters. Variants of the Gent and McWilliams (1990) scheme have been tested, with a constant GM eddy transfer coefficient, a GM eddy transfer coefficient with a stratification dependence based upon that described in Visbeck et al. (1997) , a GM eddy transfer coefficient with a mixing-length inspired energy dependence (e.g., Eden and Greatbatch, 2008;Cessi, 2008 ), and a GM eddy transfer coefficient derived from the geometric framework described by Marshall et al. (2012) . For the schemes with eddy energy dependence a parameterised equation for the domain integrated eddy energy was solved. By integrating over the domain, specific closures were derived, falling into two classes -one where the GM eddy transfer coefficient was spatially constant, and one where the GM eddy transfer coefficient had a spatial structure based upon that described in Ferreira et al. (2005) . A form with additional stratifica- Fig. 9. Contour plot of the emergent transport over ( a ) ( τ 0 , λ) space (with α = 0 . 1 ), and ( b ) ( τ 0 , α) space (with λ = 2 × 10 −7 s −1 ≈ 0 . 017 day −1 ). Regions with non-converged solutions have been greyed out. The values for the parameter sweep are given in Table 3 .  Table 3 . Regions with non-converged solutions have been greyed out. tion dependence, closer to the original proposal of Visbeck et al. (1997) , was additionally tested.
It was found that the scheme derived from the geometric framework of Marshall et al. (2012) led to almost complete emergent eddy saturation, with little or no significant dependence of the mean transport on the surface wind stress magnitude. This lack of dependence was additionally observed for a wide range of eddy energy dissipation time scales and parameterisation parameter values. Moreover, it was found that the changes to the equilibrium stratification profile with different values of peak wind stress were small (not shown). Furthermore, the dependence of the transport and other emergent quantities are consistent with the physical and mathematical arguments given in Section 2.2 . On the other hand, the use of a basic spatially and temporally constant GM eddy transfer coefficient led to a very significant dependence of the mean zonal transport with respect to the wind stress, similar to behaviour reported in low resolution ocean model tests described in Munday et al. (2013) . Variants based upon the Visbeck et al. (1997) and upon mixing length arguments were generally found to have a somewhat reduced sensitivity, but did not exhibit full eddy saturation.

Discussion and future work
This work focuses on eddy saturation, but an equally important process that has not been investigated in this work is the ability of the GM eddy transfer coefficient variants in showing eddy compensation. In particular, the extent of eddy compensation depends upon both the magnitude and the spatial structure of the eddy induced transport, and the degree to which it cancels with the local Eulerian circulation ( Meredith et al., 2012 ). The model considered in this article has no representation for ocean basins and hence is unsuitable for studying eddy compensation. An investigation into the ability of the Marshall et al. (2012) variant of the GM eddy transfer coefficient in showing emergent eddy compensation would require a more sophisticated eddy energy budgets than the one employed here, and is left as future work.
Assuming that the eddy energy is given via a parameterised eddy energy budget, the only remaining freedom in the Marshall et al. (2012) variant is in the specification of the non-dimensional geometric parameter α, as all dimensional information on the magnitude of the GM eddy transfer coefficient is already provided by the eddy energy and mean stratification. In this work α was chosen to have a constant value of 0.1, which was guided by the diagnoses of the equilibrated states in a three-layer wind forced quasi-geostrophic double gyre simulation ( Marshall et al., 2012 ) and an Eady spindown simulation of the hydrostatic primitive equations ( Bachman et al., 2017 ). In diagnostic calculations α is not a constant, and in particular α was found in Bachman et al. (2017) to vary depending on whether the system is in a linear growth phase or in later phases of the spindown evolution. It is perhaps of theoretical interest to have α evolving in time to capture the initial instability, finite-amplitude regime, and transition into an equilibrated state, although this is beyond the scope of the current work.
In this paper we have found that the functional dependence for the GM eddy transfer coefficient proposed in Marshall et al. (2012) , which incorporates energetic constraints through the solution of a parameterised eddy energy budget, yields near total emergent eddy saturation in a highly idealised configuration. While the degree of saturation in the ocean is not known, numerical models do appear to support this dynamically interesting regime, and this parameterisation variant is able to show emergent eddy saturation. A clear extension would be to implement and test this scheme in a global ocean model. Since the GM scheme is normally already built into global ocean circulation model as a core component, it would appear the main additional challenge would be (i) to add a parameterised eddy energy budget that couples with the GM scheme, and (ii) derive an appropriate form for a local parameterised eddy energy budget. The domain integrated eddy energy budget employed here is much too restrictive for use in a global ocean model. We envisage the scheme may be implemented into an operational global circulation model as follows: 1. Solve for the provisional eddy transport velocities, with a preferred vertical profile for the eddy transfer coefficient, utilising the standard GM scheme; 2. Vertically integrate the implied eddy form stress and compare with the theoretical prediction derived from the Marshall et al. (2012) geometric framework, using a prescribed nondimensional parameter α; 3. Solve for the parameterised, vertically integrated eddy energy budget, analogous to Eden and Greatbatch (2008) but for the full, rather than kinetic, eddy energy; 4. Rescale the eddy transport velocities, equivalent to rescaling the GM eddy transfer coefficient, uniformly over the vertical column such that each vertical integral of the eddy form stress matches the theoretical prediction from the Marshall et al. (2012) geometric framework.
By applying the energetic constraint in the vertical integral of the eddy form stresses, the recipe given above succeeds in retaining the positive-definite conversion of mean to eddy energy associated with the GM scheme, as well as the derived energetic constraint given in the Marshall et al. (2012) geometric framework.
In a closure for ocean turbulence one must typically tune the closure parameters in order to match a desired large-scale or mean state of interest. However for many key questions in physical oceanography, it is not only the mean state itself, but also the sensitivity of that mean state to external changes, which is of interest. This is, for example, critical to the understanding of the long time response of the ocean and broader climate system to long term forcing changes. The Gent-McWilliams closure is now a key component in large scale climate relevant ocean modelling, but it has been found that existing variants of the scheme in wide use, in particular with a constant Gent-McWilliams eddy transfer coefficient, do not yield accurate representations of ocean transport sensitivities with respect to changed in wind forcing (e.g., Farneti and Gent, 2011;Gent and Danabasoglu, 2011 ). This work provides the first evidence that the phenomenon of eddy saturation may be captured without major changes to the existing Gent-McWilliams closure, simply by employing the Marshall et al. (2012) form for the GM eddy transfer coefficient, derived from first principles with no tunable dimensional parameters, coupled with a parameterised eddy energy budget. A proposal on how this scheme may be implemented into a global circulation model via the addition of a parameterised eddy energy equation has been given here. Investigations into implementing this into a general circulation model, as well as theoretical developments for a parameterised eddy energy budget, are under investigation.
with D being a weak derivative, d is the dimensionality of the domain, 1 ≤ r , q ≤ ∞ , j / m ≤ a ≤ 1, and s > 0 is arbitrary. This does not cover some exceptional cases, though they are not of interest here. The constants C 1, 2 only depend on the domain and the choice of the parameter values. For d = 2 here, taking j = 0 , p = 2 and s = 1 , it is noted that m = 1 , r = 1 and a = 1 is one option (which is a form of the Sobolev inequality; e.g., Evans 1998 , Section 5.6.1), and taking s = 1 and f = E results in ∥ E ∥ L 2 ≤ C 1 ∥ D E ∥ L 1 + C 2 ∥ E ∥ L 1 .

(B.5)
If m = 1 , r = 2 , a = 1 / 2 and q = 1 instead, then (B.6) which is analogous to the inequality of Nash (1958) . Other possibilities exist involving higher derivatives. Either way, assuming that the terms involving the derivatives are small compared to C 2 ∥ E∥ L 1 , then the relation (22) again follows, with a constant of proportionality that only depends on the domain and the parameter values chosen in the Gagliardo-Nirenberg interpolation inequality and is bounded away from zero and infinity. Again, this implies a lower bound on the weighted norm on the right hand side.