Implementation of a geometrically and energetically constrained mesoscale eddy parameterization in an ocean circulation model

The global stratification and circulation of the ocean and their sensitivities to changes in forcing depend crucially on the representation of the mesoscale eddy field. Here, a geometrically informed and energetically constrained parameterization framework for mesoscale eddies --- termed GEOMETRIC --- is proposed and implemented in three-dimensional primitive equation channel and sector models. The GEOMETRIC framework closes mesoscale eddy fluxes according to the standard Gent--McWilliams scheme, but with the eddy transfer coefficient constrained by the depth-integrated eddy energy field, provided through a prognostic eddy energy budget evolving with the mean state. It is found that coarse resolution calculations employing GEOMETRIC broadly reproduce model sensitivities of the eddy permitting reference calculations in the emergent circumpolar transport, meridional overturning circulation profile and the depth-integrated eddy energy signature; in particular, eddy saturation emerges in the sector configuration. Some differences arise, attributed here to the simple prognostic eddy energy budget employed, to be improved upon in future investigations. The GEOMETRIC framework thus proposes a shift in paradigm, from a focus on how to close for eddy fluxes, to focusing on the representation of eddy energetics.


Introduction
Accurate representation of the mesoscale eddy field and its feedback onto the mean ocean state is one of the most pressing challenges for ocean modelling, in particular in the ocean circulation models used for climate prediction, which often lack explicit representation of the mesoscale eddy field. Over the past two decades a widely adopted approach is due to Gent and McWilliams (1990, hereafter GM). The GM scheme parameterizes eddies through both a diffusion along neutral surfaces (Redi 1982) and an eddy-induced circulation that acts to flatten neutral density surfaces (Gent et al. 1995;McDougall and McIntosh 2001), thereby extracting available potential energy from the mean state. The adoption of GM immediately resolved a number of known deficiencies in ocean circulation models by removing the spurious diapycnal water mass conversions that were prevalent in the existing eddy parameterization schemes (Danabasoglu et al. 1994).
A known deficiency of the GM eddy parameterization is the very different response of the Southern Ocean circulation to changes in surface wind stress in models with GM and explicit eddies. With a spatially constant eddy diffusivity, the circumpolar transport increases with the strength of the surface wind forcing, whereas little sensitivity is observed in the equivalent models with explicit eddies (e.g., Munday et al. 2013;Farneti et al. 2015). This is known as eddy saturation (Hallberg and Gnanadesikan 2001) and was first predicted on theoretical grounds by Straub (1993). Eddy saturation is generally found in models that partially resolve a mesoscale eddy field (e.g., Hallberg and Gnanadesikan 2006;Hogg and Blundell 2006;Hogg et al. 2008;; Morrison and Hogg 2013;Munday et al. 2013;Hogg and Munday 2014) but not in models where eddies are parameterized by the GM scheme where the eddy transfer coefficient is constant in space and time (e.g., Munday et al. 2013;Farneti et al. 2015;Mak et al. 2017).
A further discrepancy between eddy permitting and coarse resolution models is the reduced sensitivity of the time-mean residual meridional overturning circulation to changing wind forcing obtained in eddy permitting models (e.g., Meredith et al. 2012;Viebahn and Eden 2012;Morrison and Hogg 2013;Munday et al. 2013;Hogg and Munday 2014;Farneti et al. 2015). This is known as eddy compensation (Viebahn and Eden 2012). Eddy compensation is less well understood than eddy saturation, depending in subtle ways on the vertical structure of the eddy response to changes in surface forcing (e.g., Morrison and Hogg 2013). The response is further complicated by the fact that the residual meridional overturning circulation is affected by, for example, bathymetric details (e.g., Hogg and Munday 2014;Ferrari et al. 2016;de Lavergne et al. 2017). Generally, it is found that eddy permitting calculations are strongly eddy saturated and partially eddy compensated (e.g., Munday et al. 2013;Farneti et al. 2015).
In contrast, eddy saturation and compensation are not well represented in models that parameterize eddies through GM with an eddy transfer coefficient that is constant in space and time.
Partial eddy saturation and compensation can be obtained when the eddy transfer coefficient is allowed to vary in space and time (e.g., Gent and Danabasoglu 2011;Hofman and Morales Maqueda 2011;Farneti et al. 2015), due to the nonlinear dependence of the eddy transfer coefficient on the mean density gradients, in particular downstream of major bathymetric features. Numerous papers have attempted to derive the functional dependence of the eddy transfer coefficient on the ocean state as a function of space and time, from first principles (e.g., Treguier et al. 1997;Visbeck et al. 1997) and via diagnoses of numerical simulations (e.g., Ferreira et al. 2005;Ferrari et al. 2010;Bachman and Fox-Kemper 2013;Mak et al. 2016;Bachman et al. 2017).
In Eden and Greatbatch (2008) it was instead proposed that relating the eddy transfer coefficient to the eddy kinetic energy through a mixing length argument (see also Cessi 2008, Marshall and Adcroft 2010and Jansen and Held 2013. This approach requires solving for the eddy kinetic energy through a prognostic eddy energy budget. More recently, Marshall et al. (2012) have developed a new framework, here termed "GEOMETRIC", in which the inferred GM eddy transfer coefficient is entirely determined by the total eddy energy, the stratification, and an unknown non-dimensional parameter that is bounded in magnitude by unity. The predicted eddy transfer coefficient broadly agrees with that diagnosed in eddy resolving calculations .
Moreover, the GEOMETRIC eddy transfer coefficient leads to eddy saturation when implemented in an idealized two-dimensional model of the Antarctic Circumpolar Current .
The aims of this work are to: 1. implement GEOMETRIC in a three-dimensional ocean circulation model; 2. diagnose the extent to which GEOMETRIC reproduces eddy saturation and eddy compensation as obtained in the eddy permitting calculations; 3. explore the spatial variations of the eddy energy in the eddy permitting calculations and the extent to which these are reproduced with GEOMETRIC.
The article proceeds as follows. Section 2 outlines the GEOMETRIC approach and section 3 discusses the associated parameterized eddy energy budget. Implementation details relating to the parameterization schemes considered in this article are given in Section 4 Results from idealized channel and sector configurations are detailed in Section 5 and 6 respectively, with the model numerics described within the sections. The article summarizes and concludes in Section 7, and discusses further implementation challenges and research directions.

GEOMETRIC
GEOMETRIC ("Geometry of Ocean Mesoscale Eddies and Their Rectified Impact on Climate") represents a framework for parameterizing mesoscale eddies that preserves the underlying symmetries and conservation laws in the un-averaged equations of motion. GEOMETRIC was originally derived under the quasi-geostrophic approximation (Marshall et al. 2012) although elements of the framework generalize to the thickness-weighted averaged primitive equations (Maddison and Marshall 2013).
There are two fundamental ingredients in GEOMETRIC: 1. representation of the eddy-mean flow interaction through an eddy stress tensor, which can be bounded in terms of the total eddy energy (Marshall et al. 2012); 2. solution of a consistent eddy energy equation (cf. Eden and Greatbatch 2008;Cessi 2008;Marshall and Adcroft 2010).
Crucially, given knowledge of the total eddy energy and mean stratification, all of the remaining unknowns are dimensionless, i.e., there is no freedom to specify dimensional quantities such as eddy length scales or eddy diffusivities.
In the simplest limit, in which the lateral eddy Reynolds stresses are neglected, GEOMETRIC reduces to GM, with the eddy transfer coefficient given by where E is the total eddy energy, N = (∂ b/∂ z) 1/2 is the buoyancy frequency, M 2 = |∇ H b| is the magnitude of the lateral buoyancy gradient, b is buoyancy and ∇ H is the horizontal gradient operator. Here the overbar represents a time filter applied at fixed height. An equivalent form of this was given in Jansen et al. (2015), obtained through combining a mixing length argument with those of Larichev and Held (1995), but with the eddy kinetic energy in place of the total eddy energy. Once the eddy energy field is known, the only freedom then is in the specification of the non-dimensional eddy efficiency parameter α, satisfying |α| ≤ 1 in the quasigeostrophic limit (Marshall et al. 2012). This α parameter can be diagnosed from eddy permitting simulations: results from wind-driven gyre calculations in a quasigeostrophic model (Marshall et al. 2012) and nonlinear Eady spindown calculations in a primitive equation model  suggest that typically α = O(10 −1 ).
The efficacy of GEOMETRIC has been established through three proofs of concept: • in the linear Eady (1949) model of baroclinic instability, an analytical test case, GEOMETRIC produces the correct dimensional energy growth rate Marshall et al. (2012); • in the fully-turbulent nonlinear Eady spin-down problem, as simulated by Bachman et al. (2017), the eddy transfer coefficient predicted by GEOMETRIC, (1), gives good agreement with those diagnosed from the numerical calculations, across four orders of magnitude of the eddy transfer coefficient; • when applied to a two-dimensional model of the Antarctic Circumpolar Current with a domainintegrated eddy energy budget , GEOMETRIC produces eddy saturation, i.e., a circumpolar volume transport that is insensitive to the surface wind stress, due to an interplay with the zonal momentum budget and eddy energy budget , the essential components of which are preserved by GEOMETRIC.

Eddy energy equation
The outstanding challenge is then to solve for the eddy energy field. Solution of a prognostic equation for the kinetic eddy energy in three dimensions has been attempted by Eden and Greatbatch (2008).
In GEOMETRIC, the total eddy energy is required. In this paper, it is proposed that the depth-integrated eddy energy is solved for, as this offers a number of advantages: (i) the conceptual and logistical simplicity of working in two rather than three dimensions; (ii) avoidance of division by zero in (1) when the isopycnals are flat; (iii) retention of desirable properties of GM such as the positive-definite sink of available potential energy (Gent and McWilliams 1990;Gent et al. 1995) that are instrumental in its robustness.
The consequence is that the eddy diffusivity is energetically constrained in the vertical integral only. Specifically, suppose the eddy transfer coefficient varies in the vertical according to where Γ(z) is a prescribed dimensionless structure function, (e.g., Γ(z) = N(z) 2 /N 2 ref , Ferreira et al. 2005). Then the proposed eddy transfer coefficient is which is to be coupled to a parameterized budget for the depth-integrated eddy energy, E dz. Rather than derive a depth-integrated eddy energy budget from first principles, which contains terms that are unknown, the following heuristic approach is taken. The primary source of eddy energy is baroclinic instability (Charney 1948;Eady 1949), and must be incorporated in a manner that is consistent with the loss of mean energy due to the slumping of density surfaces as represented by GM.
In observations, the eddy energy is observed to propagate westward, at roughly the intrinsic long Rossby phase speed (Chelton et al. 2007(Chelton et al. , 2011, with an additional advective contribution that is adequately modelled by the depth-mean flow (Klocker and Marshall 2014). In this paper, the contribution from the intrinsic long Rossby phase speed is not included, while recognising that it will be required in an eventual implementation of GEOMETRIC in a global circulation model. Previous studies indicate that the lateral redistribution of eddy energy is not required to obtain eddy saturation with GEOMETRIC Mak et al. 2017), but it will likely affect the detailed response.
A Laplacian diffusion of eddy energy is incorporated following Eden and Greatbatch (2008). While included as a stabiliser, there are indications that the use of a Laplacian diffusion corresponds to the divergence of the mean energy flux in an f -plane barotropic model of turbulence (Grooms 2015).
The dissipation of eddy energy is complicated, involving a myriad of processes.
These include: bottom drag (e.g., Sen et al. 2008); lee wave radiation from the sea floor (e.g., Naveira Garabato et al. 2004;Nikurashin and Ferrari 2011;Melet et al. 2015); western boundary processes (Zhai et al. 2010); and loss of balance (e.g., Molemaker et al. 2005). Moreover, the eddy energy dissipation through these various processes will critically depend on the partition between eddy kinetic and eddy potential energy, and the vertical structure of the eddy kinetic energy (Jansen et al. 2015;Kong and Jansen 2017). Each of these required detailed investigation. Instead, a simple approach is followed here, representing eddy energy dissipation through a linear damping at a rate λ , recognising that λ parameterizes all of the physics outlined above.
To summarize, the proposed parameterized eddy energy budget is where: u z is the depth-averaged flow; e x is the unit vector in the longitudinal direction, and c is the intrinsic long Rossby phase speed which varies with latitude (Chelton et al. 2007(Chelton et al. , 2011; κ gm M 4 /N 2 is the eddy energy source, equal to the release of mean potential energy by slumping of density surfaces via GM; λ is the linear damping coefficient for the eddy energy that could in principle be a function of space and time; η E is the coefficient of the Laplacian diffusion of the depthintegrated eddy energy.

Experimental design
As a first step, the following simplifications are made: κ gm is taken to be vertically constant (so Γ(z) ≡ 1); c is set to zero; λ is a control parameter that is a constant in space time. Three sets of experiments are considered, as described by the following.

1) GEOM loc
The first set of experiments employ GEOMETRIC locally in latitude and longitude, as detailed in Section 2 and 3, with the eddy transfer coefficient computed as (3), coupled to the parameterized eddy energy budget (4). The GEOM loc scheme is implemented wholly within the GM/Redi package within MITgcm (Marshall et al. 1997a,b). First, all spatial derivatives of the mean density field are passed through a five point smoother, and κ gm is calculated according to (3) with the smoothed M 2 /N (the vertical integral of M 2 /N is bounded below by a small number to prevent division by zero). Then, if desired as a precaution to prevent possible large eddy induced velocities, κ gm may be capped below and above by some κ min and κ max , and the GM/Redi tensor is formed and passed through a slope tapering/clipping scheme. With the resulting κ gm , the eddy energy budget (4), discretized in space by a centered second-order differencing, is time stepped with a third order Adams-Bashforth scheme (started with forward and second order Adams-Bashforth steps), again with the smoothed M 2 /N.

2) GEOM int
In the second set of experiments, GEOM int , as previously considered in Mak et al. (2017), the eddy energy budget is integrated in space. With x as longitude and y as latitude, the eddy transfer coefficient is calculated as and is coupled to the parameterized eddy energy budget given by where there are no longer advective contributions to the eddy energy tendency, and the Laplacian diffusion of eddy energy has been removed. The eddy transfer coefficient (5) is now a constant in space but may vary in time, and the eddy energy budget (6) becomes an ordinary differential equation. The GEOM int scheme was also implemented wholly within the GM/Redi package in MITgcm following analogous steps, except the eddy energy budget (6) is time-stepped by a backward Euler scheme for numerical stability. The same κ min , κ max and slope tapering scheme as GEOM loc were used.

3) CONST
Finally, a control case with the standard GM scheme and a constant, prescribed eddy transfer coefficient is considered. The emergent parameterized eddy energy does not affect any of the resulting dynamics, and the routines for time stepping the parameterized eddy energy budget are bypassed.
The coarse resolution calculations GEOM int , GEOM loc and CONST are compared to calculations from eddy permitting reference calculations (REF). To assess the performance of the parameterization variants, various diagnoses of the resulting time-averaged data are presented; unless otherwise stated, all subsequent figures and statements refer to the time-averaged data. No mixed layer schemes are employed in the calculations.
The theory behind GEOMETRIC applies to the GM eddy transfer coefficient and not to the enhanced eddy diffusion of tracers along isopycnals (e.g., Redi 1982). While GM and Redi diffusion are often implemented in the GM/Redi tensor together (e.g., Griffies 1998;Griffies et al. 1998), the corresponding coefficients need not be the same. In all calculations presented here, the Redi diffusion coefficient is prescribed to be κ redi = 200 m 2 s −1 , and the GM eddy transfer coefficient follows the prescription of GEOM int , GEOM loc or CONST as appropriate. The consequences of the parameterization choices as well as the simplifications made for this work are discussed at the end of the article.

a. Setup and diagnostics
As an extension of the f -plane, zonally-averaged channel model of the Antarctic Circumpolar Current presented in Mak et al. (2017), an idealized channel configuration on a β -plane is considered in MITgcm. The configuration is essentially a shorter version of the channel configuration reported in Munday et al. (2015) and , with no continental barriers. The domain is 4000 km long, 2000 km wide and with a maximum depth of 3000 m. The model employs a linear equation of state with temperature only, and with an implicit free surface. A ridge with height of 1500 m and width 800 km blocks f /H contours and allows for the topographic form stress to balance the surface wind stress (e.g., Munk and Palmén 1951;Johnson and Bryden 1989).
An idealized zonal wind stress of the form is imposed, where L y is the meridional width of the channel, and τ 0 is the peak wind stress. The temperature is restored to the linear profile with ∆T = 15 K on a time-scale of 10 days over the top cell of height 10 m. The vertical temperature diffusion has magnitude κ d = 10 −5 m 2 s −1 , except in a tapered sponge region to the north of width 150 km where the vertical temperature diffusion is increased sinusoidally to κ d = 5 × 10 −3 m 2 s −1 to maintain a non-trivial stratification and energize the eddies (e.g. Hogg 2010; Munday et al. 2015).
A linear bottom drag with coefficient r is applied in the deepest level above the bathymetry. The vertical domain is discretized with 30 uneven vertical levels, thinnest of 10 m at the top, down to the thickest of 250 m, and with a partial cell representation of the bathymetry. A staggered baroclinic time stepping scheme was employed. See Munday et al. (2015) and  for further model details.
For the eddy permitting reference calculations (REF), the horizontal grid spacing is uniform at 10 km. A control simulation with control peak wind stress τ 0 = τ c = 0.2 N m −2 and control bottom drag coefficient r = r c = 1.1× 10 −3 m s −1 is carried out for 400 model "years" (360 days), from which perturbation experiments at varying τ 0 and r were carried out for a further 200 model years, before averages were taken for a further 20 model years. The eddy permitting calculation employs the full Leith viscosity (e.g., Fox-Kemper and Menemenlis 2008) with a coefficient of 2.
For the coarse resolution models, the horizontal grid spacing is mostly at 100 km, except at the northern boundary where the grid spacing is 50 km so as to have at least three grid points over the northern sponge region. A control GEOM int calculation with τ c and λ = λ c = 10 −7 s −1 (consistent with observation-constrained estimates in the Southern Ocean, Marshall & Zhai, pers. comm.; see also Melet et al. 2015) is first carried out over 500 model years. Perturbation experiments in GEOM int , GEOM loc and CONST at varying τ 0 and λ (rather than r) are then restarted and carried out for a further 300 model years, and averages are taken for a further 200 model years. The coarse resolution calculations employ a harmonic friction in the momentum equation that forces the grid scale Reynolds number to be 0.0075. The calculations have α and κ 0 tuned so that the circumpolar transport at the control parameter values matches the control calculation of REF, after which they are fixed for the perturbation experiments.
Note that in REF, the control parameters are τ 0 and r, while in the coarse resolution calculations, the control parameters τ 0 and λ , and r is fixed. The relevant parameter values are documented in Table 1.
Several diagnostics are computed to compare mean properties of the parameterization variants GEOM int , GEOM loc and CONST, against the reference calculation REF. The total circumpolar transport given by where L x is the length of the circumpolar channel, and (·) denotes a time filter performed at fixed height. Another is the thermal wind transport, given by the thermal wind velocity is given by assuming that u therm (z = −H) = 0, with ρ obtained from the temperature via the linear equation of state. Then, the complement transport to thermal wind transport is defined to be Note that the definition of the transport decomposition employed here differs from that employed in Munday et al. (2015). There, the bottom flow transport is defined to be equal to the bottom flow multiplied by the depth (thus sensitive to bottom flow details), and the baroclinic transport then is the remaining component of the total transport. It should be noted that the diagnosed values of the baroclinic transport as defined in Munday et al. (2015) and the thermal wind transport defined in (11) differ only very slightly in these channel experiments. Finally, following the definition of Gnanadesikan (1999) (see also Abernathey and Cessi 2014), a thermocline location is obtained by computing essentially a center-of-mass calculation for the vertical coordinate z, and this quantity is averaged over the northern sponge region where the thermocline is deepest. This provides a measure of the model stratification, with a deeper thermocline (more negative z therm ) expected to correlate with increased thermal wind transport.

b. Results
The diagnosed results are presented in Figure 1. As a summary, in this channel set up, the total transport of REF decreases with increasing wind, and increases with increased linear bottom drag. The total transport is composed principally of transport due to thermal wind. The complement component to the thermal wind transport increases with increased wind, and decreases slightly with increased bottom linear drag. The changes in the thermal wind transport are reflected in the resulting z therm , where a deeper thermocline corresponds to a larger spatial extent of the thermal wind. With this in mind, the CONST calculations display the opposite sensitivity to REF in terms of the dependence of the total circumpolar transport, thermal wind transport and the resulting thermocline location on the peak wind stress magnitude. On the other hand, both GEOM int and GEOM loc capture the changes in thermal wind transport and thermocline location with increasing wind stress. While there is a degree of tuning for the range of λ exployed here, the GEOM int and GEOM loc simulations at these values of λ results in similar trends to REF for the total and thermal wind transport.

1) VARYING WIND EXPERIMENTS
First, it is interesting to note that, even in REF, the total transport decreases with increasing wind, and the transport is non-zero at zero wind. The latter is due to the northern sponge region with enhanced vertical temperature diffusivity, which acts to maintain a stratification at depth and, together with surface restoring of temperature, results in tilting isopycnals and thus a thermal wind transport (e.g., Morrison and Hogg 2013). In this model, the model thermocline becomes shallower with increasing wind. As a result, the geostrophic flow occupies a smaller volume even though the peak geostrophic flow speed may be larger, thus resulting in a smaller integrated thermal wind transport. The decreased thermocline depth with increasing wind is likely due to the choice of imposing the northern sponge region condition; such behavior is not observed when a fully dynamical basin sets the northern channel stratification (as in the sector profile in the next section) or when the northern boundary temperature is relaxed to a prescribed profile (as in, e.g., Abernathey and Cessi 2014, where they employ instead a flux boundary condition at the ocean surface).
Despite the perhaps unexpected sensitivity to changing wind forcing in REF, it is encouraging to see that both the GEOM int and GEOM loc are able to reproduce the analogous sensitivities, particularly in the thermal wind transport and thermocline location diagnostic. In contrast, the standard CONST variant displays opposite sensitivity in the transport and thermocline location. Figure 2 shows the emergent zonally averaged temperature profile and the effective eddy induced overturning that counteracts the Eulerian overturning, resulting in steeper isopycnals. This leads to increased thermal wind transport, and is consistent with the diagnostics displayed in Figure 1(b, d). This feature of increased thermal wind transport is reproduced by the GEOM int and GEOM loc calculations, and is consistent with the findings of Mak et al. (2017).

