## Abstract

A framework for parameterizing eddy potential vorticity fluxes is developed that is consistent with conservation of energy and momentum while retaining the symmetries of the original eddy flux. The framework involves rewriting the residual-mean eddy force, or equivalently the eddy potential vorticity flux, as the divergence of an eddy stress tensor. A norm of this tensor is bounded by the eddy energy, allowing the components of the stress tensor to be rewritten in terms of the eddy energy and nondimensional parameters describing the mean shape and orientation of the eddies. If a prognostic equation is solved for the eddy energy, the remaining unknowns are nondimensional and bounded in magnitude by unity. Moreover, these nondimensional geometric parameters have strong connections with classical stability theory. When applied to the Eady problem, it is shown that the new framework preserves the functional form of the Eady growth rate for linear instability. Moreover, in the limit in which Reynolds stresses are neglected, the framework reduces to a Gent and McWilliams type of eddy closure where the eddy diffusivity can be interpreted as the form proposed by Visbeck et al. Simulations of three-layer wind-driven gyres are used to diagnose the eddy shape and orientations in fully developed geostrophic turbulence. These fields are found to have large-scale structure that appears related to the structure of the mean flow. The eddy energy sets the magnitude of the eddy stress tensor and hence the eddy potential vorticity fluxes. Possible extensions of the framework to ensure potential vorticity is mixed on average are discussed.

## 1. Introduction

The ocean is populated by an intense geostrophic eddy field with a dominant energy-containing scale on the order of 100 km at midlatitudes. To obtain plausible turbulent cascades of dynamic and passive tracers, it is necessary for models to resolve spatial scales of at least an order of magnitude higher, with grid spacings of 10 km if not finer (Hecht and Smith 2008). Thus, ocean climate models are unlikely routinely to resolve geostrophic eddies for the foreseeable future and further development and validation of improved parameterizations remains an important task. Moreover, development and validation of improved eddy parameterizations is an excellent strategy for testing and further advancing our knowledge of how geostrophic ocean eddies impact the large-scale circulation.

A popular parameterization of geostrophic eddies is the scheme of Gent and McWilliams (1990). This scheme is motivated by the physical principle that eddies flatten neutral density surfaces adiabatically and do so via extraction of potential energy from the mean state, modeled by the inclusion of an eddy-induced overturning circulation acting to flatten neutral surfaces (Gent et al. 1995), equivalent to a diffusion of the height of neutral surfaces. Despite many notable model improvements offered by the Gent and McWilliams closure (Danabasoglu et al. 1994; Gent 2011), the influence of eddy Reynolds stresses is neglected. This results in a number of processes being poorly modeled by the scheme, including freely decaying turbulence and the establishment of Fofonoff gyres (Cummins 1992; Wang and Vallis 1994), strong recirculations around topographic features (Bretherton and Haidvogel 1976), and the formation of jets as a result of up-gradient momentum fluxes (Starr 1968; McIntyre 1970) (also see Eden 2010, for a recent attempt to parameterize up-gradient momentum fluxes).

An alternative concept to the Gent and McWilliams closure is that geostrophic eddies mix potential vorticity along neutral density surfaces (Green 1970; Rhines and Young 1982b). The idea that eddies mix potential vorticity in this manner arises from the result that fluid parcels materially conserve their potential vorticity in the absence of forcing and dissipation and therefore might be expected to mix the potential vorticity along neutral surfaces like a passive tracer. This has led to the idea that the eddy flux of potential vorticity can be parameterized as being directed down the mean gradient potential vorticity gradient,

Here, *Q* is the Rossby–Ertel potential vorticity (e.g., Pedlosky 1987; Vallis 2006), **u** is the horizontal velocity, *κ* is an eddy diffusivity (normally assumed positive definite), and **∇**_{b} represents the component of the gradient operator directed along a buoyancy surface; primes and overbars represent transient and time-mean components, where the mean is defined through a suitable filtering operator, here evaluated on a buoyancy surface.

However, because potential vorticity is a fundamental dynamical tracer and because the dynamics are constrained by conservation principles, constraints upon *κ* are implicit in such a down-gradient potential vorticity closure. Specifically, although such a down-gradient closure leads to a consistent eddy enstrophy budget, it does not in general lead to a consistent eddy energy budget or satisfy relevant momentum constraints. Exclusion of these constraints can therefore result in physically inconsistent flows. In the following, we discuss these shortcomings by considering two simple examples, in which the eddy potential vorticity fluxes are strongly constrained by the need to conserve energy and momentum, respectively. For a more general discussion, the reader is referred to Salmon (1998, chapter 6).

### a. Conservation of energy

In Adcock and Marshall (2000), an idealized flow over topography is presented, sketched schematically in Fig. 1. The configuration consists of an inverted reduced-gravity model with a moving abyssal layer overlaying a seamount. Complete homogenization of the abyssal potential vorticity requires the layer interface between the abyssal and upper layers to rise completely over the seamount. Such motion is forbidden because of energetic constraints: such a final state has a higher energy than the initial state. Rather, Adcock and Marshall (2000) conclude, based on further numerical integrations that in this case potential vorticity is partially mixed, with the degree of mixing constrained by the initial system energy. Similar behavior was first described in a barotropic model by Bretherton and Haidvogel (1976).

