## Abstract

The global stratification and circulation, as well as their sensitivities to changes in forcing, depend crucially on the representation of the mesoscale eddy field in a numerical ocean circulation model. Here, a geometrically informed and energetically constrained parameterization framework for mesoscale eddies—termed Geometry and Energetics of Ocean Mesoscale Eddies and Their Rectified Impact on Climate (GEOMETRIC)—is proposed and implemented in three-dimensional channel and sector models. The GEOMETRIC framework closes eddy buoyancy 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 models employing GEOMETRIC display broad agreement in the sensitivity of the circumpolar transport, meridional overturning circulation, and depth-integrated eddy energy pattern to surface wind stress as compared with analogous reference calculations at eddy-permitting resolutions. Notably, eddy saturation—the insensitivity of the time-mean circumpolar transport to changes in wind forcing—is found in the coarse-resolution sector model. In contrast, differences in the sensitivity of the depth-integrated eddy energy are found in model calculations in the channel experiments that vary the eddy energy dissipation, attributed to the simple prognostic eddy energy equation employed. Further improvements to the GEOMETRIC framework require a shift in focus from how to close for eddy buoyancy fluxes to the representation of eddy energetics.

## 1. 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 modeling, especially in the ocean circulation models used for climate prediction, which generally lack explicit representation of the mesoscale eddy field. Over the past two decades, a widely adopted approach for parameterizing the missing eddy fluxes has been that due to Gent and McWilliams (1990, hereafter GM). The GM scheme parameterizes eddies through both a diffusion of tracers along neutral density 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 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 existing GM-based eddy parameterizations is the very different response of the Southern Ocean circulation to changes in surface wind stress in models employing GM as compared with models with explicit eddies. Coarse-resolution models employing existing GM-based parameterizations are generally found to be more sensitive than eddy-permitting models to changing surface winds, though models employing eddy transfer coefficients that vary in three spatial dimensions are less sensitive than those that employ an eddy transfer coefficient that is varying in two spatial dimensions or is spatially constant (e.g., Farneti et al. 2015). One observed phenomenon is that the circumpolar transport increases with the strength of the surface wind forcing in coarse-resolution models employing existing GM-based parameterizations, whereas little sensitivity is observed in the equivalent models with explicit eddies (e.g., Munday et al. 2013; Farneti et al. 2015). This lack of sensitivity is known as *eddy saturation* (Hallberg and Gnanadesikan 2001; Tansley and Marshall 2001) and was first predicted on theoretical grounds by Straub (1993). Eddy saturation is generally found in models that at least partially resolve a mesoscale eddy field (e.g., Hallberg and Gnanadesikan 2006; Hogg and Blundell 2006; Hogg et al. 2008; Farneti and Delworth 2010; Farneti et al. 2010; Morrison and Hogg 2013; Munday et al. 2013; Hogg and Munday 2014) but not in models in which eddies are parameterized by GM-based schemes (e.g., Munday et al. 2013; Farneti et al. 2015).

A further discrepancy between eddy-permitting and coarse-resolution models employing existing GM-based parameterizations 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 reduced sensitivity in eddy-permitting models 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 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. On the other hand, partial eddy saturation and eddy compensation can be obtained in a model that parameterizes eddies 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). The reader is referred to the work of Farneti et al. (2015) for a recent comprehensive comparison of global circulation ocean models at coarse resolutions and their ability to reproduce eddy saturation and eddy compensation.

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 through diagnoses of numerical simulations (e.g., Ferreira et al. 2005; Ferrari et al. 2010; Bachman and Fox-Kemper 2013; Bachman et al. 2017). On the other hand, through a mixing length argument, Eden and Greatbatch (2008) proposed an eddy transfer coefficient that is related to the eddy kinetic energy (see also Cessi 2008; Marshall and Adcroft 2010; Jansen and Held 2014). This approach requires solving for the eddy kinetic energy through a prognostic eddy energy budget.

Recently, Marshall et al. (2012) have developed a new energetically constrained eddy parameterization framework, here termed Geometry and Energetics of Ocean Mesoscale Eddies and Their Rectified Impact on Climate (GEOMETRIC). The eddy forcing in the momentum equation may be described in terms of an eddy flux tensor, whose entries may be written in terms of geometric parameters that depend on the eddy kinetic and eddy potential energy. A bound of the tensor in the quasigeostrophic limit in terms of the total (i.e., kinetic and potential) eddy energy results in an inferred GM eddy transfer coefficient that is entirely determined by the total eddy energy, the stratification, and an unknown nondimensional parameter that is bounded in magnitude by unity.

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 spindown problem, as simulated by Bachman et al. (2017), the diagnosed eddy transfer coefficient from the numerical calculations is consistent with the eddy transfer coefficient predicted by GEOMETRIC across four orders of magnitude of the eddy transfer coefficient.

When applied to a two-dimensional model of the Antarctic Circumpolar Current with a domain-integrated eddy energy budget (Mak et al. 2017), GEOMETRIC produces eddy saturation, that is, a circumpolar volume transport that is insensitive to the magnitude of surface wind stress. This is due to an interplay between the zonal momentum budget and eddy energy budget (Marshall et al. 2017), the essential ingredients of which are preserved by GEOMETRIC.