3) OTHER EMERGENT QUANTITIES
The emergent eddy energy level and GM eddy transfer coefficient κ gm have also been diagnosed. Figure 3 shows the domain-averaged eddy energy E and domainaveraged GM eddy transfer coefficient κ gm . While the GEOM int and GEOM loc calculations have a value of the total eddy energy from the parameterized eddy energy budget, the depth-integrated value of the total eddy energy 1 for REF and CONST is calculated by diagnosing the sum of the depth-integrated (specific) eddy kinetic energy Here, u = u+u ′ and z = z+z ′ , where the time filter on the latter is to be carried out in density co-ordinates, so z ′ is the deviation from the mean isopycnal height. The thicknessweighted averaging is carried out here with the layers package in MITgcm (e.g., Abernathey et al. 2011). For this channel configuration with a linear equation of state for temperature, the calculation is carried out in temperature co-ordinates, with temperature referenced to the top model level at the surface, with binning over 81 discrete layers between −4 • C and 16 • C, equally spaced at 0.25 • C. In the eddy permitting calculations, it is the EPE contributions that dominate, accounting for around 90% of the total eddy energy; at the highest wind forcing, EKE accounts for about 12% of the total eddy energy, and decreases to about 5% for the largest value of linear bottom drag coefficient employed. For GEOM int and GEOM loc , the emergent E increases with increasing wind stress, though not necessarily at the same rate as REF. The rate of increase for GEOM int is slightly sub-linear, as opposed to the predicted linear scaling given in Mak et al. (2017). The increase in κ gm is also slightly sub-linear, consistent with the behavior of E . More variation is shown in GEOM loc in both the resulting E and κ gm levels, though the increase roughly follows that of GEOM int . It should be noted that while κ gm ≤ κ max in GEOM loc , locally κ max does get applied to the emergent κ gm albeit in a small region of the domain (where the parameterized eddy energy is large, see Figure 4c, d). At the lower peak wind stress values, the emergent eddy energy value from GEOM int and GEOM loc is much smaller than REF, which is consistent with the circumpolar transport of GEOM int and GEOM loc being larger than REF in Figure 1(a) as a result of the reduced emergent κ gm . For CONST, the diagnosed E is roughly two orders of magnitude smaller, and so appears almost on the axes in this plot with linear scales.
1 More precisely, the specific total eddy energy, with units of m 2 s −2 . For changing dissipation, while the sensitivity of the emergent E with changing bottom drag coefficient r in REF is consistent with the eddy calculation reported in , and the sensitivity of the emergent E is consistent with the GEOM int calculations reported in Mak et al. (2017), these display opposite sensitivity to each other. The resulting sensitivity of κ gm in the coarse resolution calculations is consistent with the decreasing E , as well as the circumpolar transport increasing, though the cause and effect is more convoluted (see discussion in Mak et al. 2017). This discrepancy indicates that the difference between changing λ and r, while broadly agreeing in other diagnostics, is more subtle in the resulting eddy energetics. This discrepancy is discussed at the end of this article.
The advantage of GEOM loc over GEOM int is the ability to provide a horizontal spatial structure through the emergent depth-integrated eddy energy field. Figure 4 shows the depth-averaged eddy energy field, together with the transport streamfunction of REF and GEOM loc for the control case, the large wind case and the large dissipation case. The eddy energy is mostly concentrated downstream of the ridge (located in the region −400 km ≤ y ≤ 400 km). The eddy energy signature extends further with increased wind. It is particularly noteworthy that the parameterized eddy energy (which strongly correlates with the emergent κ gm ) is able to capture aspects of the eddy energy signature displayed by REF. The emergent values may differ but additional tuning may be done to provide a better match of the magnitude. An interesting observation is that the eddy energy is extended too far to the east, which may be remedied by eddy energy propagation westward at the long Rossby phase speed, a feature not included here. However, it is noted that eddies are observed to propagate eastward at the long Rossby wave phase speed within the core of the Antarctic Circumpolar Current (Klocker and Marshall 2014), so the effect of including the additional propagation is unclear.
For completeness, the eddy energy field at the large dissipation value is also included. At larger r, the dominant contribution of the eddy energy in REF comes from the EPE. On the other hand, increasing λ in GEOM loc appears to instead concentrate the eddy energy around the ridge, with an increase in the magnitude over the ridge. The signature pattern is not entirely different to the control case and in fact resembles well the general EKE pattern of REF (not shown).