In Eden and Greatbatch (2008), a modification to the Gent and McWilliams closure is formulated with such an energy constraint integrated into the formulation. The thickness diffusivity is estimated from the eddy energy using mixing length theory, and the eddy energy is carried as a prognostic variable by solving a parameterized eddy kinetic energy budget. In Marshall and Adcroft (2010), it is shown that an energetically constrained down-gradient potential vorticity closure leads to a parameterized analog of Arnold’s first stability theorem (Arnold 1965).

### b. Conservation of momentum

Consider the specific case of a zonal, periodic channel of uniform depth, as sketched in Fig. 2. Utilizing an appropriate *δ*-sheet definition of potential vorticity to account for surface buoyancy anomalies, it can be shown that, in a quasigeostrophic framework (Bretherton 1966),

Here, *q* is the quasigeostrophic potential vorticity; *υ* is the meridional velocity component; *x*, *y*, and *z* are zonal, meridional, and vertical Cartesian coordinates; and the average is evaluated at constant *z*, consistent with the quasigeostrophic approximation. This result is clear when one considers the meridional potential vorticity flux as the divergence of the Eliassen–Palm flux (Andrews and McIntyre 1976), subject to zero normal eddy flux on domain boundaries. Hence, eddies can only redistribute but not create momentum.^{1}

Now suppose that the eddy flux of potential vorticity is parameterized through a down-gradient closure. Then,

unless very strong constraints are imposed on *κ* (Green 1970; Marshall 1981; White and Green 1982; Killworth 1997). In particular, *κ* > 0 implies that must change sign within the domain; or, conversely, a mean flow for which is positive or negative definite requires that *κ* = 0. This is no more than a heuristic derivation of the Charney–Stern stability criterion (Charney and Stern 1962), once again indicating that the constraints imposed by conservation principles are closely related to stability criteria in classic stability theory (for related discussions, see Vallis 2006, chapter 6; Wood and McIntyre 2009).

### c. Aims of this paper

In this paper, we formulate a parameterization framework whereby momentum and energy constraints are satisfied by construction. Specifically, we seek to relate the eddy stress tensor, whose divergence is the eddy potential vorticity flux, to the total (kinetic plus potential) eddy energy. Subject to appropriate boundary conditions on the components of the eddy stress tensor, this satisfies the integral momentum constraint by construction and can satisfy energy conservation subject to solving for a parameterized eddy energy budget. This manuscript represents an initial contribution toward the longer-term goal of developing an eddy parameterization that mixes potential vorticity while conserving energy and momentum.

The paper proceeds as follows. In section 2, we present an overview of the new framework and list the benefits of parameterizing the eddy shape and orientation, rather than the eddy potential vorticity fluxes directly. In section 3, we derive a bound for a norm of the eddy stress tensor expressed in terms of the total eddy energy. In section 4, this bound is used to reexpress the stress tensor in terms of the eddy energy and dimensionless parameters relating to the eddy shape and orientation. In section 5, we show that the functional form of the Eady growth rate is retained by this new framework. We further show that a down-gradient buoyancy closure in the new framework yields a variant of the Gent and McWilliams scheme that is similar in form to that proposed by Visbeck et al. (1997). In section 6, we present results from a three-layer baroclinic quasigeostrophic double-gyre simulation and show that the nondimensional parameters are of tractable order (i.e., that they are not vanishingly small) and have large-scale structure. This confirms that the eddy energy sets the magnitude of the eddy stress tensor. The paper concludes with a summary in section 7, including a discussion of how one might enforce potential vorticity mixing within this framework.

## 2. The new framework

We seek a framework for parameterizing eddy potential vorticity fluxes in a manner that is consistent with conservation of both energy and momentum. We initially discuss the more general residual-mean equations, where the eddy forcing appears as the divergence of a stress tensor. We then restrict our analysis to quasigeostrophic dynamics, which, although requiring further assumptions for its applicability, provides a tractable framework in which to progress.

### a. Residual-mean equations

Our starting point is the equations of motion written in residual-mean form. The key idea, following Greatbatch and Lamb (1990), Wardle and Marshall (2000), and Ferreira and Marshall (2006), is to write the equations in terms of residual-mean variables, such that all of the eddy forcing appears on the right-hand side of the momentum equation through a term involving only the eddy potential vorticity flux.

Recently, Young (2012) has shown that the residual-mean formulation requires far fewer assumptions than previously suspected (also see Andrews 1983, for the zonally averaged case). Assuming only that the buoyancy increases monotonically with height, Young shows that the Boussinesq primitive equations can be rewritten in residual-mean form,