However, GEOMETRIC has thus far not been implemented and tested in a three-dimensional ocean circulation model; this is the primary aim of the present study. First, we implement GEOMETRIC in an idealized three-dimensional channel, extending the calculations of Mak et al. (2017) in a two-dimensional channel model. Cases with both integrated and spatially varying parameterized eddy energy budgets are considered, comparing the results with those of eddy-permitting calculations. Second, we implement GEOMETRIC in a sector model with a basin and Southern Ocean reentrant channel, supporting an interhemispheric meridional overturning circulation in addition to a circumpolar current. The key advantage of the channel integrations is that we can afford to compare with eddy-permitting “model truths” than in the sector integrations, the latter taking far longer to equilibrate. Moreover, the channel model displays an interesting inverse sensitivity of thermal wind circumpolar transport to wind stress, which has not been previously documented but is reproduced by the GEOMETRIC parameterization. The key advantage of the sector integrations is that we are able, for the first time, to assess the extent to which GEOMETRIC is able to capture eddy compensation.

The article proceeds as follows. Section 2 outlines the GEOMETRIC approach and discusses the associated parameterized eddy energy budget. Implementation details relating to the parameterization schemes considered in this study are given in section 3. Results from GEOMETRIC are first presented for the channel model in section 4, followed by results in the sector model in section 5. The article concludes in section 6 with a summary of key findings and discussion of outstanding implementation challenges and future research questions.

## 2. GEOMETRIC