Sector configuration
To compare the characteristics of GEOM int and GEOM loc in a more complex setting, a sector configuration with a re-entrant channel connected to an ocean basin is employed. Besides a circumpolar transport, this configuration allows the possibility of a latitudinally extended residual meridional overturning circulation (RMOC). A growing number of analyses and results from eddy permitting numerical models suggests that while the circumpolar transport is largely insensitive to changes in wind forcing, the RMOC shows some sensitivity to changes in wind forcing (e.g., Hogg et al. 2008;Farneti and Gent 2011;Gent and Danabasoglu 2011;Meredith et al. 2012;Morrison and Hogg 2013;Munday et al. 2013;Farneti et al. 2015), i.e., numerical ocean models are expected to be largely eddy saturated, and partially eddy compensated. A sector configuration allows for study of whether the GEOM int and GEOM loc have the potential in reproducing both eddy saturation and eddy compensation effects, in a more complex and realistic setting.

a. Setup
The sector configuration detailed in Munday et al. (2013) was employed. As a brief summary, the domain spans 60 • S to 60 • N in latitude, with a re-entrant channel from 60 • S to 40 • S, connected to a narrow basin of 20 • in longitude. The model employs the Jackett and McDougall (1995) nonlinear equation of state. The depth is 5000 m everywhere except for a 1 • wide ridge of height 2500 m located on the eastern side of the channel (or 1 grid box for the 2 • coarse resolution calculations), which blocks the f /H contours (see Fig. 1 of Munday et al. 2013). An idealized wind forcing centered just north of the channel of the form τ s = τ 0 sin 2 (π(y + 60)/30), if y < −30, 0, otherwise, is imposed, where τ 0 is the peak wind stress and y is the latitude in degrees. On the surface, the temperature and salinity is restored to 2 T = T S + ∆T sin(π(y + 60)/120), y ≥ 0, T N + (∆T + T S − T N ) sin(π(y + 60)/120), y < 0,  Munday et al. (2013).
In this instance, the eddy permitting reference calculation (REF) has a 1/6 • horizontal grid spacing. The control simulation is taken to have control peak wind stress τ 0 = τ c = 0.2 N m −2 and control diapycnal diffusivity of κ d = κ d,c = 3 × 10 −5 m 2 s −1 . The perturbation experiments at varying τ 0 and now κ d (rather than r in the channel calculations) were restarted from perturbed states reported in Munday et al. (2013) for a further 10 years for extra diagnostics (again, a model year is 360 days). The eddy permitting reference calculations employ a biharmonic dissipation in the momentum equation that maintains a grid scale Reynolds number to be 0.15. A spatially and temporally constant GM eddy transfer coefficient of κ gm = κ 0 = 0.26 m 2 s −1 is employed.
For the coarse resolution models, the horizontal spacing is 2 • , with the CONST, GEOM int and GEOM loc variants considered; note that the domain of integration in GEOM int is taken over the whole domain rather than, for example, just over the circumpolar region. An initial calculation was first restarted from the 2 • simulation at the control parameter value of Munday et al. (2013), which is a CONST calculation with κ 0 = 1000 m 2 s −1 , for a further 1000 model years as a CONST calculation but with κ 0 = 1500 m 2 s −1 . Then perturbation experiments were carried out for a further 1800 years, following by a time-averaging  Table 2.
Two diagnostics are computed to test first for the eddy saturation property; again, the diagnostic quantities are time-averaged unless otherwize stated. The total circumpolar transport is calculated as in (10). Similar to equation (14), a pycnocline depth diagnostic is obtained by computing the pycnocline location quantity and averaging over the region between 30 • S and 30 • N, which roughly gives an estimate of the pycnocline of the ocean, avoiding the southern and northern regions where deep mixed layers may bias the results.
b. Results Figure 5 shows the aforementioned diagnostics at varying wind and dissipation values. To summarize, for varying wind, the eddying calculation REF possesses a circumpolar transport that displays weak dependence  on the peak wind stress and may be called eddy saturated. The pycnocline location is also only weakly dependent on varying peak wind stress. For increasing diapycnal diffusivity, the circumpolar transport increases and the pycnocline depth increases (more negative z pyc ). Assuming again that the circumpolar transport is dominated by thermal wind transport and that isopycnals are essentially pinned at the outcropping regions, increase in pycnocline depth are linked directly to increased circumpolar transport via increasing the tilt of isopycnals.
With this in mind, CONST is categorically not eddy saturated, displaying large sensitivity of the circumpolar transport and pycnocline depth to changing wind forcing. On the other hand, both the circumpolar transport and pycnocline location of GEOM int and GEOM loc display relative insensitivity with changing peak wind stress, which is far more consistent with the REF case. As in the channel configuration, λ in the GEOM int and GEOM loc calculations increases the circumpolar transport and pycnocline depth, much like increasing κ d in REF.