Here, **u** and *w* are the horizontal and vertical velocities, **∇** is the horizontal gradient operator, **∇**_{3} · is the three-dimensional divergence operator, *D*/*Dt* = ∂/∂*t* + **u** · **∇** + *w*∂/∂*z* is the Lagrangian time derivative, *f* is the Coriolis parameter, *p* is pressure, *ρ*_{0} is the reference density, **F** represents explicit mechanical forces, **k** is the unit vertical vector, *b* is the buoyancy, and represents explicit buoyancy forcing. Eddy forcing is provided by the divergence of a (3 × 2) eddy stress tensor^{2 }**E** consisting of the three-dimensional Eliassen–Palm vectors in each of its columns (see section 2c). The main virtue of this formulation is that each of the variables in Eqs. (3)–(6) are defined as residual-mean quantities where the eddy forcing is confined to the right-hand side of the horizontal momentum equation and nonresidual-mean variables do not appear. The details of the averaging in (3)–(6) are complicated, and the reader should refer to Young (2012) for full details.

For the remainder of this paper, we restrict our attention to quasigeostrophic dynamics. We make the approximation that the eddy flux of Rossby–Ertel potential vorticity along neutral density surfaces is closely related to the eddy flux of quasigeostrophic potential vorticity along height surfaces (following Charney and Stern 1962; Treguier et al. 1997), such that

(see appendix), where

is the quasigeostrophic potential vorticity; *ψ* is the geostrophic streamfunction defined such that *u* = −∂*ψ*/∂*y*, *υ* = ∂*ψ*/∂*x*, and *f* = *f*_{0} + *βy*, where *f*_{0} and *β* are constants; and . In (7) and all subsequent expressions, the average is taken at constant height.

The quasigeostrophic approximation has a number of obvious disadvantages, including the neglect of finite variations in bottom topography, the importance of which has been emphasized in section 1b. Nevertheless, we hope that the relation between the eddy forcing in the residual-mean primitive equations and in the quasigeostrophic equations presents a clear road map for applying our ideas to the residual-mean primitive equations in future manuscripts.

### b. Conservation of energy

Having written the equations in residual-mean form, it is straightforward to evaluate the conversion of energy between the mean flow and eddies. For quasigeostrophic dynamics, the eddy energy equation is

where

is the eddy energy and is an operator representing eddy energy dissipation due to bottom drag (e.g., Sen et al. 2008), loss of energy to internal waves (e.g., Molemaker et al. 2005; Polzin 2008; Marshall and Naveira Garabato 2008; Nikurashin and Ferrari 2010), western boundary dissipation (e.g., Dewar and Hogg 2010; Zhai et al. 2010), and other processes. The term in brackets on the left-hand side of (9) represents internal energy fluxes: for example, due to westward Rossby propagation (Chelton et al. 2007).

Crucially, the first term on the right-hand side of (9) represents the conversion of energy from the mean flow to the eddies; a corresponding term appears in the mean energy equation, such that the total energy (mean plus eddy) is conserved aside from explicit energy sources and sinks. Notwithstanding the challenges of parameterizing the internal eddy energy fluxes and eddy energy dissipation, here we assume that (9) is solved for the total (kinetic plus potential) eddy energy and we use this to constrain the magnitude of the eddy stress tensor, as defined in the following subsection. Note that Eden and Greatbatch (2008), in contrast, solve a prognostic eddy kinetic energy equation and use this to set the value of the eddy diffusivity, but further approximations must be applied in constructing an eddy kinetic energy budget.

### c. Eddy stress tensor

Through (7), the residual-mean eddy force, or equivalently the quasigeostrophic eddy potential vorticity flux, can be expressed as the divergence of an eddy stress tensor (Plumb 1986). The stress tensor **E** has columns equal to the Eliassen–Palm flux vectors and takes the form

where

and where *f*_{0} is the mean value of the Coriolis parameter and is the buoyancy frequency. Here, *M* and *N* represent the eddy Reynolds stresses associated with lateral momentum transfer,^{3 }*P* is the eddy potential energy. Here, *R* and *S* are proportional to the components of the eddy buoyancy flux; under the quasigeostrophic approximation, these are in turn equal to the interfacial eddy form stress associated with vertical momentum transfer (Munk and Palmén 1951; Johnson and Bryden 1989; Aiki and Richards 2008). Note that *P* yields a purely divergent residual-mean force, equivalent to a rotational potential vorticity eddy flux, and hence the stress due to *P* has no influence on the resulting dynamics.

The terms involving *M* and *N*, as well as the related barotropic “E vector,” are discussed in Hoskins et al. (1983). The lower row (the buoyancy fluxes) contains exactly those terms that are parameterized in Gent and McWilliams (1990) and Greatbatch and Lamb (1990).

The surface and bottom boundary conditions for the eddy form stresses are, neglecting variations in bottom topography (consistent with the quasigeostrophic approximation), *R* = *S* = 0. On lateral boundaries, *M* = −*K* cos2*φ*_{0} and *N* = *K* sin2*φ*_{0}, where *φ*_{0} is the angle at which the boundary is oriented with respect to the *x* axis and *K* is the eddy kinetic energy,

Note that, with no-slip boundaries, which are common in ocean models, *M* = *N* = 0.

Our proposal is to parameterize the four dynamically relevant components of the eddy stress tensor, *M*, *N*, *R*, and *S*, rather than the eddy potential vorticity flux directly. Note that, although this requires parameterization of four eddy fluxes (*M*, *N*, *R*, and *S*), this is no more than in the original equations of motion (two eddy Reynolds stresses and two eddy buoyancy fluxes), though two more than in the residual-mean equations (two eddy potential vorticity fluxes).