GEOMETRIC is a framework for parameterizing mesoscale eddies that preserves the conservation laws in the unaveraged equations of motion via closures that preserve the tensorial properties and symmetries possessed by the eddy flux tensor [see also related ideas in (see also related ideas in Marshall and Adcroft 2010). GEOMETRIC was originally derived under the quasigeostrophic 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:

representation of the eddy-mean flow interaction through an eddy stress tensor that can be bounded in terms of the total eddy energy in the quasigeostrophic limit (Marshall et al. 2012) and

solution of a consistent eddy energy equation accounting for parameterized and resolved dynamical processes and their role in supplying or removing eddy energy from the relevant length scales (cf. Eden and Greatbatch 2008; Cessi 2008; Marshall and Adcroft 2010; Eden et al. 2014).

In the simple limit in which the lateral eddy Reynolds stresses are neglected and the eddy buoyancy fluxes are closed as in the GM scheme (i.e., ), GEOMETRIC reduces formally to GM but as a consequence of ingredient 1, with the eddy transfer coefficient given by

Here, *E* is the total eddy energy, is the buoyancy frequency, is the magnitude of the lateral buoyancy gradient, *b* is buoyancy, and is the horizontal gradient operator. The overbar represents a time mean (over many mesoscale eddy turnover times) at fixed height, and in the context of a coarse-resolution model, may be interpreted as the buoyancy field resolved by the numerical model (see, e.g., McDougall and McIntosh 2001; Ferreira et al. 2005). An equivalent form of (1) but with the eddy kinetic energy in place of the total eddy energy was given in Jansen et al. (2015), obtained through combining a mixing length argument along with scalings derived in Larichev and Held (1995).

Crucially, given knowledge of the total eddy energy and the stratification profile, the remaining unknown *α*, satisfying in the quasigeostrophic limit (Marshall et al. 2012), is dimensionless; there is no freedom to specify dimensional quantities such as eddy length scales or eddy diffusivities. The eddy efficiency parameter *α* depends on the geometry of eddy fluxes, as documented in Marshall et al. (2012), and in principle varies in space and time (e.g., Stewart et al. 2015; Youngs et al. 2017). Results from eddy-permitting wind-driven gyre calculations in a quasigeostrophic model (Marshall et al. 2012) and eddy-resolving nonlinear Eady spindown calculations in a primitive equation model (Bachman et al. 2017) suggest that *α* is typically . Further diagnoses to constrain the functional dependence of *α* on the ocean state is a subject for further investigation; for the present work, *α* is taken to be a prescribed constant.

The remaining challenge in implementing GEOMETRIC is then to address ingredient 2, that is, to solve for the eddy energy field. Mak et al. (2017) implemented GEOMETRIC with a domain-integrated eddy energy budget while recognizing that, to delineate different dynamical regimes in more complex numerical ocean models, the eddy energy, and thus the associated parameterized eddy energy budget, should vary spatially. Solution of a prognostic equation for the eddy kinetic energy in three dimensions has been attempted by Eden and Greatbatch (2008). It is proposed here that the depth-integrated total eddy energy is solved for, as this offers both the conceptual and logistical simplicity of working in two rather than three dimensions and also avoids division by zero in (1) when the isopycnals are flat at some depth (but not if the isopycnals are flat throughout the water column).

### Eddy energy budget

For this study, GEOMETRIC is implemented in a three-dimensional model with a total eddy energy budget that varies in space and time, in contrast to the work of Mak et al. (2017), where GEOMETRIC was implemented in a two-dimensional model with a domain-integrated eddy energy budget. For the reasons mentioned above, a budget for the depth-integrated eddy energy is considered. Rather than derive a depth-integrated eddy energy budget from first principles, which contains terms that are nontrivial to parameterize (see, e.g., Eden and Greatbatch 2008), a heuristic approach is taken. We propose to use the following parameterized eddy energy budget:

Here, is the depth-averaged flow, is the unit vector in the longitudinal direction, *c* is the intrinsic long Rossby phase speed that varies with latitude, *λ* is a linear damping coefficient for the eddy energy that could in principle be a function of space and time, and is the Laplacian diffusion coefficient of the depth-integrated eddy energy. A description of the terms in (2) governing the time evolution of the eddy energy is given below.

#### 1) Advection

In observations, the eddy energy is seen to propagate at the velocity of the depth-mean flow, Doppler shifted by the intrinsic long Rossby phase speed (Klocker and Marshall 2014). In this study the contribution from the intrinsic long Rossby phase speed is not included (i.e., ) while recognizing that this may be required in an eventual implementation of GEOMETRIC in a global circulation model where it may affect the representation of western boundary currents (Chelton et al. 2007, 2011). Previous studies indicate that the lateral redistribution of eddy energy is not required to obtain eddy saturation with GEOMETRIC (Marshall et al. 2017; Mak et al. 2017), but it may affect the detailed response.

#### 2) Source

In general the sources of eddy energy depend on multiple instability types associated with the ocean state. For the present study it is assumed that the primary source of eddy energy is associated with baroclinic instability (Charney 1948; Eady 1949). The term in (2) represents the loss of available potential energy due to the slumping of neutral density surfaces as represented by GM; see Marshall et al. (2017) and Mak et al. (2017) for further details. Sources of eddy energy from other instabilities may also be included in future parameterizations but are not considered in this present study.

#### 3) Dissipation

The dissipation of eddy energy is complicated, involving a myriad of processes. These include bottom drag (e.g., Sen et al. 2008; Klymak 2018), lee wave radiation from the seafloor (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 ingredients requires detailed investigation. Instead, a simple approach is followed here, representing eddy energy dissipation through a linear damping at a rate *λ*, recognizing that *λ* parameterizes all of the physics outlined above.

#### 4) Diffusion

A Laplacian diffusion of eddy energy is incorporated following Eden and Greatbatch (2008). There are indications that the use of a Laplacian diffusion is an appropriate model of the divergence of the mean energy flux in an *f*-plane barotropic model of turbulence (Grooms 2015).

A consequence of choosing to solve for the depth-integrated total eddy energy is that the eddy transfer coefficient is energetically constrained in the vertical integral only. Specifically, suppose the eddy transfer coefficient varies in the vertical according to

where is a prescribed dimensionless structure function, (e.g., ; Ferreira et al. 2005). Then the proposed eddy transfer coefficient is (see, e.g., Mak et al. 2017)

The parameterization scheme proposed here is thus energetically constrained in the vertical and varies horizontally. Moreover, since the eddy buoyancy fluxes are still closed according to GM, the present form of GEOMETRIC retains the desirable properties of the GM parameterization, such as the positive-definite sink of available potential energy (Gent and McWilliams 1990; Gent et al. 1995), that are instrumental in its robustness.

## 3. Parameterization implementation

To assess the performance of the proposed GEOMETRIC parameterization scheme outlined above, calculations employing models with eddies parameterized by different schemes are compared with calculations from equivalent models at eddy-permitting resolutions. The parameterization schemes employed in this study are as follows: GEOM_{loc}, the GEOMETRIC parameterization scheme outlined in the previous section, with a depth-integrated but horizontally varying eddy energy equation; GEOM_{int}, the GEOMETRIC parameterization scheme outlined in Mak et al. (2017), with a domain-integrated eddy energy equation; and the standard GM scheme with a prescribed that is constant in time (CONST). For this study, the following further simplifications are made: is taken to be vertically constant such that , the intrinsic Rossby speed *c* is set to zero, and the linear damping rate of eddy energy *λ* is assumed to be a constant in space and time. The latter assumption represents a gross simplification of the physical processes’ contribution to the eddy energy sink, but we note, as an aside, the recent work of Klymak (2018) suggesting that the eddy energy sink may be linear in the bottom eddy kinetic energy.

### a. GEOM_{loc}

The first set of experiments employ GEOMETRIC locally in latitude and longitude, as detailed in section 2, with the eddy transfer coefficient computed as in (4), coupled to the parameterized eddy energy budget in (2). The GEOM_{loc} scheme is implemented wholly within the GM/Redi package in MITgcm (Marshall et al. 1997a,b), building upon the existing implementation of Visbeck et al. (1997). First, is passed through a five-point smoother; this may not be necessary but is done to conform to the MITgcm implementation of the Visbeck et al. (1997) scheme. Then, is calculated according to (4) with the smoothed (with the value of bounded below by a small number to prevent possible division by zero). In this present study, is capped below and above by a and to maintain a minimum level of parameterized eddy activity and to prevent very large eddy-induced velocities that could potentially lead to numerical instability. The GM/Redi tensor is formed and passed through a slope-tapering/slope-clipping scheme to switch off GM when isopycnal slopes become steep near outcropping regions, where mixed layer dynamics are expected to dominate, and to prevent large eddy-induced velocities.

The eddy energy budget in (2), discretized in space by centered second-order differencing, is time stepped with a third-order Adams–Bashforth scheme (started with forward Euler and second-order Adams–Bashforth steps) with the smoothed and after capping of . In this study, , , and the eddy energy diffusion coefficient is employed, the latter two chosen empirically for numerical stability. The Gerdes et al. (1991) slope-tapering scheme is chosen with the maximum slope parameter set to be , guided by the coarse-resolution setup in Munday et al. (2013).

### b. GEOM_{int}

For the second set of experiments, GEOM_{int}, the eddy energy budget is integrated over the whole domain. With *x* as longitude and *y* as latitude, the eddy transfer coefficient is calculated as

and it is coupled to the parameterized eddy energy budget given by (Mak et al. 2017)

The eddy transfer coefficient in (5) is now a constant in space but may vary in time, and the eddy energy budget in (6) becomes an ordinary differential equation. The GEOM_{int} scheme is also implemented wholly within the GM/Redi package in MITgcm following analogous steps, except the eddy energy budget in (6) is time stepped by a backward Euler scheme, in anticipation of potential numerical instabilities associated with an explicit time-stepping scheme. The same , , and slope-tapering scheme as GEOM_{loc} are used.

### c. CONST

Finally, a control case with the standard GM scheme and a constant prescribed eddy transfer coefficient,

is considered. The 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 reference calculations at eddy-permitting resolutions (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.

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 are not the same (e.g., Abernathey et al. 2013). In all calculations presented here, the Redi diffusion coefficient is prescribed to be , and the GM eddy transfer coefficient follows the prescription of GEOM_{int}, GEOM_{loc}, or CONST as appropriate.

In the coarse-resolution calculations, the parameters *α* and *λ* in GEOM_{int}, GEOM_{loc}, and in CONST are tuned for the control calculation such that the diagnosed time-averaged circumpolar transports correspond roughly to that in REF. Two sets of perturbation experiments are carried out, one set with varying wind stress at fixed dissipation and another with varying eddy energy dissipation at fixed wind stress, with each employing the same parameters as in the control. Following Marshall et al. (2017), “varying eddy energy dissipation” here is understood to mean varying the linear bottom drag for the eddy-permitting calculations and varying the parameterized linear eddy energy dissipation rate *λ*. The set of varying bottom drag calculations in the sector configuration, however, was too computationally costly, and thus the varying eddy energy dissipation experiments in the sector consist only of those varying *λ* in GEOM_{int} and GEOM_{loc} but not in REF. In this study we have not attempted to tune the parameters such that both the diagnosed circumpolar transport and domain-integrated eddy energy levels are matched to REF.

## 4. Channel configuration

### 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), a three-dimensional idealized channel configuration on a *β* plane is considered. The configuration is essentially a shorter version of the channel configuration reported in Munday et al. (2015) and Marshall et al. (2017), with no continental barriers. The domain is 4000 km long and 2000 km wide, with a maximum depth of 3000 m. The model employs a linear equation of state with temperature only and an implicit free surface. A ridge with a height of 1500 m and width of 800 km blocks contours and allows for the topographic form stress to balance the surface wind stress (Munk and Palmén 1951); the reader is referred to Fig. 3 of Munday et al. (2015) for a schematic of the bathymetry.

An idealized zonal wind stress of the form

is imposed, where is the meridional width of the channel and is the peak wind stress. The temperature is restored to the linear profile

with , on a time scale of 10 days over the uppermost cell of thickness 10 m. No surface mixed layer scheme is employed. The vertical temperature diffusivity is , except in a region of width 150 km to the north of the model domain where the vertical temperature diffusivity is tapered to to maintain a nontrivial 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, varying in thickness from 10 m at the surface to 250 m in the abyss, and with a partial cell representation of the bathymetry. A staggered baroclinic time-stepping scheme is employed. See Munday et al. (2015) and Marshall et al. (2017) for further model details.

For the eddy-permitting reference calculations termed REF, the horizontal grid spacing is uniform at 10 km. A control simulation with control peak wind stress and control bottom drag coefficient is carried out for 400 model “years” (each of 360 days), from which perturbation experiments are carried out for a further 220 model years, with time-averaged diagnostics computed in the final 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 models with parameterized eddies, the horizontal grid spacing is 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 region with enhanced vertical temperature diffusivity. A control GEOM_{int} calculation with and (consistent with observation-constrained estimates in the Southern Ocean; Zhai and Marshall 2018, manuscript submitted to *Geophys. Res. Lett*.; see also Melet et al. 2015) is first carried out over 500 model years. Perturbation experiments in GEOM_{int}, GEOM_{loc}, and CONST are then restarted and carried out for a further 500 model years, with time averages taken in the final 200 model years. The coarse-resolution calculations employ a harmonic friction in the momentum equation that forces the gridscale Reynolds number to be 0.0075.

Sensitivity experiments are then carried out in which either the magnitude of wind stress or the eddy energy dissipation (linear bottom drag for REF and *λ* for GEOM_{int} and GEOM_{loc}) are varied. 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 REF. The diagnosed total circumpolar transport is given by

where is the length of the circumpolar channel and denotes a time average performed at fixed height. Further, to assess the degree that the transport is set by the geostrophic motions, the thermal wind transport, given by

is diagnosed, where the thermal wind velocity is given by

with obtained from the temperature via the linear equation of state assuming that . Finally, following the definition of Gnanadesikan (1999) (see also Abernathey and Cessi 2014), a thermocline depth is diagnosed by computing

This is essentially a “center of mass” calculation for the temperature anomalies relative to the seafloor, and this quantity is averaged over the northern region with enhanced vertical temperature diffusivity where the thermocline is deepest, serving as a proxy for the thermocline in the basins to the north of the Southern Ocean. A deeper thermocline is expected to correlate with increased thermal wind transport, and thus total transport.

### b. Summary of key results

The key results, examining the effects of varying wind stress and eddy energy dissipation, are presented in Fig. 1. As a summary, in this channel setup, the total transport of REF decreases with increasing wind stress and increases with increased linear bottom drag. The total transport is composed principally of transport due to thermal wind, as seen in Figs. 1c and 1d. The changes in the thermal wind transport are reflected in the resulting , 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 all three diagnostics. On the other hand, both GEOM_{int} and GEOM_{loc} capture the changes of REF in all three diagnostics with increasing wind stress. The GEOM_{int} and GEOM_{loc} simulations at these choices of *λ* show similar trends to REF in the diagnostics presented.

### c. Detailed results from perturbation experiments

#### 1) Varying wind stress experiments

First, it is interesting to note that, even in REF, the total transport decreases with increasing wind, and the transport is significant even at zero wind stress. The latter is due to the enhanced vertical temperature diffusivity near the northern boundary, 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., Hogg 2010; Munday et al. 2011). In this model, the thermocline becomes shallower with increasing wind. As a result, the geostrophic flow occupies a smaller region even though the peak geostrophic flow speed may be larger, resulting in a smaller integrated thermal wind transport. The decreased thermocline depth with increasing wind is partially due to the choice of imposing high vertical diffusivity near the northern boundary; such behavior is not observed when a fully dynamical basin sets the northern channel stratification (as in the sector configuration in the next section) or when the northern boundary temperature is relaxed to a prescribed profile [as in, e.g., Abernathey and Cessi (2014), though they employ 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 GEOM_{int} and GEOM_{loc} are able to reproduce the analogous sensitivities, particularly in the thermal wind transport and thermocline depth diagnostics. The agreement between the results with GEOM_{int} and GEOM_{loc} and those with REF is less satisfactory at lower winds where the thermocline is deeper and, correspondingly, the transport is noticeably larger; the causes of these discrepancies remain to be investigated further. In contrast, the standard CONST calculations display opposite sensitivity in the transport and thermocline depth. Figure 2 shows the zonally averaged temperature profile and zonal flow of the eddy-permitting calculation and coarse-resolution calculations. The GEOM_{int} and GEOM_{loc} calculations are able to capture the changes in the stratification displayed by REF, especially in the upper ocean, in terms of the morphology and location of the temperature contours. An examination of the absolute difference in zonally averaged zonal velocity (not shown) shows the largest discrepancies lie within the high vertical temperature diffusivity region, where the coarse-resolution calculations generally have weaker zonal mean flows.

#### 2) Varying eddy energy dissipation experiments

With increased bottom drag, the total transport of REF increases, consistent with the results of Marshall et al. (2017). The rationale is that increased eddy energy dissipation requires steeper isopycnals for the eddy energy to be replenished through baroclinic instability. This leads to increased thermal wind transport, and is consistent with the diagnostics displayed in Figs. 1b, 1d, and 1f. 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).

### d. Impact on the diagnosed eddy energy and

The eddy energy and GM eddy transfer coefficient are also diagnosed. Figure 3 shows the domain-averaged eddy energy and domain-averaged GM eddy transfer coefficient . 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 (more precisely, the specific total eddy energy with units of ) for REF is calculated by diagnosing the sum of the depth-integrated (specific) eddy kinetic energy,

and the depth-integrated (specific) eddy potential energy (see, e.g., Vallis 2006, chapter 3),

Here, , where the overbar is a time average at fixed height and prime denotes deviations from that average, and , where the overbar denotes a time average at fixed density, and so denote the deviations from the mean isopycnal height. The latter time averaging is carried out 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 coordinates, with temperature referenced to the top model level at the surface and binning over 81 discrete layers between and , equally spaced at . 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 stress 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. While there are still deviations from the time mean that could be used to form an eddy energy in CONST, this is by construction small (typically three orders of magnitude smaller than REF), and thus the diagnosed eddy energies for CONST have been omitted from the diagram.

For GEOM_{int} and GEOM_{loc}, the resulting increases with increasing wind stress, though not necessarily at the same rate as REF. The rate of increase for GEOM_{int} is slightly sublinear, as opposed to the predicted linear scaling given in Mak et al. (2017). The increase in is also slightly sublinear, consistent with the behavior of . More variation is shown in GEOM_{loc} in both the resulting and levels, though the upward trend roughly follows that of GEOM_{int}. It should be noted that while in GEOM_{loc}, locally does get applied to the emergent for the simulations at larger wind forcing (e.g., Fig. 4f). At the lower peak wind stress values, the eddy energy value from GEOM_{int} and GEOM_{loc} is much smaller than REF. The corresponding is also smaller in GEOM_{int} and GEOM_{loc}, consistent with the diagnosed circumpolar transport being larger than REF in Fig. 1a.

For changing eddy energy dissipation, while the sensitivity of the diagnosed with changing bottom drag coefficient *r* in REF is consistent with the eddy-permitting calculations reported in Marshall et al. (2017), and the sensitivity of in GEOM_{int} is consistent with the GEOM_{int} calculations reported in Mak et al. (2017), these sensitivities are opposite to each other. The resulting sensitivity of in GEOM_{int} and GEOM_{loc} is consistent with the decreasing , as well as the circumpolar transport increasing, though the cause and effect is more convoluted (see the 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.

Figure 4 shows the spatially varying depth-averaged eddy energy and , together with the transport streamfunction of REF and GEOM_{loc} for the control case, the large wind stress case, and the large-eddy energy dissipation case. Generally speaking, since the models with parameterized eddies are more diffusive, the resulting flow in GEOM_{loc} possesses weaker meanders. With increasing wind, stronger recirculating gyres extend farther downstream, consistent with previous works (e.g., Nadeau and Straub 2012; Nadeau and Ferrari 2015; Munday et al. 2015). The eddy energy is mostly concentrated downstream of the ridge and extends farther east with increased wind stress. It is particularly noteworthy that the parameterized eddy energy is able to capture aspects of the eddy energy displayed by REF.

In ocean observations, mesoscale eddies within the core of the Antarctic Circumpolar Current are observed to propagate eastward at a speed consistent with advection by the background depth-mean flow, Doppler shifted by the westward propagation at the intrinsic long Rossby phase speed (Klocker and Marshall 2014). An interesting observation here in this channel model is that the eddy energy is extended too far to the east in GEOM_{loc}. The inclusion of westward propagation by the intrinsic long Rossby wave speed may remedy this deficiency by offsetting the contribution from the background mean flow advection.

In the control and large-eddy energy dissipation case in GEOM_{loc}, the resulting signature largely follows the eddy energy signature, with a slightly poleward shift, possibly in relation to the outcropping structure, and has a peak of value of around , though the domain-averaged value is much lower (see Figs. 3b,d). For the large wind stress case, the resulting is very large (hitting the imposed cap, ), and the peak regions are significantly shifted poleward of the peak eddy energy regions.

For completeness, the eddy energy fields at the large eddy energy dissipation values are 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 eddy energy pattern is not entirely different from the control case and in fact resembles well the general EKE pattern of REF (not shown), possibly indicating that the scheme as implemented is able to better capture changes in EKE pattern.

## 5. Sector configuration

A sector with a reentrant channel connected to an ocean basin is considered to allow for the possibility of an interhemispheric 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 Delworth 2010; Farneti et al. 2010; Farneti and Gent 2011; Gent and Danabasoglu 2011; Meredith et al. 2012; Morrison and Hogg 2013; Munday et al. 2013; Farneti et al. 2015). A sector configuration allows for study of whether the GEOM_{int} and GEOM_{loc} have the potential to reproduce both eddy saturation and eddy compensation in a more complex and realistic setting.

### a. Setup and diagnostics

The sector configuration detailed in Munday et al. (2013) is employed. As a brief summary, the domain spans from to in latitude, with a reentrant channel from to , connected to a narrow basin of in longitude. The model employs the Jackett and McDougall (1995) nonlinear equation of state. The depth is 5000 m everywhere except for a wide ridge of height 2500 m located on the eastern side of the channel (or one grid box for the coarse-resolution calculations), which blocks the contours (see Fig. 1 of Munday et al. 2013). An idealized wind forcing centered just north of the channel of the form

is imposed, where is the peak wind stress and *y* is the latitude in degrees. On the surface, the temperature and salinity are restored to^{1}

and

with and , over a time scale of 10 and 30 days, respectively. No mixed layer scheme is employed. The vertical domain is discretized with 42 uneven vertical levels, varying in thickness from 10 m at the surface to 250 m in the abyss. All other details are as reported in Munday et al. (2013).

In this instance, the eddy-permitting reference calculation REF has a horizontal grid spacing. The control simulation is taken to have control peak wind stress and diapycnal diffusivity of . Calculations are restarted from the perturbed states reported in Munday et al. (2013) and are integrated for a further 10 “years” (each of 360 days) to calculate the EPE as in (15). Only the perturbation experiments with varying (at fixed ) are carried out. The eddy-permitting reference calculations employ a biharmonic dissipation in the momentum equation that maintains a gridscale Reynolds number of 0.15. A spatially and temporally constant GM eddy transfer coefficient of is employed to conform to the calculations reported in Munday et al. (2013).

For the coarse-resolution calculations, the horizontal spacing is , with the CONST, GEOM_{int}, and GEOM_{loc} eddy closures considered. An initial calculation is first restarted from the simulation at the control parameter value of Munday et al. (2013), which is a CONST calculation with , integrated for a further 1000 model years as a CONST calculation but with . Then perturbation experiments are carried out for 2000 years, with time-averaged diagnostics computed over the last 200 years of the simulations. The control linear eddy energy dissipation coefficient is chosen to be , as in the channel calculations.

Varying wind stress and varying eddy energy dissipation experiments are carried out but the latter only for the GEOM_{int} and GEOM_{loc} calculations owing to computational constraints. The relevant parameter values are documented in Table 2.

The total circumpolar transport is again diagnosed as in (10). Similar to (13), a pycnocline depth diagnostic is obtained by computing the pycnocline depth:

and averaging over the region between 30°S and 30°N.

The RMOC is diagnosed via the MITgcm layers package (Abernathey et al. 2011) as

where *x* is the longitude, *υ* is the resolved meridional velocity, is the eddy-induced meridional velocity associated with the GM scheme (and is zero if the GM scheme is not active), is the thickness, and the time average is carried out in density coordinates. For this sector model with a nonlinear equation of state, the diagnoses are carried out in potential density coordinates. The potential density *ρ* is referenced to the 30th model level (at around 2000-m depth), is the potential density value at the bottom of the domain, and the binning is over 241 discrete layers between 1031 and , equally spaced at .

### b. Response of the circumpolar transport

Figure 5 shows the diagnosed circumpolar transport and pycnocline depth for varying wind stress and eddy energy dissipation values. To summarize, for varying wind stress, the eddying calculation REF possesses a circumpolar transport that displays weak dependence on the peak wind stress and may be described as eddy saturated. The pycnocline depth is also only weakly dependent on varying peak wind stress. Assuming again that the circumpolar transport is dominated by thermal wind transport and noting that isopycnals are essentially pinned at the outcropping regions, increases in pycnocline depth are linked directly to increased circumpolar transport via increases in the tilt of isopycnals and thermal wind balance.

With this in mind, the CONST calculations are 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 depth in GEOM_{int} and GEOM_{loc} display weak sensitivity to changing wind stress, far more consistent with the REF case. It is interesting to note that at slightly lower winds, the pycnocline depth of GEOM_{int} and GEOM_{loc} increases, with a corresponding signal in the diagnosed circumpolar transport, much like the channel configuration (see Fig. 1e). Increasing *λ* in the GEOM_{int} and GEOM_{loc} calculations increases the circumpolar transport and pycnocline depth. Again, the rationale is that increased eddy energy dissipation requires steeper isopycnals for the eddy energy to be replenished through baroclinic instability, leading to larger circumpolar transport through thermal wind balance.

### c. Response of the meridional overturning circulation

For the varying wind stress experiments, while the GEOM_{int} and GEOM_{loc} calculations are eddy saturated, the associated sensitivity in the RMOC remains to be investigated. The diagnosed RMOCs for varying wind stress are shown in Fig. 6. Focusing first on the control case for REF (Fig. 6b; cf. Fig. 8c of Munday et al. 2013), it may be seen that the RMOC consists of two main cells: (i) an upper positive cell that represents the model analog of North Atlantic Deep Water (NADW) downwelling in the Northern Hemisphere, upwelling in the Southern Ocean and returning northward in surface layers; (ii) a lower negative cell that represents the model analog of Antarctic Bottom Water (AABW), established by the convective activity occurring in the southern edges of the domain, spreading northward at depth, upwelling, and returning southward. Additionally, there is an Antarctic Intermediate Water (AAIW) negative cell, located slightly north of the NADW upwelling region, characterized by shallow convection.

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 Figs. 6e, 6h, and 6k for GEOM_{int}, GEOM_{loc}, and CONST, respectively. The main differences arise in the lack of an excursion of the RMOC above the time- and 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, on the other hand, likely depend on both the eddy induced circulation and convective processes; a discussion of the latter difference is deferred to the discussion section.

When varying wind stress, the changes in the RMOC displayed by REF are largely matched by GEOM_{int} and GEOM_{loc}. With no wind forcing, the NADW positive cell is approximately of the same magnitude and with similar extents into the Southern Hemisphere. With large wind stress forcing, increases in magnitude and extent of 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 improvements over CONST; where the latitudinal extent of the NADW with zero wind forcing differs significantly from REF, there is increased noise in the AABW cell, and the NADW cell spans over a smaller set of water mass classes with large wind stress forcing. The enhanced level of noise in and just north of the channel region in CONST coincides, and is consistent, with increased convective activity in the same regions, where the prescribed is overwhelmed by the strong Eulerian overturning cell, leading to steep isopycnals and increased convective activity that is absent in REF. Of course, if the initial is higher in CONST, then the noise in the RMOC may be reduced, although control calculations will become detuned.

### d. Impact on the diagnosed eddy energy and

Figure 7 shows the domain-averaged eddy energy and domain-averaged GM eddy transfer coefficient 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 values of EKE and EPE for REF are comparable in magnitude at control peak wind stress, but EKE dominates especially in the circumpolar region at large wind stress. For GEOM_{int} and GEOM_{loc}, increases approximately linearly with increasing wind stress, consistent with the prediction given in Mak et al. (2017). The increase in for GEOM_{int} and GEOM_{loc} is consistent with the increase in eddy energy. The resulting for GEOM_{loc} is smaller since is small over the basin but can be locally large in the channel; for the large wind stress, can locally reach in the model’s circumpolar current (see Fig. 8d).

With increasing eddy energy dissipation, increasing *λ* results in decreased in GEOM_{int} and GEOM_{loc}, consistent with the findings of the channel configuration and the results in Mak et al. (2017). The equivalent experiments have not been performed for REF owing to computational constraints.

Finally, Fig. 8 shows the depth-averaged total eddy energy field and the field, with the transport streamfunction for REF and GEOM_{loc}, for the control case and the large wind stress case. For REF, the EKE and EPE contributions to the total eddy energy are comparable, with the EKE contribution going from around 20% at zero wind to 80% at the largest wind forcing. For the control case, it is observed that, as in the channel setting, the pattern of the parameterized eddy energy in GEOM_{loc} resembles the diagnosed total eddy energy from REF around the circumpolar current region and also in the Northern Hemisphere downwelling region. For both REF and GEOM_{loc}, the eddy energy is large on the western part of the channel, decreasing to the east. Again, the parameterized eddy energy in GEOM_{loc} is more extended to the east than in REF. The resulting broadly correlates with the eddy energy, which in turn correlates with regions of intense baroclinicity.

At large wind stress forcing, a recirculation region is seen north of the circumpolar current in both REF and GEOM_{loc}. The eddy energy is large in the circumpolar current, with correspondingly large , hitting the cap in some places. Notably, the REF calculation displays substantially larger eddy energy values even within the basin compared to GEOM_{loc}. 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 obtained in REF, something that is not reproduced in GEOM_{loc}.

## 6. Discussion and concluding remarks

This article has described the implementation of Geometry and Energetics of Ocean Mesoscale Eddies and Their Rectified Impact on Climate (GEOMETRIC) in a three-dimensional primitive equation ocean model. The GEOMETRIC framework utilizes the Gent–McWilliams eddy parameterization but with the eddy transfer coefficient prescribed as , derived through rigorous mathematical bounds (Marshall et al. 2012) and with a linear dependence on the total (kinetic and potential) eddy energy. The eddy transfer coefficient is coupled to a parameterized budget for the depth-integrated total eddy energy [instead of an eddy kinetic energy equation in three dimensions as in Eden and Greatbatch (2008)]. Done this way, the feedback of mesoscale eddies is still parameterized as a slumping of neutral density surfaces via an eddy-induced circulation as in the Gent–McWilliams parameterization, but the eddy transfer coefficient becomes energetically constrained in the vertical integral and varies in the horizontal. The coarse-resolution calculations utilizing 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 stress, the coarse-resolution sector model employing GEOMETRIC is eddy saturated, and, furthermore, the sensitivity of the residual meridional overturning circulation to surface wind stress forcing is able to reproduce some of the eddy compensation behavior obtained in the eddy-permitting reference calculations.

On the other hand, this study has highlighted several subtleties that need to be addressed. The following discussions will focus on details of the parameterization, but it is recognized that 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; Holmes et al. 2018) 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 stress forcing, there are features that are at odds with the reference calculation, notably in the strength and extent of the modeled AABW. A candidate in improving the RMOC is to incorporate a vertically varying eddy response. While this study presents results for a vertically uniform eddy transfer coefficient , 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 to have a large impact on the model response (e.g., Danabasoglu et al. 2008; Farneti et al. 2015).

A set of calculations with the structure function (Ferreira et al. 2005) was carried out, with , as in Ferreira et al. (2005), with the magnitude of bounded above by 1. 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 that other model aspects are also reproduced. For example, in the sector configuration, an initial set of calculations with allowed to go to zero in the ocean interior resulted in a shutdown of the latitudinally extended RMOC. The reason is that, near the interface between the channel and the basin, the ocean is strongly stratified near the surface and weakly stratified in the interior, so takes small values in the interior. This means that while the eddy response is surface intensified, the Eulerian overturning acts unopposed in the interior, resulting in substantial changes to the interior stratification.

An interhemispheric RMOC for the control simulation is recovered in sample calculations with a larger imposed and/or a lower bound on (e.g., as in Danabasoglu and Marshall 2007). Other choices of vertical structure are possible (e.g., Ferrari et al. 2008, 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. A comprehensive investigation of the RMOC response under GEOMETRIC requires careful consideration of the vertical variation of the eddy transfer coefficient, among other details, and is deferred to a future study.

While slumping of isopycnals in baroclinic instability and eddy-induced stirring along isopycnals [as parameterized by Gent and McWilliams (1990) and Redi (1982), respectively] are often implemented together (e.g., Griffies 1998; Griffies et al. 1998), in this study is fixed to be a constant in space and time while follows the GEOMETRIC prescription. Changing 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 (e.g., Pradal and Gnanadesikan 2014; Abernathey and Ferreira 2015). It is noted here that diagnoses of isopycnal mixing in numerical simulations appear to show to be varying vertically and depending linearly on the eddy energy (e.g., Abernathey et al. 2013; Abernathey and Ferreira 2015). A potentially promising approach is to relate to as suggested in Smith and Marshall (2009) but is beyond the scope of this work.

As discussed in the text, while eddy saturation is not expected to depend to leading order on the lateral redistribution of eddy energy (Mak et al. 2017), other details such as the model RMOC and western boundary currents may do so. In the present implementation of GEOMETRIC, eddy energy is advected by the depth-mean flow only, and the resulting eddy energy signature is generally found to have a more eastward extension in GEOM_{int} than the corresponding eddy-permitting calculation. While the magnitude of eddy energy diffusion will have a role in redistributing the parameterized eddy energy, an obvious question is whether inclusion of a westward advective contribution at the long Rossby phase speed (consistent with Chelton et al. 2007, 2011; Zhai et al. 2010; Klocker and Marshall 2014) can remedy the overly eastward extension of the eddy energy signature. Taking the linear eddy energy damping rate employed here at [the dissipation rate inferred from Southern Ocean observations by Zhai and Marshall (2018, manuscript submitted to *Geophys. Res. Lett*.)] and a propagation speed of , one can anticipate that Rossby propagation might displace eddy energy features westward by around in longitude at the midlatitudes. So, in practice, this effect may not be sufficient to explain the observed discrepancies obtained in the circumpolar channel in the present study; however, it is likely significant for eddy hot spots adjacent to western boundary currents (Zhai et al. 2010; Zhai and Marshall 2013). The inclusion of westward propagation by mesoscale eddies is a further subject deserving further 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; Klymak 2018), 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), and 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 is chosen to represent the collective effect of the aforementioned processes. With this choice, it is found that coarse-resolution models with GEOMETRIC are able to reproduce the broad sensitivities of the circumpolar transport and pycnocline depth obtained in the eddy-permitting reference for varying wind stress and eddy energy dissipation.