1) VARYING WIND EXPERIMENTS
While GEOM int and GEOM loc are eddy saturated, the associated sensitivity in the RMOC remains to be investigated. The RMOC may be diagnosed via the MITgcm layers package (Abernathey et al. 2011). The RMOC streamfunction is diagnosed as with x is the longitude, h = (∂ ρ/∂ z) −1 is the thickness, and the time filter is carried out in density co-ordinates. For this sector model with a nonlinear equation of state, the calculations are carried out in potential density coordinates. The potential density ρ is referenced to the 30 th model level (at around 2000 m depth), and the binning is over 241 discrete layers between 1031 and 1037 kg m −3 , equally spaced at 0.025 kg m −3 . For a simulation with an explicit representation of the eddy field, the RMOC streamfunction encapsulates both contributions from the Eulerian overturning circulation and eddy-induced transport. For all the simulations, since there is a parameterized component of the circulation via parameterized eddy induced velocity, the additional component needs to be added in (this is very weak in REF since κ 0 was chosen to be very small). The diagnosed RMOCs for varying wind forcing are shown in Figure 6. Focusing first on the control case for REF (Figure 6b; cf. Figure 8c of Munday et al. 2013), it may be seen that the RMOC consists of two main cells: (i) an upper positive cell that is the model analogue of the North Atlantic Deep Water (NADW), downwelling in the northern hemisphere and upwelling in the model Southern Ocean region; (ii) a lower negative cell that is the model analogue of the Antarctic Bottom Water (AABW), controlled largely by the convective activity occurring in the southern edges of the domain. Additionally, there is an Antarctic Intermediate Water (AAIW) negative cell slightly north of the NADW upwelling region, characterized by shallow convection.
Excursions above the surface potential density contour represents significant eddy density transport giving rise to locally higher densities in time/space. For the control wind forcing, the global morphology of the RMOC appears to be well captured in all the coarse resolution calculations, as seen in Figure 6(e, h, k) for GEOM int , GEOM loc and CONST respectively. The main differences arise in the excursion of the RMOC above the time-zonal-mean surface density in the north and in the details of the AABW negative cell. The former is because there are no explicit mesoscale eddies in the coarse resolution calculations. The latter is likely much more subtle since this involves convective processes responsible for the formation of AABW, as well as the vertical response of the eddy field via supplying the warm, salty NADW water to be transformed into AABW, and in setting the extent of the AABW cell via the eddy induced circulation.
When varying wind forcing, the changes in the RMOC displayed by REF are largely matched by GEOM int and GEOM loc . At no wind forcing, the NADW positive cell is approximately of the same magnitude and with similar extents into the southern hemisphere. At large wind, the increase in magnitude and extent in both the NADW positive cell and AABW negative cell are seen. Both GEOM int and GEOM loc struggle to reproduce the latitudinal extent and the strength of the AABW negative cell. However, both GEOM int and GEOM loc certainly appear to provide improvement on CONST, where the latitudinal extent of the NADW at zero wind forcing differs significantly from REF, and increased noise in the AABW cell and a NADW cell spanning over a smaller set of water mass classes at large wind forcing. This enhanced level of noise in and just north of the channel region in CONST coincides with increased convective activity in the same regions, where the prescribed κ gm = κ 0 is overwhelmed by the strong Eulerian overturning cell, leading to steep isopycnals and increased convective activity that is absent in REF.