There are number of potential advantages of this approach:

It follows from (7), (11), and the boundary conditions that appropriate momentum constraints are satisfied. For example, in a zonal channel (2) is satisfied.

It is easily ensured that the parameterized eddy potential vorticity flux preserves the “tensorial properties” and symmetries of the original eddy potential vorticity flux (cf. Popovych and Bihlo 2011, manuscript submitted to

*J. Math. Phys.*).If we neglect the Reynolds stresses,

*M*and*N*, then this formulation reduces to parameterizing the eddy form stress, as in Gent and McWilliams (1990) and Greatbatch and Lamb (1990). Thus, the framework is a natural one to use for extending Gent and McWilliams to include contributions from Reynolds stresses.The two columns of the eddy stress tensor are the Eliassen–Palm vectors, associated with the flux of wave activity, in turn related to the propagation of eddy energy (Eliassen and Palm 1961; Andrews and McIntyre 1976). These may provide useful information for parameterizing the fluxes of eddy energy.

The eddy energy provides a rigorous upper bound on a norm of the eddy stress tensor (detailed in section 3). Note that it is the total eddy energy that is bounded, not the eddy kinetic energy used by Eden and Greatbatch (2008) for which there is no conservation law.

This upper bound allows the eddy stress tensor to be rewritten, with no loss of generality, in terms of the eddy energy, two eddy anisotropy parameters, and trigonometric factors involving three eddy orientations. Assuming that the eddy energy is available from the solution of a parameterized eddy energy equation, all of the remaining unknowns are nondimensional and bounded in magnitude by unity. All dimensional freedom in the parameterization problem is therefore encapsulated in the eddy energy (section 4).

The eddy orientations have a strong connection with classical stability theory, the eddies extracting energy from the mean flow when the eddies lean against the mean shear, and conversely returning energy to the mean flow when the eddies lean with the mean shear (section 4).

As a result of preserving the symmetries of the original eddy potential vorticity flux, the framework is able to retain some results from classical stability theory, such as the functional form of the Eady growth rate (section 5) and, with a further constraint to ensure potential enstrophy is dissipated, Arnold’s first stability theorem (section 7).

Application to the Eady problem suggests a connection between the proposed framework and the Visbeck et al. (1997) form of the Gent and McWilliams eddy closure, in the limit of negligible Reynolds stresses (section 5).

As with any new approach, there are also disadvantages. Here, the most significant disadvantage is the loss of explicit potential vorticity mixing. On the other hand, for the reasons outlined in section 1, we argue that simple-minded down-gradient diffusive closures have fundamental limitations that, despite four decades of research (since Green 1970), remain far from being resolved, if indeed this is possible (see Ringler and Gent 2011, for a related discussion). In section 7, we propose a resolution to this particular issue by solving a prognostic budget for the eddy potential enstrophy, in precisely the same manner as the eddy energy, and using this budget to ensure that potential vorticity is mixed on average through explicit dissipation of eddy potential enstrophy.

## 3. Bounded norm of the eddy stress tensor

We now proceed to bound the magnitude of the eddy stress tensor (11) by the eddy energy. This will allow us, in section 4, to rewrite the eddy stress tensor in terms of the eddy energy and bounded nondimensional parameters.

First, we note that the components of the eddy stress tensor, *M* and *N*, are bounded by

where *K* is the eddy kinetic energy (Hoskins et al. 1983). To prove (13), note that the result is an equality at any instant in time, before averaging, then apply the triangle inequality to the integral in the averaging operator. Second, we note that *P* is exactly the eddy potential energy. Third, again through application of the triangle inequality, we have

The latter inequality follows by noting *P* = *E* − *K* and maximizing the resultant quadratic.

Summing each of the above results, we obtain the final result,

where *E* is the total eddy energy as defined in Eq. (10). The left-hand side of (15) represents a weighted^{4} Frobenius norm of the eddy stress tensor and bounds its magnitude. Note that the bound is expressed in terms of the total eddy energy and not simply the eddy kinetic energy used to set the magnitude of the eddy diffusivity in Eden and Greatbatch (2008). Physically, a finite eddy buoyancy flux requires both eddy kinetic energy, for finite **u**′, and eddy potential energy, for finite *b*′.

## 4. Eddy flux angles and eddy anisotropies

Without any loss of generality, the bound (15) allows us to rewrite the components of the eddy stress tensor for an arbitrary eddy field in the form

Here, we have introduced two eddy anisotropies, 0 ≤ *γ _{m}* ≤ 1 and 0 ≤

*γ*≤ 1; two eddy flux angles, 0 ≤

_{b}*φ*≤

_{m}*π*and −

*π*≤

*φ*≤

_{b}*π*; and an eddy energy partition angle, 0 ≤

*λ*≤

*π*/2. Thus, eddy fluxes are recast in terms of six parameters, five of which are independent. Five of the parameters are dimensionless and bounded, with the only remaining dimensional parameter being the eddy energy