On the other hand, the sensitivity of the domain-averaged eddy energy magnitude, while reasonable in the varying wind stress experiments, is at odds in the varying dissipation experiments in the channel configurations. Further investigation is required to reproduce the eddy energetic sensitivities displayed in eddy-permitting reference calculations. Additionally, while the Reynolds stresses have been neglected such that it is only buoyancy fluxes that have been closed, the inclusion of Reynolds stresses are known to be important for shaping the mean flow of inertial jets (e.g., Hughes and Ash 2001; Li et al. 2016; Tamarin et al. 2016) and for flows over variable bottom topography (e.g., Wang and Stewart 2018). The inclusion of a closure of Reynolds stresses within the GEOMETRIC framework is not pursued here but clearly is an important topic for future investigation.

In closing, with the important caveat that there are many details that can be improved upon, the results of this study lend further support to the GEOMETRIC framework as a viable parameterization scheme to better represent mesoscale eddies in coarse-resolution models, such as reproducing more accurately the response of the large-scale ocean state with explicit eddies to changes in forcing. 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 and Marshall 2018, manuscript submitted to *Geophys. Res. Lett*.), idealized turbulence models (Grooms 2015, 2017), and ocean-relevant simulations (e.g., Stewart et al. 2015; Youngs et al. 2017) 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 underlines the need to shift the focus from how to close for eddy buoyancy fluxes, to also developing improved representations of the eddy energetics and associated eddy energy pathways.

## Acknowledgments

This work was funded by the UK Natural Environment Research Council Grant NE/L005166/1 and NE/R000999/1 and utilized the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). DRM is supported by the U.K. Natural Environment Research Council (ORCHESTRA; Grant NE/N018095/1). The authors thank the two referees for comments that greatly improved the presentation of the article. The lead author thanks Gurvan Madec for discussions relating to energetic pathways. All data used herein are archived in the Edinburgh DataShare service (at https://doi.org/10.7488/ds/2297).

## REFERENCES

*Ocean Modeling in an Eddying Regime*,

*Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 319–337, https://doi.org/10.1029/177GM19.

*Atmospheric and Oceanic Fluid Dynamics.*Cambridge University Press, 745 pp.

## Footnotes

© 2018 American Meteorological Society.

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).

^{1}

Correcting a typographical error in (3) of Munday et al. (2013).