2) VARYING DISSIPATION EXPERIMENTS
Increasing diapycnal diffusivity κ d increases the rate of water mass transformation, which deepens the pycnocline, thus leading to a larger region with thermal wind transport, consistent with the diagnoses of results shown in Figure 5(b, d). While increasing λ can reproduce sensitivities in the circumpolar transport, this is not the case in the resulting RMOC streamfunction, displayed in Figure 7 at the large dissipation scenario (10 × κ d,c for REF, 1.5 × λ c for GEOM int and GEOM loc ). At such a large κ d for REF, the increased rate of water mass transformation results in a RMOC with a NADW positive cell that is latitudinally confined, since the water mass is transformed within the basin before it can upwell in the re-entrant channel. On the other hand, a latitudinally extended RMOC is still somewhat maintained at large λ for GEOM int and GEOM loc . The appearance of noise in the AABW cell may be attributed to the fact that the value of the emergent eddy energy and thus κ gm has decreased (cf. Figure 8b, d), and an imbalance in the eddy induced and Eulerian overturning leads to increased convective activity. Figure 8 shows the domain-averaged eddy energy E and domain-averaged GM eddy transfer coefficient κ gm for varying input parameters, diagnosed as in the channel configuration (now with potential density instead of temperature as the gridding field when using the layers package). In this particular instance, the diagnosed domain-averaged value of EKE and EPE for REF is roughly equal in magnitude at control peak wind stress, but EKE becomes dominant especially in the channel region at large wind stress. For GEOM int and GEOM loc , E increases with increasing wind, at a roughly linear rate, which is consistent with the prediction given in Mak et al. (2017). Note that E for REF is increasingly super-linearly (cf., Munday et al. (2013) but for the domain-averaged EKE). The diagnosed E for CONST is typically two order of magnitudes smaller and appears on the axes in this plot with linear scales. The increase in κ gm for GEOM int and GEOM loc is consistent with the increase in eddy energy. The emergent κ gm for GEOM loc is smaller since κ gm is small over the basin; locally in the channel however κ gm can be large via a large local eddy energy, and κ max takes over in the model Antarctic Circumpolar Current (cf. Figure 9c, d).