*E*. In the limit of plane waves, it follows that

*γ*=

_{m}*γ*= 1 and

_{b}*φ*=

_{m}*φ*=

_{b}*φ*′, the angle of phase propagation.

The anisotropy parameter *γ _{m}* measures the relative magnitudes of the Reynolds stresses and eddy kinetic energy,

Similarly, the anisotropy parameter *γ _{b}* measures the relative magnitude of the eddy buoyancy fluxes to the (square root of the) kinetic and potential energies,

One can further define a “net anisotropy” parameter,

where 0 ≤ *γ _{e}* ≤ 1. This quantifies the extent to which the bound (15) constrains the magnitude of the eddy fluxes. It follows that

*γ*,

_{e}*γ*, and

_{m}*γ*are related through

_{b}However, *γ _{e}* contains no additional information to that already contained within

*γ*,

_{m}*γ*, and

_{b}*λ*.

The angles *φ _{m}* and

*φ*measure the dominant orientation of the eddy momentum and buoyancy fluxes,

_{b}Note that there is a *π* rotation symmetry in *M* and *N* as opposed to a 2*π* rotation symmetry in *R* and *S*, consistent with the symmetry of the original eddy fluxes. Thus, parameterizing these angles ensures that the components of the eddy stress tensor retain the correct rotational symmetries. The angle *λ* measures the partitioning of eddy energy between kinetic and potential forms,

For an elliptical horizontal eddy, shown in Fig. 3a, it follows that *φ _{m}* is the angle of eddy orientation and

*γ*is a measure of the eddy eccentricity in the horizontal,

_{m}where *c*_{1} and *c*_{2} are as shown in Fig. 3a. This interpretation of *φ _{m}* and

*γ*has previously been derived in Hoskins et al. (1983).

_{m}The terms *γ _{b}* and

*λ*measure the vertical shape of the eddies. The buoyancy fluxes can be written equivalently as

where *γ _{t}* and

*φ*are defined such that

_{t}and where it is easily shown that 0 ≤ *γ _{b}* ≤

*γ*≤ 1. Then, for a vertically planar elliptical eddy, shown in Fig. 3b, it follows that

_{t}*φ*is the angle of eddy tilt in the vertical and

_{t}*γ*is a measure of the eccentricity,

_{t}where *d*_{1} and *d*_{2} are as shown in Fig. 3b.

This geometric interpretation of the eddy fluxes also provides a strong physical connection with linear stability theory (Pedlosky 1987, chapter 7), illustrated schematically in Fig. 4. Eddies extract energy from the mean flow when they lean into the shear, consistent with instability; conversely, eddies return energy to the mean flow when they lean with the shear, consistent with stability. A related discussion of observed and modeled eddy Reynolds stresses decelerating and accelerating the Kuroshio is given in Waterman et al. (2011, particularly their Fig. 24); also see Greatbatch et al. (2010) for an illustration of eddy Reynolds stresses accelerating the separated Gulf Stream.

## 5. Application to the Eady model

In the next two sections, we apply the new framework to two limiting cases: (i) baroclinic instability in a channel and (ii) fully developed turbulence in a closed, wind-driven basin.

Following Eady (1949), we consider a uniform baroclinic shear in an infinite zonal channel,

The Coriolis parameter is assumed constant (i.e., *f* = *f*_{0}), and hence the potential vorticity is uniform at each level, *Q* = *f*_{0}, except at the upper and lower boundaries, where we interpret the buoyancy gradients as a *δ* sheet of potential vorticity anomalies following the approach of Bretherton (1966). Eady theory predicts that the energy of the most unstable mode on this background shear grows at a rate,

Note that this is twice the growth rate of the eddy streamfunction and other linear eddy quantities.

### a. Eady growth rate in the new parameterization framework

We first show that the functional form of the Eady growth rate is automatically retained by the new parameterization framework. The global eddy energy budget is

The terms involving *M* and *N* integrate out through the boundary conditions. Assuming only that *S* is directed down the mean buoyancy gradient, necessary for the eddy energy to grow, and

from (16), where *α* = *γ _{b}* sin

*φ*sin2

_{b}*λ*≤ 1 is an undetermined nondimensional parameter, we find

where .

Thus, without specifying the form of the eddy closure for the components of the stress tensor, other than that the eddy buoyancy flux is on average down the mean buoyancy gradient, this new approach is guaranteed to preserve the functional form of the Eady growth rate, with an upper bound that is within a factor of 0.62 of the growth rate actually obtained. This is in marked contrast to previous down-gradient eddy closures in which arbitrary, dimensional eddy diffusivities need to be specified.^{5}

### b. A down-gradient buoyancy closure in the new parameterization framework

As shown above, the eddy Reynolds stresses have no integral effect on the eddy energy growth in the Eady problem, effectively reducing to the new approach to that employed by Gent and McWilliams (1990) and subsequent extensions. We now consider a general down-gradient buoyancy closure and infer an effective eddy diffusivity. Writing

and applying a down-gradient buoyancy closure,

it follows that

where *α* = *γ _{b}* sin2