3) OTHER EMERGENT QUANTITIES
With increasing dissipation, again there is further suggestion that while increasing λ results in decreased E in GEOM int and GEOM loc , consistent with the findings of the channel configuration and the results in Mak et al. (2017), it does not capture the changes displayed by changing κ d in REF.
Finally, Figure 9 shows the emergent depth-averaged total eddy energy field and the transport streamfunction for REF and GEOM loc at the control case, the large wind case, and the large dissipation case. For REF, EKE and EPE contributions to the total eddy energy are roughly equal, with the EKE contribution going from around 20% at zero wind to 80% at the largest wind forcing, and staying around 60% for changing diapycnal diffusion. Comparing with GEOM int , observe that, like the channel setting, the pattern of the emergent parameterized eddy energy -again strongly correlating with the emergent κ gm -resembles the diagnosed total eddy energy from REF around the channel region, and also in the northern hemisphere downwelling region at the control case and large wind case. Both regions possess steep isopycnals At eddy permitting resolutions, eddies generated from the channel as well as the northern sinking region may travel into the basin that, together with the presence of waves, will contribute to the eddy energy signature seen in REF. While it is reassuring to see that GEOM loc is performing in the regions where baroclinic instability is expected to strong, it remains a theoretical and modelling challenge to represent such advective effects.
Increasing diapycnal diffusion results in a larger eddy energy signature in the basin. These changes however are not captured by GEOM loc when increasing λ , again demonstrating the discrepancy between the two dissipations.

Discussion and concluding remarks
This article has outlined and described the implementation of GEOMETRIC ("Geometry of Ocean Mesoscale Eddies and Their Rectified Impact on Climate") in a three dimensional primitive equation ocean model. The GEOMETRIC recipe utilizes the Gent-McWilliams formulation but with the eddy transfer coefficient κ gm = αE(N/M 2 ), derived through rigorous mathematical bounds (Marshall et al. 2012), and with a linear dependence on the total eddy energy. This is coupled to a parameterized budget for the depth-integrated total eddy energy budget (cf. Eden and Greatbatch 2008). Done this way, the parameterization of mesoscale eddies is still through an induced adiabatic stirring as in the Gent-McWilliams scheme, but becomes energetically constrained in the vertical and varies in the horizontal through the emergent eddy energy signature.
The coarse resolution calculations utilising variants of the GEOMETRIC parameterization presented here are able to capture the bulk model sensitivities of corresponding reference calculations with an explicit mesoscale eddy field. In particular, for varying wind forcing, the coarse resolution sector model employing GEOMETRIC is eddy saturated and, furthermore, the resulting residual meridional overturning circulation also bears remarkable resemblance to the eddy permitting reference calculation, showing potential for reproducing eddy compensation.
On the other hand, this work has highlighted several subtleties, in particular with respect to eddy energy dissipation, that need to be addressed. The following discussion will focus on details of the parameterization, but it is recognized that, for example, other model details such as bathymetry play a central role in shaping the RMOC (e.g., Hogg and Munday 2014;Ferrari et al. 2016;de Lavergne et al. 2017) and will also affect the overall model response.
While the calculations with GEOMETRIC appear to capture the bulk morphological changes of the RMOC over changing wind forcing, there are features that are at odds with the reference calculation, notably in the strength and extent of the model AABW. A candidate in improving the emergent RMOC response is to incorporate a vertically varying eddy response. While this article presents results for a vertically uniform eddy transfer coefficient (Γ(z) ≡ 1), it has long been recognized that the eddy transfer coefficient should vary in the vertical (e.g., Ferreira et al. 2005). Further, since the eddy activity is expected to be strongest near the surface, the treatment of the mesoscale parameterization scheme near the ocean surface is likely going to have a large impact on the model response (e.g., Danabasoglu et al. 2008).
A set of calculations with the structure function Γ(z) = N 2 /N 2 ref (Ferreira et al. 2005) was carried out. While the associated coarse resolution calculations following the GEOMETRIC prescription captures the sensitivity in the circumpolar transport (and, in particular, is eddy saturated in the sector model), care needs to be taken so other model aspects are also reproduced. For example, in the sector configuration, an initial set of calculations with 0 ≤ Γ(z) ≤ 1 results in a shutdown of the latitudinally extended RMOC. The reason for this is that the resulting eddy response, while surface intensified, shuts off near the interface between the channel and the basin, and thus the Eulerian overturning acts unopposed in that region, causing the basin stratification to change substantially. Sample calculations with a larger imposed κ min and/or a lower bound on Γ(z) (e.g. Γ min = 0.1 as in Danabasoglu and Marshall 2007) result in a latitudinally extended RMOC. Other choices of vertical structure are possible (e.g., Ferrari et al. 2008Ferrari et al. , 2010, which may be coupled to mixed layer schemes (e.g., Large et al. 1994) and/or slope tapering schemes (e.g., Gerdes et al. 1991), all introducing additional tuning parameters. In summary, comprehensive investigation of the RMOC response under GEOMETRIC requires careful considerations of the vertical variation of the eddy transfer coefficient, among other modelling details, and is deferred to a future study.
The theory behind GEOMETRIC addresses the slumping of density surfaces in baroclinic instability. While isopycnal slumping and eddy induced stirring (from McWilliams 1990 andRedi 1982 respectively) are often implemented together (e.g., Griffies 1998;Griffies et al. 1998), in this work κ redi was fixed to be a constant in space and time, while κ gm follows the GEOMETRIC prescription. Changing κ redi is expected to affect tracer transport and is thus of great importance in the study of the ocean's role in heat transport and carbon storage, to name a few (e.g., Pradal and Gnanadesikan 2014;Abernathey and Ferreria 2015). This is beyond the scope of this work. It is noted here that diagnoses of isopycnal mixing in numerical simulations appear to show κ redi to be varying vertically and depending linearly on the eddy energy (e.g., Abernathey et al. 2013;Abernathey and Ferreria 2015). Analogous treatment of κ redi as outlined in this article may well be appropriate.
As discussed in the text, while eddy saturation is not expected to depend to leading order on the lateral redistribution of eddy energy , other details may. In the present implementation of GEOMETRIC, eddy energy is advected by the depthmean flow only, and the emergent eddy energy signature was generally found to have a more eastward extension in the coarse resolution model than the corresponding eddy permitting calculation. Inclusion of a westward advective contribution at the long Rossby phase speed (consistent with Chelton et al. 2007, Zhai et al. 2010and Klocker and Marshall 2014 is likely a remedy for the overly eastward extension of the eddy energy signature. Taking the linear eddy energy damping rate employed here at 10 −7 s −1 (dissipation rate inferred for Southern Ocean from Zhai & Marshall, pers. comm.) and a propagation speed of 0.02 m s −1 , the extent of the energy signature is on the order of 200 km, which is approximately 2 • in longitude. So while this effect may not be so significant in the Southern Ocean, it is likely significant for western boundary currents, since the inferred dissipation rate is lower in basins (Zhai & Marshall, pers. comm.). The inclusion and investigation into representing westward propagation by mesoscale eddies is a subject of a future investigation.
Perhaps the most poorly constrained aspect of the present implementation of GEOMETRIC is the treatment of eddy energy dissipation. Dissipation of mesoscale eddy energy can be through a myriad of processes, such as bottom drag (e.g., Sen et al. 2008), lee wave radiation (e.g., Naveira Garabato et al. 2004;Nikurashin and Ferrari 2011;Melet et al. 2015), western boundary processes (Zhai et al. 2010), loss of balance (e.g., Molemaker et al. 2005), all of which vary in time, space, and magnitude. Given the overwhelming complexity and the uncertainty in representing such energy pathways, the choice of linear damping of eddy energy at a constant rate over space was chosen to represent the collective effect of the aforementioned processes in this first study of GEOMETRIC. With this choice, it was found that coarse resolution models with GEOMETRIC are able to reproduce sensitivities of the circumpolar transport and thermo/pycnocline depths of the eddy permitting reference at varying wind forcing and dissipation. On the other hand, the sensitivity of the domain-averaged eddy energy magnitude, while reasonable for varying wind forcings, is completely at odds in the varying dissipation experiments. Further investigation is required to reproduce the eddy energetic sensitivities displayed in eddy permitting reference calculations. It is anticipated here that a "book-keeping" approach, accounting for the energy pathways through explicitly represented or parameterized components of the dynamics in an ocean model, will prove to be the most fruitful approach. In the present GEOMETRIC implementation, release of potential energy is accounted for in the associated eddy energy budget, but one could envisage that a parameterization for lee wave generation from geostrophic motions (e.g., Melet et al. 2015) could serve as a spatially (and temporally) varying sink in the GEOMETRIC eddy energy budget, or the sink of energy employed in GEOMETRIC could be accounted for as a source in an energetically constrained turbulence closure scheme (e.g., Gaspar et al. 1990;Madec et al. 1998), and so forth. This proposed approach, alongside individual investigations of the individual processes leading to energy transfer between scales is well beyond the scope of this work here and is deferred to future investigations.
In closing, with the understanding that there are details that can be improved upon, the results of this work lend further support to the GEOMETRIC framework as a viable parameterization scheme that better parameterizes mesoscale eddies in coarse resolution models, such that the resulting response in the emergent mean matches more closely to models that explicitly represent mesoscale eddies. For implementation into a global circulation ocean model, the primary change required is to couple a depth-integrated eddy energy budget to the existing Gent-McWilliams module. Diagnoses of eddy energetics via observations (e.g., Zhai & Marshall, pers. comm.), idealized turbulence models (Grooms 2015(Grooms , 2017 as well as ocean relevant simulations (e.g., Stewart et al. 2015) will provide a first constraint on how to improve the representation of the advection and dissipation of eddy energy, aiding in a more accurate and useful representation of the ocean climatological response.
In terms of approach, the GEOMETRIC framework marks a shift of paradigm, from a focus on how to parameterize eddy fluxes to focusing on parameterizing the eddy energetics and the associated energy pathways.