*λ*is a nondimensional constant, bounded by unity in magnitude, and is the Richardson number. This is a general form of the Gent and McWilliams coefficient in a quasigeostrophic framework.

We can, for example, relate this form for the Gent and McWilliams coefficient to Visbeck et al. (1997). Assuming equipartition of eddy energy yields , where *u*_{eddy} is an eddy velocity. Following Green (1970) and Stone (1972), the eddy velocity can be approximated as an eddy length scale divided by the Eady growth rate, , yielding the result

This is precisely the form proposed by Visbeck et al. (1997), which drops out of the new formulation as a natural limiting case. Note, however, that because of ambiguity in the definition of the eddy length scale, the coefficient that appears here can no longer be bounded.

Of the parameters in (30), Ri depends upon the mean flow, whereas *α* and *E* are properties of the eddies. Hence a down-gradient buoyancy closure in this general form requires a parameterized eddy energy budget and a parameterization for *α*. Note, however, that Visbeck et al. (1997) find in numerical experiments that the constant in (31) has little spatial dependence. It should also be noted that, in stable regions of the flow, corresponding to regions of high Richardson number and to large Eady growth time scales, the coefficient prescribed by (30) becomes the product of a small eddy energy (as the flow is stable) but a large Eady growth time scale. This potentially complicates the direct application of this functional form of the Gent and McWilliams coefficient.

## 6. Application to wind-driven gyres

The eddy stress tensor decomposition (11) is general, but it is not self-evidently applicable as a framework for eddy parameterization. In particular, should the anisotropy parameters be small or should any of the parameters exhibit finescale spatial structure, then the potential vorticity fluxes would be highly sensitive to the details of their parameterization. In such a situation, the problem of parameterizing the geometric parameters would be highly intractable and could not be expected to capture the influence of the eddies on the mean flow.

In this section, we investigate, in an idealized oceanographic context exhibiting fully developed geostrophic turbulence, whether the nondimensional parameters have a magnitude of tractable order for the purposes of parameterization (i.e., that they have an order close to unity). We further investigate whether the parameters of the flux decomposition exhibit large-scale structure with features relating to the mean flow structure. For these purposes, we use a three-layer quasigeostrophic model to simulate a baroclinic double-gyre configuration, as described in Berloff et al. (2007).

### a. Model configuration

The three-layer quasigeostrophic equations are discretized using a potential vorticity–conserving, second-order centered finite difference formulation in space and using the leapfrog method in time with a Robert–Asselin filter of strength 0.1 in time. The simulations are conducted in a square domain −*L* ≤ *x* ≤ *L*, −*L* ≤ *y* ≤ *L* with a partial-slip boundary condition applied at all lateral boundaries. The model parameters are listed in Table 1 and correspond to a Reynolds number, defined using the upper-ocean Sverdrup velocity scale *U* = *τ*_{0}/(*ρ*_{1}*H*_{1}*Lβ*) of Re = 1600 and a Munk boundary layer scale of *δ _{M}* = (

*ν*/

*β*)

^{1/3}= 17 km. An asymmetric wind forcing is applied, equivalent to a prescribed Ekman upwelling velocity in the upper layer,

where *A* = 1.1 generates an asymmetric wind stress amplitude in each gyre, *B* = 0.2 creates a moderate tilt of the zero wind stress line to the *x* axis, and *τ*_{0} is a prescribed parameter that sets the magnitude of the wind stress.

The model is integrated using a uniform grid with 513 nodes in the meridional and zonal directions, corresponding to 7.5-km resolution and with a time step size of 1800 s. The model is integrated for 20 000 days, until a statistical equilibrium is reached (determined using measurements of the flow kinetic and potential energy), and then integrated for a further 50 000 days, with the eddy fluxes (12) computed based upon a time average over this full 50 000-day window. These eddy fluxes are decomposed into the eddy energy, the eddy anisotropies, and the eddy angles following (16). Because the buoyancy is defined at model layer interfaces and the velocity is defined within layers, the buoyancy fluxes *R* and *S* are computed within each layer by averaging the buoyancy between the two adjacent interfaces. The eddy energy is diagnosed by rewriting in the form

where this latter expression is naturally computed within each layer. The eddy potential energy is diagnosed as the residual of the eddy energy and eddy kinetic energy, *P* = *E* − *K*.

Figure 5 shows the final instantaneous transport streamfunction and potential vorticity. One sees the expected characteristic double-gyre regime, with an intense eastward jet separating two gyres and with transient eddies observed throughout the domain. Partial homogenization of mean potential vorticity in the middle layer and in the central part of the domain is also observed (Holland and Rhines 1980; Rhines and Young 1982a,b).

### b. Eddy energy

The eddy energy and eddy energy partitioning angle *λ* are shown in Fig. 6. The eddy energy is primarily concentrated in the thin eastward jet, and it decreases in the lower layers, as expected. The eddy energy partitioning indicates that the energy is approximately equipartitioned between potential and kinetic parts, with a slightly larger potential energy part in the middle larger and kinetic part in the lower layer. The normalized histogram of the partitioning angle *λ* peaks at values of 0.25*π*, 0.28*π*, and 0.16*π* in the upper, middle, and lower layers, respectively. Elevated values of *λ*, indicating relatively increased potential energy, are observed in the eastward jet and in the western boundary region.

### c. Eddy anisotropies

The momentum anisotropy *γ _{m}* and the buoyancy anisotropy

*γ*are shown in Fig. 7. Significant structure is observed in the momentum anisotropy parameter

_{b}*γ*, with larger values of

_{m}*γ*observed toward the eastern boundary. A thin zonal line of low momentum anisotropy is observed near the northern and southern boundaries. The buoyancy anisotropy

_{m}*γ*also exhibits some significant structure. In particular, two clear tracks of increased buoyancy anisotropy are observed on the flanks of the eastward jet in the upper and middle layers. Somewhat decreased values of the buoyancy anisotropy are observed in the recirculation regions. The normalized histograms of

_{b}*γ*and

_{m}*γ*in each layer are shown in Fig. 8. The histogram of

_{b}*γ*exhibits a clear peak in the upper layer at

_{m}*γ*≈ 0.15. The values of

_{m}*γ*are observed to increase in the lower layers. The boundary condition for

_{m}*γ*on lateral boundaries leads to a second peak at

_{m}*γ*= 1 in all layers. The histogram of

_{m}*γ*also exhibits a clear peak in the upper layer at

_{b}*γ*≈ 0.15. The values of

_{b}*γ*are observed to decrease somewhat in the lower layers.

_{b}### d. Eddy flux angles

For completeness, the eddy flux angles *φ _{m}* and

*φ*are shown in Fig. 9. Both

_{b}*φ*and

_{m}*φ*show significant variation in the region of the eastward jet and toward the western boundary. We do not consider the detailed structure of

_{b}*φ*and

_{m}*φ*here, except to note that the structure in the flux angles is large scale and related to the structure of the eddy anisotropies. The physical processes controlling these structures will be studied in detail in a subsequent manuscript.

_{b}### e. Implications for parameterization

These simulations represent an initial step toward assessing the utility of the eddy flux decomposition (16). Hence, we have not examined correlations between the parameters, either between each other or with the mean flow. We have also limited attention to the long-term average; one might expect low-frequency variability of the eddy statistics. Crucially, however, the anisotropy parameters are not vanishingly small—rather they are on the order of ~0.1 or larger—and all geometric parameters exhibit large-scale spatial structure with identifiable features that are related to the mean flow. These are crucial requirements for the eddy flux decomposition (16) to serve as a useful framework for eddy parameterization.

## 7. Concluding remarks

In this manuscript, we have proposed a new framework for parameterizing eddy potential vorticity fluxes in ocean circulation models. The framework involves rewriting the residual-mean eddy force, or equivalently the eddy potential vorticity flux, as the divergence of an eddy stress tensor. We have shown that a norm of the eddy stress tensor is bounded by the eddy energy. This allows the elements in the eddy stress tensor to be rewritten in terms of (i) the eddy energy, for which we propose solving prognostic equation for the eddy energy (following Eden and Greatbatch 2008), and (ii) nondimensional parameters that are bounded in magnitude by unity and have geometric interpretation in terms of eddy shapes and orientations. Numerical calculations suggests that the eddy anisotropies are typically on the order of 0.1 in fully developed geostrophic turbulence, indicating that the bound provides a useful constraint on the magnitude of the eddy fluxes.

It is important to emphasize that the new framework does not represent an eddy closure. We have merely moved the need for parameterization from one set of quantities (eddy potential vorticity fluxes) to another (eddy anisotropies and tilts). Nevertheless, assuming that the eddy energy is available at each time level, the fact that the outstanding quantities requiring parameterization are all nondimensional and bounded exerts strong constraints on the resultant eddy potential vorticity fluxes, ensuring that relevant energy and momentum constraints are preserved. We have also demonstrated strong connections between the new framework and classical stability theory.

### a. Potential vorticity mixing

The most important physics missing from the framework is the principle that eddies mix potential vorticity along neutral density surfaces. Thus, we need to ensure that any parameterization of the eddy anisotropies and tilts is consistent with potential vorticity mixing and does not permit potential vorticity unmixing. An obvious approach, analogous to the method for conserving total energy (Eden and Greatbatch 2008), is to solve a prognostic equation for the eddy potential enstrophy,

where

is the eddy potential enstrophy and is an suitable dissipation operator. Because potential vorticity mixing is equivalent to dissipation of eddy potential enstrophy, this allows eddy potential vorticity to be mixed in a controlled manner, provided that we have some mechanism to turn off the eddy fluxes when the potential enstrophy vanishes. Moreover, it is easily shown, following the same approach discussed in Marshall and Adcroft (2010), that the eddy framework then preserves an analog of Arnold’s first stability theorem (Arnold 1965): when *dq*/*dψ* is everywhere positive in free, steady flow, eddy energy can only grow at the expense of eddy potential enstrophy and the flow is stable in the sense of Lyapunov.

Finally, on a more speculative note, in Fig. 10 we show the eddy potential enstrophy and its ratio with the eddy energy. This ratio is equivalent to mean squared wavenumber in eddy energy space. In layer 1, there is close, though not exact, correspondence between the structure of this ratio and that of the eddy anisotropies plotted in Fig. 7, and this ratio, although a similar correspondence is not found in layers 2 and 3, in particular in the pool of nearly homogenous potential vorticity where the potential enstrophy is small. In future work, we will investigate the extent to which the ratio of the eddy potential enstrophy and eddy energy might constrain the eddy anisotropies or otherwise.

### b. Parameterization of eddy tilts

The numerical diagnostics suggest that equipartitioning of eddy energy provides an adequate zero-order parameterization of the vertical eddy tilts. For the horizontal eddy tilts, there are a number of potential approaches.

The first is to use classical stability theory to select eddy angles that maximize the growth of eddy energy, on the assumption that these are the eddies that emerge at finite amplitude. The second is to select eddy angles that maximize mixing of potential vorticity, although preliminary results suggest that such an approach requires calculation of high-order derivatives and is impractical.

Alternatively, to the extent that the eddies can be modeled as linear waves, one might use ray-tracing theory (e.g., Bühler 2009) to determine how the eddies are deformed by the background shear. Intriguingly, the ray-tracing formulation is also closely related to the propagation of wave activity, in turn related to the propagation (and growth) of eddy energy. It is therefore possible that a ray-tracing approach could be used simultaneously to solve for the propagation and growth of eddy energy as well as the deformation of eddies by the mean flow.

### c. Divergent eddy potential vorticity fluxes

The new framework relates the total eddy potential vorticity flux to the divergence of the eddy stress tensor. However, it is only the divergent component of the eddy potential vorticity flux that influences the evolution of the mean flow; the rotational flux results in a divergent force on the right-hand side of the momentum Eq. (3) that projects on the pressure gradient and has no influence on the evolution of the flow (in the three-dimensional, nonquasigeostrophic residual-mean equations, the divergent and rotational fluxes must be defined through a three-dimensional decomposition; cf. Marshall and Pillar 2011).

An important issue is therefore the extent to which the new framework provides useful constraints on the magnitude and structure of the divergent, rather than full, eddy potential vorticity fluxes. Previous diagnostics of rotational and divergent eddy potential vorticity fluxes suggest that the latter are typically an order of magnitude smaller (Marshall and Shutts 1981; Roberts and Marshall 2000). The new framework includes an explicit rotational potential vorticity flux through the contribution of the eddy potential energy, but it is likely that the remaining terms also project significantly onto the rotational eddy potential vorticity flux.

### d. Application to the primitive equations

Finally, although we have focused on quasigeostrophic eddy potential vorticity fluxes in this manuscript, for reasons of analytical tractability, the fact that eddy fluxes of Rossby–Ertel potential vorticity appear on the right-hand side of the residual-mean Eq. (3) provides a clear road map for applying any parameterization developed using the new framework directly to the primitive equations. Although this will require further assumptions, the use of an eddy stress tensor means that it should be straightforward to enforce any momentum constraints. We also anticipate that the application of the new framework to the primitive equations may be of value for studying the influence of variable bottom topography on the eddy potential vorticity fluxes and their forcing of mean flows.

## Acknowledgments

We wish to thank Geoff Vallis and an anonymous reviewer for exceptionally constructive and thorough reviews that led to a much improved manuscript. We also thank Bill Young and Xiaoming Zhai for additional insightful comments. Financial support was provided by the UK Natural Environment Research Council (NE/H020454/1).

### APPENDIX

#### Quasigeostrophic Residual-Mean Equations

The time-filtered quasigeostrophic momentum and buoyancy equations can be written as

where *ζ* is the relative vorticity and subscripts (…)_{g} and (…)_{ag} indicate geostrophic and ageostrophic components.

Now, define the nondivergent eddy bolus velocity (Gent et al. 1995),

This allows the momentum and buoyancy equations to be rewritten in residual-mean form,

where *K* and *P* are the eddy kinetic and potential energies. Explicit eddy forcing appears only in the momentum equation, through the eddy potential vorticity flux. The ageostrophic velocity, (**u**_{ag}, *w*), is replaced by a residual ageostrophic velocity, (**u**_{ag} + **u***, *w* + *w**), and the ageostrophic pressure is replaced by a modified pressure, *p*_{ag} + *ρ*_{0}(*K* − *P*), neither of which affects the evolution of the geostrophic flow.

## REFERENCES

## Footnotes

^{1}

One may also consider this constraint as a conservation of angular momentum; conservation of zonal momentum in the rotating frame is equivalent to conservation of angular momentum in an inertial frame.

^{2}

Young has pointed out that **E** should be defined as two separate vectors, because **E** does not, in general, transform as a tensor (W. R. Young 2011, personal communication). However, in the quasigeostrophic limit we consider for the remainder of this paper, **E** does transform as a tensor under horizontal transformations and might therefore be termed a “quasigeostrophic tensor.”

^{3}

Note that the signs of the fluxes *M* and *N* are defined inconsistently in the literature.

^{4}

Weighting of vertical components by factors of occurs frequently in quasigeostrophic theory.