A deformation-based parametrization of ocean mesoscale eddy reynolds stresses

Ocean mesoscale eddies strongly affect the strength and variability of large-scale ocean jets such as the Gulf Stream and Kuroshio Extension. Their spatial scales are too small to be fully resolved in many current climate models and hence their effects on the large-scale circulation need to be parametrized. Here we propose a parametrization of mesoscale eddy momentum ﬂuxes based on large-scale ﬂow deformation. The parametrization is argued to be suitable for use in eddy-permitting ocean general circulation mod- els, and is motivated by an analogy between turbulence in Newtonian ﬂuids (such as water) and laminar ﬂow in non-Newtonian ﬂuids. A primitive-equations model in an idealised double-gyre conﬁguration at eddy-resolving horizontal resolution is used to diagnose the relationship between the proposed closure and the eddy ﬂuxes resolved by the model. Favourable correlations suggest the closure could provide an appropriate deterministic parametrization of mesoscale eddies. The relationship between the closure and different representations of the Reynolds stress tensor is also described. The parametrized forcing pos- sesses the key quasi-geostrophic turbulence properties of energy conservation and enstrophy dissipation, and allows for upgradient ﬂuxes leading to the sharpening of vorticity gradients. The implementation of the closure for eddy-permitting ocean models requires only velocity derivatives and a single parameter that scales with model resolution.


Introduction
Ocean mesoscale eddies are found throughout the world ocean, and are observed to be especially vigorous in regions of strong western boundary jets such as the Gulf Stream and Kuroshio Extension, as well as throughout the Antarctic circumpolar current. The strength and variability of these ocean jets is strongly enhanced by upgradient momentum fluxes due to mesoscale eddies ( Starr, 1968 ). In order for ocean general circulation models (GCMs) to accurately simulate the mean state and variability of these jets, the effects of mesoscale eddies on the large-scale flow must be represented. Unfortunately, mesoscale eddies have spatial scales in the range 10-100 km, which is too fine to be resolved by the ocean GCMs currently used in coupled climate models and so their effects must be parametrized.
There are broadly two classes of mesoscale eddy parametrization. In GCMs with very coarse horizontal resolution, ≈100 km or coarser, mesoscale eddies are not resolved at all and their effects are parametrized using a scheme that mimics their cumulative effect on the resolved scales of motion, as done by the Gent-McWilliams scheme that represents the flattening of isopycnal surfaces resulting from baroclinic instability ( Gent and McWilliams, 1990 ). In GCMs with slightly finer horizontal resolution (e.g. 50 km), often referred to as eddy-permitting models, mesoscale eddies are partially resolved, but can behave unrealistically because the eddy scale is close to the model grid scale. In this case, the goal of an eddy parametrization is to improve the representation of the eddy variability that is already partially present in the model, and potentially to improve the mean state. Since many of the next generation of coupled climate models will make use of eddy-permitting ocean GCMs, and eddy-resolving GCMs will remain too computationally expensive to be widely used in the near future, it is important to develop mesoscale eddy parametrizations that are suitable for models with eddy-permitting spatial resolutions.
Turbulent mesoscale eddies act to transfer energy from small to large spatial scales as occurs in 2D turbulence, resulting in selforganised upgradient momentum fluxes driving large-scale ocean jets ( Charney, 1971 ). For a numerical model to resolve this be- haviour, the important upscale and downscale transfers resulting from the nonlinear (i.e. advective) terms in the equations of motion must be adequately represented. Eddy-permitting models are by definition truncated somewhere in the middle of the range of scales over which such transfers occur, which implies that the effects of parameterized eddy viscosity may be felt strongly within this range because the eddy scales are close to the model gridscale. While it may be necessary for numerical stability to assume that eddies behave diffusively at the smallest resolved scales, eddy viscosity is equivalent to downgradient flux and therefore cannot represent the mean-flow forcing effects of upgradient eddy fluxes.
A more physically consistent parametrization of eddies should account for both upgradient and downgradient momentum fluxes. It should respect the quasi-geostrophic turbulence properties that there are net upscale transfers of energy and net downscale transfers of enstrophy, with eddy-mean interaction leading to the maintenance of large-scale jets. The parametrization should also depend on the variability present in the model, as it is intended for use in eddy-permitting models, and ideally this dependence should have as few as possible adjustable parameters. Our aim here is to design a parametrization that incorporates these properties.
Approaches to eddy parametrization are varied and are often explored in the context of idealised models (e.g. Frederiksen and Davies, 1997;Eden and Greatbatch, 2008;Grooms et al., 2013;Berloff, 2015a;2015b ). Berloff (2005) showed that upgradient eddy fluxes in a quasi-geostrophic model could be modelled statistically by fitting autoregressive processes to the eddy statistics obtained from an eddy resolving run, yielding an improved jet in a low-resolution model. Jansen and Held (2014) and Jansen et al. (2015) corrected spurious dissipation of kinetic energy by diagnosing the energy lost to gridscale viscosity and re-injecting this energy at larger spatial scales, so as to preserve energy conservation by mimicking the upscale cascade of energy that is expected in 2D turbulence. Porta Mana and Zanna (2014) designed an eddy parametrization by determining a function of the coarse-grained flow in a high-resolution quasi-geostrophic model that correlated well with the eddy forcing to serve as the basis for a stochastic parametrization that depends on the resolved scales. Such varied approaches all attempt to represent the rectified effects of upgradient momentum fluxes and energy backscatter (i.e. upscale energy transfer).
The approach of Porta Mana and Zanna (2014) was based on assuming that a more general stress-deformation relation than the standard eddy viscosity may apply to the flow of turbulent fluids. Rivlin and Ericksen (1955) showed how a fluid stress tensor can depend generally on the gradients of velocity as well as acceleration and so-called higher-order accelerations. The total stress tensor can be expressed as a summed series of tensors that depend on the strain tensor, with the strain tensor itself being the first term in the series. Hence truncation of the series retaining only its first term yields a diffusive forcing in the momentum and vorticity equations, and this type of fluid is termed a Newtonian fluid ( Slemrod, 1999 ). Fluids for which further terms in the series contribute to the stress are termed non-Newtonian, and we assume here as Porta Mana and Zanna (2014) did that retaining these further terms in the series provides a way to model the turbulent stress. This approach to parametrizing turbulence is not new ( Rivlin, 1957;Crow, 1968;Lumley, 1970;Meneveau and Katz, 20 0 0 ); the novelty of our study is in the application of this approach to the quasi-geostrophic turbulence that characterises oceanic mesoscale eddies, specifically the eddy Reynolds stresses due to those eddies.
In assuming a non-Newtonian stress-strain relation to be valid for the turbulent flow, we are not assuming the ocean to be an actual non-Newtonian fluid (water is Newtonian). An example of an actual non-Newtonian flow is a fluid containing polymers, i.e. macro-molecules comprised of long chains of thousands of linked molecules ( Spurk, 1997 ). Such fluids are known to show dependence of the viscosity on the fluid shear, analogously to how the ( Smagorinsky, 1963 ) parametrization scheme represents turbulence with a viscosity coefficient that depends on the fluid deformation. The additional stress tensor terms for such fluids represent the rectified effects of the macro-molecule dynamics; for example, shear-thinning behaviour, where the fluid viscosity decreases at increased shear, can result from macro-molecules being more likely to get tangled together at lower shears. By attempting to model turbulence using similar mathematical expressions, we presume that turbulent eddies represent a "microstructure" in the fluid having some net effect on the "macro" flow that we wish to model. The heuristic justification for this approach is therefore similar to that of eddy viscosity (in which an analogy is made between random molecular motions and turbulent fluid motions). Just as laboratory experiments are used to validate stress-strain relations for actual non-Newtonian fluids, here we validate our approach using results from eddy-resolving numerical simulations of an ocean primitive equations model in an idealised configuration.
The paper is structured as follows. In Section 2 we consider possible stress tensors applicable to the ocean problem, and consider the effects of one proposed parametrization on momentum, vorticity, energy and enstrophy budgets. Section 3 then uses diagnostics from numerical simulations with a primitive equations model to evaluate the vorticity forcing by the proposed parametrization. In Section 4 we discuss how our approach relates the Reynolds stress tensor of the eddying flow. Conclusions are given in Section 5 .

Theory
Our approach is to parametrize the horizontal eddy Reynolds stress tensor directly. The momentum and vorticity forcings due to parametrized eddies are then calculated from the stress tensor as flux divergences, and the parametrization satisfies conservation constraints if the relevant flux components vanish on the boundaries (e.g. Marshall et al., 2012 ). The parametrized Reynolds stress tensor will depend on the spatial gradients of velocity, and will be concisely expressed in terms of the flow deformation and the vorticity, allowing for an intuitive description of its behaviour. We will argue that the proposed stress tensor provides a parametrization that has the desirable quasi-geostrophic turbulence properties of allowing for upgradient fluxes and enstrophy dissipation.
The structure of the section is as follows. Section 2.1 introduces the deformation tensors on which the stress tensor may depend. Section 2.2 describes the momentum and vorticity forcing obtained from the proposed parametrization, and Section 2.3 considers its effects on energy and enstrophy. Section 2.4 then discusses its behaviour in more qualitative terms.

Form of the stress tensor
Consider a general stress tensor, T , representing the effects of eddies on a two-dimensional (2D) fluid flow -that is, T is a parametrization for the horizontal eddy Reynolds stress tensor. The stress tensor divergence ∇ · T , where ∇ ≡ ( ∂ x , ∂ y ), is a vector whose components give the forcing due to parametrized eddies in the corresponding vector components of the momentum equation: the x -component of ∇ · T appears on the right-hand side (RHS) of the prognostic equation for the zonal velocity, u , and the y -component of ∇ · T appears on the RHS of the equation for the meridional velocity, v . We assume an expression for T can be found that approximates the effects of small-scale turbulence on the large-scale flow. With respect to the assumption of 2D flow, al-though quasi-geostrophic turbulence (Charney, 1971) is "quasi-2D", we focus here on the simpler barotropic case.
Following Rivlin and Ericksen (1955) , we assume T can be expressed as a linear combination of tensors whose components are functions of the spatial gradients of velocity and their time derivatives. The velocity gradient tensor is (1) whose symmetric and antisymmetric parts, S = 1 2 (∇u + ∇u T ) and W = 1 2 (∇u − ∇u T ) , are referred to as the strain and vorticity tensors, respectively. Since T is symmetric 1 it can be written in the form where a, b, c are scalar valued functions. We will use the 2D matrices to provide a convenient basis for representation of any 2D matrix (e.g. Waterman and Lilly, 2015 ). Since R is antisymmetric it does not appear in T , and the matrices R + and R × represent reflection in the x -axis and reflection in an axis oriented at 45 °to the x -axis, respectively. Using these matrices, S and W can be written as where ζ = v x − u y and σ = u x + v y are the vorticity and divergence, respectively. Stretching deformation and shearing deformation are denoted by ˜ D = u x − v y and D = u y + v x , respectively. Nonzero values of the fields ˜ D and D represent changes in the shape of a fluid parcel. For ˜ D > 0 , a parcel is stretched in the x -direction and contracted in the y -direction, and vice versa for ˜ D < 0 . The effect of D is similar, but the stretching is aligned with axes rotated 45 °from the x, y -axes. At any given point in the flow it is possible to choose a reference frame, referred to here as the deformation frame , that yields a simpler picture of the deformation: by rotating to new coordinate axes x , y at an angle γ = 1 2 tan −1 (D/ ˜ D ) with respect to the original x -axis, we obtain D = 0 and ˜ D = δ > 0 in the new frame. Here δ = ˜ D 2 + D 2 is the total deformation, and like ζ and σ it is invariant under rotation of the coordinate axes. In the deformation frame, the fluid parcel is stretched in the x direction and contracted in the y direction, and the magnitude of the stretching is δ, which gives the fractional rate of change in the parcel's aspect ratio (analogously to σ giving the fractional rate of change of its area). The x -axis in the deformation frame is referred to as the axis of dilation . Spensberger and Spengler (2014) provide further description of flow deformation, including situations of interest for atmospheric synoptic-scale flows. Note also that the Smagorinsky (1963) viscosity coefficient depends on the flow deformation.
In 2D or 3D, the Rivlin and Ericksen (1955) tensors are given by . . . 1 The Reynolds stress tensor, which T is intended to parametrize, is symmetric and hence so should be T . In general a fluid stress tensor is required to be symmetric to satisfy angular momentum conservation; see Griffies (2004 , Sec. 17 is the material derivative. Eddy viscosity is obtained by setting T = T (ν) , where where in general ν can be a fourth-order tensor (e.g. Griffies, 2004 , Sec. 17.4) but here we take it to be a constant coefficient. Eq.
(9) yields downgradient fluxes of momentum and vorticity if ν > 0, since T ( ν) leads to a diffusive forcing in the equations of motion for both momentum and vorticity. Nondiffusive or upgradient effects can be obtained by using A 2 and subsequent terms, and we wish to determine whether these tensors are useful contributions to the Reynolds stress tensor for the ocean turbulence closure problem: can the effects of unresolved or partially resolved mesoscale eddies on the larger-scale flow be usefully modelled by including the terms κ m A m for m > 1 in T ? Here the κ m are assumed, for simplicity and also guided by the results of Porta Mana and Zanna (2014) and Zanna et al. (2016) , to be constant coefficients that play a role analogous to the eddy viscosity coefficient ν (i.e. κ 1 ); the functional dependence of T on the flow and on spatial location is to be captured by the A m , not by the κ m . The key physical requirement is that inclusion of such terms should provide upgradient momentum fluxes that rectify the unrealistically weak and quiescent large-scale jets typically found in coarse-resolution and eddy-permitting ocean models. We consider here the part of A 2 that does not involve a time derivative, the term ∇ u T A 1 + A 1 ∇ u . This choice is motivated partly by simplicity (for the same reason we do not consider the A m for m > 2), and also by the desire to construct a parametrization that is a function of the instantaneous state of the flow. Including the time derivative term would introduce an explicit dependence on previous times, which could perhaps be a model for "memory" effects in the flow due to sub-gridscale eddies. There is evidence of the usefulness of this approach ( Porta Mana and Zanna, 2014;Zanna et al., 2016 ), but it is not the focus of the current study. Considering then Restricting to the 2D case and using the expressions defined in Eq.
(4) , we obtain leading to with S W − W S = 1 2 This tensor, multiplied by some constant coefficient, is a useful candidate for a parametrized Reynolds stress tensor because it depends on the vorticity and the flow deformation and, unlike S itself, does not yield only downgradient fluxes, as will be shown below. In contrast, S 2 ( Eq. (13 )) appears less useful, since the term involving I contributes no forcing of vorticity since it projects onto the pressure. The term σ S could yield upgradient or downgradient flux depending on the sign of σ as well as the sign of the constant coefficient that would multiply the term, but it will vanish for incompressible non-divergent 2D flow. Although σ = 0 for the barotropic case under consideration here, we will retain σ in our expressions in anticipation of possible future extension of the framework to include quasi-2D effects. The expression S W − W S has arisen in other approaches to modelling the Reynolds stress. Pope (1975 , hereafter P75) considered the possible forms a Reynolds stress tensor depending on mean-flow gradients could take, and described this approach as a generalised "effective viscosity" 2 . A maximum of ten candidate tensors, T m ( m = 1 , 2 , . . . , 10 ), were found, of which S is the first and S W − W S the second (see P75, Sec. 3, for the complete list). If we apply the 2D representations of S and W ( Eq. (4) ) to the P75 expressions, and ignore terms involving I because they will not affect the vorticity forcing, then the T m can be shown to be equivalent to (omitting any constant coefficients) All of the T m are functions of either S or S W − W S , and if σ = 0 then only five of the T m need be considered ( m = 1 , 2 , 6 , 7 , 8 ). Another turbulence closure in which S W − W S appears is the "nonlinear gradient model" approach (Leonard 1974; Meneveau and Katz, 20 0 0), in which a Taylor series expansion using resolved-flow gradients is used to model the turbulent eddy statistics. It will be shown in Section 4 that the Taylor series can be expressed in terms of S W − W S .

Momentum and vorticity forcing
By using a stress tensor to define the parametrization, the forcing for either the horizontal momentum or vertical vorticity (i.e. v x − u y ) equations is readily obtained. This makes the parametrization suitable for models that evolve momentum, such as ocean GCMs that integrate the primitive equations. Since the vorticity forcing carries over as a forcing in the quasi-geostrophic potential vorticity equation ( Vallis, 2006 , Sec. 5.4.2), the parametrization could also be used in idealized quasi-geostrophic models. However it would only represent the horizontal Reynolds stress due to eddies, i.e. the forcing due to horizontal eddy momentum fluxes, ignoring the eddy buoyancy fluxes.
The momentum flux divergence corresponding to Eq. (2) is the 2D vector i.e. the forcing of zonal momentum is a x + b x + c y and of meridional momentum is a y − b y + c x . Using the eddy viscosity form given by Eq. (9) , The vorticity flux associated with the above momentum forcing ( Eq. (17) ) is 2 Despite the name "effective viscosity", diffusion is not the only possible behaviour arising from a stress tensor that depends on mean-flow gradients.
where the operator ∇ ⊥ is defined as For the stress tensor of Eq. (2) this yields The forcing in the vorticity equation due to T is given by the curl of the momentum forcing, ∇ · T ( Eq. (17) ), equivalent to the divergence of Eq. (20) , where the operator has been defined, and the interchangeability of the order of curl and divergence operators in Eq. (21) is due to T being symmetric. The operators appearing in Eq. (21) , 2 and 2 ∂ xy , could be described as "generalized Laplacians" following Waterman and Lilly (2015) 3 . Eq. (21) is equivalent to the convergence of the vorticity flux, −∇ · F . The isotropic part of the stress tensor, a I , does not contribute to the vorticity forcing because it appears as a gradient, ∇a , in the momentum equations (like the pressure). For the eddy viscosity form Using Eq. (15) , we define a candidate parametrized Reynolds stress tensor where κ is a constant coefficient, analogous to ν in T ( ν) . In SI units, κ has units of m 2 . Then a = 0 , b = κζ D, c = −κζ˜ D , giving from Eq.
(17) the momentum tendencies ∂u ∂t and from Eq. (21) the vorticity tendency ∂ζ ∂t Recalling that ζ = v x − u y , ˜ D = u x − v y and D = u y + v x , the parametrized momentum and vorticity forcings ( Eqs. (24) and ( 25 )) depend on spatial gradients of velocity up to second and third order, respectively, and hence the order of differentiation in both cases is the same as that of Laplacian eddy viscosity. The above expression can be further simplified. Expanding the derivatives in Eq. (25) gives and hence for incompressible 2D flow ( σ = 0 ) the vorticity tendency corresponding to T ( κ) simplifies to ∂ζ ∂t Since Eqs. (25) and (27) involve only second-order spatial derivatives, the order of differentiation is no higher than that of basic eddy viscosity (i.e. of the forcing term ν∇ 2 ζ ). The derivatives of ζ are also multiplied by first-order derivatives of the velocity, which introduces slightly more complication, but no more so than the well-established ( Smagorinsky, 1963 ) scheme in which the eddy viscosity coefficient is proportional to δ. Since higher derivatives emphasise smaller scales of motion, there is some justification in thinking of Eq. (27) as representing the effects of a large-scale velocity field (via ˜ D and D ) acting to deform a smaller-scale vorticity field (via 2 ζ and 2 ζ xy ); Section 2.4 gives more discussion. For the variants on the above given by Eq. (16) , the forcing would be given by similar expressions with additional terms arising from the flow-dependent coefficients of S or S W − W S . The form given by Eq. (27) lends itself to the following interpretation of the forcing as a modified viscosity. Let the total parametrized vorticity tendency be ∂ζ ∂t (ν,κ ) and then consider the flow in the deformation frame (as defined in Section 2.1 ), which is permissible because it can be shown that the operator ˜ D 2 + 2 D∂ xy appearing in Eq. (27) is invariant under rotation, as is the vorticity itself. 4 Then Eq. (28) in this frame becomes ∂ζ ∂t (ν,κ ) Recalling that δ > 0 (by definition of the deformation frame), the total parametrized forcing then resembles a kind of anisotropic viscosity, with different viscosity coefficients acting in the directions parallel and perpendicular to the axis of dilation. Whether the viscosity is enhanced or reduced in either direction depends on the magnitude of the deformation, as well as the sign and magnitude of κ.

Energy and enstrophy
Since the oceanic mesoscale flow is quasi-2D, we expect energy to be cascaded to large scales, and enstrophy to small scales, on average (Charney, 1971;Vallis, 2006). These cascades result from the nonlinear advective terms in the equations of motion, and the goal of a parametrization for the Reynolds stresses is to accurately model the behaviour of these terms and implied cascades. Hence it is desirable for the parametrization to dissipate enstrophy, so as to mimic the loss of small-scale enstrophy that in reality ultimately occurs due to molecular viscosity. In models, eddy viscosity will tend to dissipate enstrophy at small scales, but it will also tend to give spurious loss of energy at small scales, and various parametrization schemes have been designed to counteract this effect (e.g. Eden and Greatbatch, 2008;Jansen and Held, 2014 ). Hence it is also desirable that the parametrization conserves total energy, so as not to contribute to this spurious loss.
A stress tensor T expressed in the form given by Eq. (2) can be shown to give a forcing in the 2D energy equation of 4 Note from Eq. (26) that the full expression for the forcing, for the case when σ = 0, is invariant under rotation since it is the curl of a vector. The ∇ 2 operator, ζ and σ are invariant, hence ζ ∇ 2 σ is invariant, as is ∇ζ · ∇σ since it is the cosine of the angle between ∇ζ and ∇σ , multiplied by | ∇ζ || ∇σ |, all of which are invariant. This invariance is of course required; since ζ is invariant under rotation, it cannot be possible to obtain a different value for ∂ ζ / ∂ t by rotating the coordinate axes. consistent with eddy viscosity being a sink of energy (for ν > 0).

For the proposed parametrization T ( κ) the energy tendency is
the non-flux energy tendency term vanishes 5 . Since the energy forcing by T ( κ) is the divergence of a flux, its domain-integrated value depends on the flux at the boundaries. If the boundary flux is zero then the parametrization will conserve energy, which is the case for both no-slip boundary conditions ( u, v = 0 on the boundary) and free-slip conditions ( ζ = 0 on the boundary) 6 . Because the parametrization conserves energy, its effect should be to redistribute the resolved eddy energy, acting against dissipation. The enstrophy tendency corresponding to a stress tensor T is given by where Z ≡ 1 2 ζ 2 is the enstrophy, and the vorticity flux and vorticity tendency corresponding to T are F and −∇ · F , respectively. Only the term F · ∇ζ contributes to the global enstrophy budget.
The vorticity flux for T ( κ) using Eq. (20) is given by where the second line follows from Eq. (4) and the identity Using the identity D Dt which follows from ∇u T = S − W and the antisymmetry of W , the enstrophy tendency from the non-flux term is then ∂Z ∂t Consider adiabatic motion, for which the quasi-geostrophic equation of motion in the 2D case is D (ζ + βy ) /Dt = 0 . For motion on scales smaller than the Rhines scale Z/β, where Z is a typical magnitude of the vorticity, the β-term becomes small, so that D ζ / Dt ≈0. With β = 2 × 10 −11 m −1 s −1 and Z = 10 −5 s −1 (about 10% of a typical midlatitude value of f ) then Z/β = 500 km, which is considerably larger than the small-scale eddies whose effects we wish to parametrize. In this case, with Dζ /Dt = 0 the non-flux contribution to the enstrophy budget, Eq. (37) , is proportional to the parcel-following rate of change of the (squared) magnitude of the vorticity gradient. If κ > 0 then the effect of T ( κ) is to locally remove enstrophy when | ∇ζ | is increasing -that is, when vorticity gradients are sharpening. Since increasing sharpness of vorticity gradients is associated with enstrophy cascading to smaller scales  ( Weiss, 1991;Vallis, 2006 ), this behaviour is consistent with the small scales providing a sink of enstrophy. Although a downscale enstrophy cascade is expected to occur at scales of motion that are effectively adiabatic (i.e. in the absence of direct large-scale forcing, as would occur in the ocean interior), any practicable numerical model must include some dissipation at small scales leading to diabatic terms. Consider a spatial domain within which enstrophy dissipation is assumed to occur: DZ/Dt < 0 , where the overbar represents the spatial average at a fixed time. We then have DZ/Dt = ζ Dζ /Dt < 0 , so that ζ and D ζ / Dt are anticorrelated over this domain, which means that ∇ζ and ∇D ζ / Dt point in opposite directions. (Since D ζ / Dt is equal to the total of all forcings acting on the vorticity, this is simply saying that where ζ > 0 the forcing is acting to reduce ζ , and vice versa.) Hence 2 κ∇ ζ · ∇ Dζ /Dt < 0 if κ > 0, making this term in Eq.
(37) an overall enstrophy sink, consistent with the initial assumption that there is small-scale dissipation that removes enstrophy.
On the other hand, κ < 0 would make this term into an overall enstrophy source, which is inconsistent with the initial assumption.

Flux behaviour
Here we consider the qualitative behaviour of the proposed parametrization T ( κ) on the vorticity and momentum budget. From Eq. (34) , ignoring the ∇ ⊥ ζ term that does not affect the vorticity forcing, assuming σ = 0 , and using Eq. (4) , the vorticity flux is If the coordinate axes are rotated to the deformation frame (as defined in Section 2.1 ) then the vorticity flux simplifies to The conservation law for a tracer whose density is A and whose flux is F A , without sources or sinks, is This states that the tendency of A due to transport by the flux F A is −∇ · F A , so that the tracer accumulates in regions where the flux F A converges (i.e., ∇ · F A < 0) and vice versa. Writing Eq. (25) in this form using the vorticity flux in the deformation frame gives ∂ζ ∂t (κ )   (i.e. we consider the interaction between different scales of motion). In Fig. 1 a, a deforming flow, δ, is schematically represented by streamlines and velocity vectors. Fig. 1 b depicts the sense of F (κ ) = κδR + ∇ζ (represented by large arrows) resulting from δ acting on a small patch of vorticity (whose ∇ζ is shown by thin arrows), for κ > 0. The vorticity flux removes positive vorticity from the regions to the left and right of the patch, and adds positive vorticity to the regions above and below the patch, which will tend to strengthen vorticity gradients in one direction and weaken them in the other, consistent with the general behaviour indicated by Eq. (29) . Hence the effect of F ( κ) is qualitatively consistent with our goal that the parametrization should be able to not only weaken vorticity gradients (as eddy viscosity already does) but also to strengthen them, as required to maintain and strengthen jets. Due to the nonlinearity of the vorticity flux, Eq. (41) represents interaction between different scales. To illustrate this using a 2D Fourier expansion of the streamfunction, ψ = k ψ k (t) exp (i k · x ) , where ψ k ( t ) is the time-dependent amplitude of the spatial mode of wavenumber k , from Eq. (27) it can be shown that the vorticity tendency for a single Fourier mode with wavenumber k + p is where k = (k, l) , p = (p, q ) are the spatial wavevectors of the ζ and δ terms in Eq. (41) , respectively. Nonzero vorticity tendency requires that k = p , i.e. ∇ · F ( κ) represents interaction between different scales of motion. Transfer of information between scales is required in order for eddies to produce upgradient fluxes that reinforce large-scale jets. It is well known that eddies tilting "with the shear", as if advected by a sheared mean flow, act to flux momentum into the core of a jet and remove it from the jet flanks ( Vallis, 2006;Dritschel and McIntyre, 2008 ). The parametrization T ( κ) requires that some deformation be present in order to give a nonzero forcing, e.g. as in eddy-permitting models. What forcing would result from T ( κ) for the case of a "banana-shaped" eddy, i.e. an eddy tilting with the shear? Fig. 2 schematically shows a tilted eddy ( Fig. 2 a) and the resulting forcing due to T ( κ) when averaged zonally ( Fig. 2 b-d).
Since the eddy structure is given, the Reynolds stress can be calculated and is seen in Fig. 2 b-d to agree qualitatively with the parametrization. In particular, Fig. 2 b shows that eastward acceleration is obtained in the centre of the eddy and westward acceleration on its flanks, as desired. This suggests that T ( κ) can produce reasonable results, such as upgradient momentum fluxes and jet sharpening, if its input contains partially resolved eddies.
The qualitative agreement in Fig. 2 is not sensitive to the precise details of the eddy shape, so long as the eddy tilts with the shear. The reason for agreement is found by plotting the x, y structure of u, v , ζ , ˜ D and D (not shown): the spatial phase relationship between u and v is similar to that between ζ and ˜ D , leading to similar correlations between them when zonally averaged. This suggests that the impact of T on a large-scale zonal jet would be qualitatively similar to that of the true Reynolds stresses, provided the resolved-flow spatial structures of ζ , ˜ D , D are similar enough to the spatial structures of the unresolved eddies. This argument is essentially the "self-similarity" hypothesis used to justify the nonlinear gradient model, mentioned in Section 2.1 and to be discussed further in Section 4 . In the next section we examine whether numerical simulations support this hypothesised similarity.

Model description
We use the MITgcm, a hydrostatic primitive equations ocean GCM ( Marshall et al., 1997 , see http://mitgcm.org for model source code, documentation, and further information.). The model setup mimics that of idealised double-gyre experiments commonly performed with quasi-geostrophic (QG) models; specifically, parameters are chosen to match as closely as possible those used in Porta Mana and Zanna (2014) , albeit with increased vertical resolution suitable for the primitive equations.
The model spatial domain extends 40 0 0 km zonally and meridionally, and 4 km deep with a flat bottom. Planetary rotation is given by Coriolis parameter f = 10 −4 s −1 on a beta plane with β = 2 × 10 −11 m −1 s −1 (and Cartesian x, y, z coordinates are used).
The spatial resolution is 7.5 km horizontally, and 44 levels vertically with vertical level thickness increasing with depth as is conventional. Explicit dissipation mechanisms are Laplacian horizontal viscosity of 100 m 2 s −1 , vertical viscosity 10 −5 m 2 s −1 , and linear bottom drag with coefficient 1 . 2 × 10 −4 m s −1 . Lateral boundary conditions are free-slip, there is no rigid lid, and sea-surface temperature is relaxed with a 360-day relaxation timescale toward a annual mean target state varying linearly from 25 °C at the southern boundary to 15 °C at the northern boundary. The model timestep is 300 s. A linear thermodynamic equation of state is used, with no salinity. Forcing is by a wind stress that is constant in time and mimics the extratropical atmospheric jet structure, with a midlatitude maximum of 0 . 3 N m −2 and a southwestto-northeast tilt.

Diagnostics
Model output was saved at 10-day intervals, and in some instances at 3-day or 1-day intervals to test the sensitivity of the diagnostics to the sampling frequency. Unless otherwise stated, the 10-day data are analysed. Integrations were begun from a state of rest and required approximately 3-5 years to spin up, as indicated by time series of total horizontal kinetic energy. Integrations were ended sometime between 10 and 20 years, and diagnostics are calculated using data from the statistical steady state.
To diagnose the eddy forcing, we use coarse-graining to calculate the difference between the eddy forcing in a high-resolution model and the eddy forcing that would be expected in a hypothetical model having lower resolution. Let the equation of motion for vorticity in the high resolution model be The subscript f on all quantities denotes the fine-grain model: ζ f is defined on the high-resolution ( x, y ) grid, spatial derivatives (∇ f , ∇ 2 f ) are calculated using this grid, and ν f is the constant viscosity coefficient. The imposed wind stress curl forcing is F f and the only parametrization of eddies is the viscous term, ν f ∇ 2 f ζ f . The Coriolis parameter is f 0 + βy where f 0 and β are constants.
The fine-grain model adequately resolves the eddy forcing -in other words, we regard the fine-grain model as representing the "true" effect of the eddy forcing on the larger-scale mean flow.
We then coarse-grain Eq. (43) by averaging over some area containing more than one of the gridboxes of the fine-grain model: The overbar denotes the coarse-graining operation, i.e. horizontal spatial averaging procedure. It does not affect the temporal or vertical discretization.
The ideal behaviour of a hypothetical lower-resolution model is that it mimics (statistically) the evolution given by ∂ ζ f /∂t. Let this lower resolution model evolve according to where the subscript c refers to the coarse-grain grid on which fields are defined. The viscosity of a lower resolution model is necessarily larger than the smallest value allowable in the fine-grain model: ν c > ν f . The additional term S ζ is a parametrization of the eddy forcing, intended to capture the missing effects of the unresolved eddies and thereby yield a more realistic mean flow. The explicit form of S ζ , as a function of the coarse-grain flow, is at this point unknown. We use the fine-grain simulation, coarsegrained according to Eq. (44) , to diagnose S ζ so as to compare it with the forcing given by any proposed parametrization that can be expressed as a function of the coarse-grain flow, such as T ( κ) .
Since we want the evolution of ζ c to mimic that of ζ f , let ∂ ζ c /∂ t = ∂ ζ f /∂t. Then combining Eqs. (44) and (45) , and then with v f = v c , and if the imposed large-scale forcing is smooth enough that F f = F c , the diagnosed eddy source function for vorticity is then The same procedure can be followed for other flow fields. The corresponding expressions for momentum (recall from Section 3.1 that u, v are the actual prognostic variables in our model) are Here additional terms involving β and the pressure gradient arise that did not appear in the corresponding expression for vorticity ( Eq. (47) ). Diagnosis of S ζ has the advantage that the eddy source term is simpler, not requiring computation of the extra terms appearing in the expressions for S u and S v , but the disadvantage that an extra order of differentiation is required (exacerbating numerical inaccuracies). Due to the considerable simplification provided by S ζ , and also motivated by the fundamental role played by vorticity in geostrophic and 2D turbulence, we will consider only the vorticity source term. (Further discussion regarding the order of differentiation is given in Section 4.2 .)

Simulations
Time-mean barotropic streamfunction (calculated from the vertically integrated horizontal velocity) is shown in Fig. 3 . The windforced double gyre structure is present, with a strong jet separating from the western boundary and extending into the ocean interior. Corresponding coarse-resolution experiments, with 30 km horizontal grid spacing (i.e. four times coarser than the model used in Fig. 3 ) and adjusted horizontal viscosity and timestep exhibit a substantially weaker time-mean jet (not shown). The coarse-resolution experiments show considerably less variability than their high-resolution counterparts: the weak jet tends not to meander, and the time-mean picture is more representative of the actual jet structure at any given instant. These problems with coarse-resolution ocean models -i.e. models too coarse to adequately resolve geostrophic turbulence, which occurs at scales of 10 to 100 km, comparable to the deformation radius -are well known, and motivate development of parametrizations for the effect of unresolved eddies on the large-scale flow.  terior ocean flow 7 . Turbulent eddies are evident throughout much of the domain, with the strongest eddies occurring near the western boundary (a nonlinear contour scale is used in Fig. 4 so that the weaker eddies elsewhere in the domain can also be visualized). A strong jet separating from the western boundary, near y ≈20 0 0 km, is readily apparent. The horizontal gradient of vorticity, ∇ζ , is sharpest at the jet, where its spatial structure includes variations on fine scales of order ≈50 km, a few times larger than the horizontal gridscale of the model (7.5 km). This fine-scale structure tends to be sharpest within about 10 0 0 km of the western boundary ( x < 10 0 0 km). The jet then becomes more diffuse as it extends further eastward (between x ≈10 0 0 km and 30 0 0 km; Fig. 3 ). Qualitatively this behaviour resembles that of the Gulf Stream or Kuroshio extension (e.g. Jan et al., 2015 ). Fig. 5 shows the diagnosed eddy source function for vorticity, S ζ , for the same day as shown in Fig. 4 . A coarse graining of n = 4 is used to calculate S ζ , where n represents the length of a side of a square coarse-grain gridbox in units of the fine-scale grid length (7.5 km; hence n = 4 corresponds to box sides of 30 km length, representing a 16-to-1 loss of information between the fine and coarse grids). From Eq. (47) , calculation of S ζ also requires choosing an appropriate value for ν c , which should correspond to a typical horizontal viscosity that would be feasible for a model with a horizontal resolution equivalent to the coarse-grain gridbox size. We use ν c = n ν f , where ν f = 100 m 2 s −1 is the viscosity in the high resolution model ( Section 3.1 ). Hence for n = 4 , corresponding to (30 km) 2 coarse-grain gridbox, ν c = 400 m 2 s −1 . However, since the advective terms in Eq. (47) are considerably larger than the viscous terms for reasonable values of ν c , the results are not sensitive to this choice.
As expected, Fig. 5 shows that S ζ is strongest in the region where the fine-scale variability of the jet is strongest, close to its separation from the western boundary. Since S ζ indicates the dif- ference between the actual eddy forcing and the forcing calculated using coarse-grained quantities, large values at the jet reflect the fact that coarse-graining smooths out the sharp gradients that give rise to a substantial portion of the advective term, u · ∇ζ . In regions where sharp gradients are less common, such as the more quiescent eastern part of the model domain, S ζ is much smaller due to the coarse-grained flow being a reasonably accurate representation of the actual flow (note the nonlinear contour scale).

Parametrization diagnostics
To test whether a proposed parametrization can represent the effects of partially resolved eddies, we evaluate its skill in predicting S ζ . Fig. 6 a shows the correlation between S ζ and the vorticity forcing that would result from a stress tensor T = S W − W S ( Eq. (15) , or equivalently, Eq. (23) divided by the constant 2 κ), where T is calculated from the coarse-grained fields. The correlation is positive over most of the domain, with some spatial variation in its strength. A similar result with only slightly weaker values is obtained if the rank correlation coefficient 8 is used instead, indicating the result is not strongly sensitive to the presence of extreme values (of S ζ or of the parametrized forcing). Importantly, regions of weaker correlation do not include the jet region that is of prime interest.
The correlations shown in Fig. 6 a indicate that a portion of the variance of S ζ is captured by S W − W S , which suggests that a parametrized Reynolds stress tensor with this dependence on the flow gradients -i.e. T ( κ) , as defined by Eq. (23) -could form the deterministic part of a stochastic parametrization. It is also possible that the correlation is degraded by the two orders of spatial differentiation required to calculate vorticity tendency from a stress tensor (see Eq. (21) ). An implementation of the parametrization in this model, in which velocity is prognostic, would use the 8 The Spearman rank correlation is the linear correlation of the ranks of the data in the two time series -i.e. for each time series, the data is sorted to determine the rank (1,2,... , n ) of each of the n data points. The purpose is to reduce the influence of outlier data points on the result ( Wilks, 2006 ). momentum forcing calculated from the parametrized stress tensor, which for T ( κ) is given by Eq. (24) . We have also repeated the correlation calculation of Fig. 6 a using data from the same high-resolution idealized QG model run as used in Porta Mana and Zanna (2014) (recall from Section 3.1 that our primitive equations model setup imitates the QG model setup used in that study). The resulting correlation map is similar, with slightly higher correlation values in the jet region (not shown); this may be consistent with the QG model generally exhibiting smoother flow structures than the primitive equations model, which as the vorticity snapshot of Fig. 4 indicates can be prone to noisiness (likely due to internal waves) near the gridscale. Fig. 6 b shows the regression coefficient corresponding to the correlation coefficient shown in Fig. 6 a. Since T (κ ) = 2 κ ( S W − W S ) , the regression coefficient is an estimate of the value of 2 κ.
(Equivalently, Fig. 6 b shows the regression of S ζ against the vorticity forcing that would be obtained by dividing Eq. (25) by 2 κ.) Fig. 6 b shows roughly constant values of 2 κ ≈600 to 800 km 2 over the meandering jet region ( x ≈ 0 -10 0 0 km, y ≈ 20 0 0 -30 0 0 km), and the correlation coefficient ( Fig. 6 a) is also fairly consistent over this region. Some larger regression coefficients occur elsewhere in the domain, many of which appear noisy and in some cases correspond to regions of weaker correlation.
One of the goals of our proposed parametrization is that spatial variations in parametrized eddy forcing should be a function of the resolved flow via only the deformation tensor on which T ( κ) depends, with κ being a fixed constant. That is, when implementing the parametrization, κ should not be prescribed to vary temporally and spatially; rather, variations in parametrized forcing should arise spontaneously from the deformation tensor, S W − W S , which describes how the forcing depends on the resolved flow up to a constant multiplicative factor. The spatial variations of the regression coefficient seen in Fig. 6 b are encouraging in this respect as they are small and the coefficient κ is mostly single signed. Accordingly, we will assume κ to be constant and interpret the spatial variations in Fig. 6 b as representing a combination of sampling variability and physics not captured by T ( κ) .
We can then estimate 2 κ by averaging the regression coefficient across the horizontal domain. Fig. 7 a shows the domainaveraged regression coefficient as a function of the area of the coarse-grained gridbox, which is ( n x ) 2 km 2 where x = 7 . 5 km is the grid scale of the fine-resolution model (e.g. for the n = 4 case shown in Fig. 6 b, the gridbox area is 900 km 2 ). Box-whiskers show the median, upper and lower quartiles, and 5th and 95th percentiles of the spatial distribution of the regression coefficient (i.e. the distribution corresponding to the variability illustrated by Fig. 6 b for the case n = 4 ). It is seen that 2 κ scales linearly with the coarse-grained gridbox area, and this is illustrated more explicitly in Fig. 7 b, which shows the same information as Fig. 7 a but with the regression coefficient for each n divided by ( n x ) 2 km 2 .
This suggests that a scale-aware value of κ can be specified, and that a reasonable estimate is κ ≈ (n x ) 2 / 2 = X 2 / 2 , where X is the horizontal grid scale of a coarse-resolution model. These diagnostic results also indicate κ > 0, in agreement with the analysis of Section 2.3 showing that positive κ is required for the parametrized forcing to behave as an enstrophy sink and to mimic the upgradient momentum flux resulting from a banana-shaped eddy.

Eddy geometry
Here we compare our proposed parametrization to the eddy geometry framework for describing the eddy stress that is explored  (a) For coarse-grainings ranging from n = 2 to 12 times the 7.5 km resolution of the high-resolution model (i.e. from 15 to 90 km), the regression coefficient (which is shown in Fig. 6 b for the case n = 4 ) is averaged over the horizontal domain, excluding the region within 100 km of the boundaries (so as to avoid possible edge artefacts). The domain-averaged regression coefficient is then plotted against the area of the coarse-grained gridbox (black line and dots). The grey box-whiskers indicate the range of regression coefficients found over the horizontal domain. The box shows the lower and upper quartiles and the median (i.e. 25th, 50th and 75th percentiles of the data) and the whiskers show the 5th and 95th percentiles. (b) As in (a), but with the regression coefficients scaled by the area of the coarse-grained gridbox.
in a number of studies (e.g. Hoskins et al., 1983;Marshall et al., 2012;Waterman and Hoskins, 2013;Waterman and Lilly, 2015 ). The mean-flow forcing due to eddies is given by the Reynolds stress tensor where the bar denotes the averaging over eddy phase that is used to define the mean flow. Previous studies ( Hoskins et al., 1983;Waterman and Lilly, 2015 ) have indicated that C can be separated into isotropic and anisotropic parts such that where K = 1 2 ( u u + v v ) is the eddy kinetic energy, M = 1 2 ( u u − v v ) is the excess eddy kinetic energy in the x direction over the y direction, and N = u v is the covariance between fluctuations in the two directions. The second expression rewrites the anisotropic part of the tensor in terms of a magnitude L = M 2 + N 2 and direction θ = 1 2 tan −1 (N/M) . The notation corresponds to Waterman and Lilly (2015) , with L referred as the ellipse anisotropy, and θ the orientation of the principal axes of C , which are the coordinate axes in which the off-diagonal components of C are zero (i.e. a rotated frame in which N = 0 ).
The decomposition in terms of anisotropy and orientation parameters L and θ is valid for any tensor of the above form, which from Section 2.1 can be written generally as a I + bR + + cR × . Therefore we can write a I + bR + + cR × = a I + b 2 + c 2 cos 2 φ sin 2 φ sin 2 φ − cos 2 φ where φ = 1 2 tan −1 (c/b) . For the Reynolds stress tensor we have a = K, b = M, c = N. From Eq. (23) , T (κ ) The similarity between the expressions for C and T ( κ) can be made more explicit by expressing the stretching and shearing de-formations in terms of the total deformation, δ ≡ ˜ D 2 + D 2 , as ˜ D = δ cos 2 γ , D = δ sin 2 γ , where γ = 1 2 tan −1 (D/ ˜ D ) is the angle between the axis of dilation and the x -axis (as defined in Section 2.1 ). From Eq. (23) we can then express T ( κ) as Using the fact that sin (φ + π 2 ) = cos φ and cos (φ + π 2 ) = − sin φ for any angle φ, Eq. (53) would be equivalent 9 to the anisotropic part of Eq. (51) if L = κ| ζ | δ and θ = γ ± π 4 , with the positive sign corresponding to ζ > 0 and negative sign to ζ < 0.
Vorticity evolution due to any stress tensor is governed by its anisotropic part (see Eq. (21) ). Hence for the Reynolds stress tensor it is the polarization of the eddy kinetic energy, as represented by the ellipse anisotropy L , that is important rather than the total eddy kinetic energy, K ( Waterman and Lilly, 2015 ) 10 . For T ( κ) , which has no isotropic part, κ| ζ | δ plays the role of L . A possible physical interpretation is that, given the presence of large-scale flow deformation δ, the "missing" amount of polarized eddy kinetic energy is proportional to | ζ | δ, representing the effects of the elongation of vorticity patches by the anisotropic flow.
Assuming no spatial variation in κ, the gradient in the first term is just ∇( ζ δ) and the whole expression is then proportional to a single constant κ. Analogously to Eq. (55) and the interpretation suggested above, Eq. (56) can be described as the vorticity flux associated with spatial gradients in both (1) the amount of polarized eddy kinetic energy, hypothesized to be proportional to ζ δ, and (2) the polarization direction, related to the orientation of the axis of dilation.
9 It is the anisotropic part of −C to which T ( κ) is equated. This is because the momentum tendency due to C as defined by Eq. (50) , which is a conventional way to define the Reynolds stress tensor, is −∇ · C, but we have defined the sign of our parametrized stress tensor so that its momentum tendency is ∇ · T ( κ) . 10 In the principal axes frame of C , for which the off-diagonal component of C is zero, we have L = √ M 2 = 1 2 ( u u − v v ) . Since L represents the difference between the eddy kinetic energies in the two orthogonal directions of the principal axes frame, we refer to it as the polarization of the eddy kinetic energy.

Nonlinear gradients
Here we describe another way in which C and T ( κ) are related. When coarse graining the flow in the high-resolution model, the eddies are defined as the deviation from the coarse-grained average. We may then approximate the eddying flow using a Taylor series expansion, where u x , etc, are the spatial gradients of the coarse-grained flow. The Reynolds stress tensor may then be written, before averaging, as Taking the mean over a square coarse-graining box of grid cells, x 2 = y 2 = l 2 where l is a length scale characterizing the size of the coarse-graining box, and xy = 0 . Then From Section 2.1 , S 2 − W 2 = 1 4 (δ 2 − σ 2 + ζ 2 ) I + σ S , which yields no contribution to the vorticity tendency from the I -term, and for incompressible flow ( σ = 0 ) no contribution from the S -term. The contribution of Eq. (59) to the vorticity budget should then vary as S W − W S , which is the same deformation tensor on which T ( κ) depends ( Eq. (23) ). Furthermore, Eq. (59) indicates that the value of κ should scale with the area of the coarse-graining box, l 2 , which is consistent with the behaviour described in Section 3.3 ( Fig. 7 ) if l ≈ n x . Hence the tensor S W − W S should approximate the diagnosed behaviour of the Reynolds stress tensor, which can be tested by directly comparing tensor components calculated from the coarsegrained flow. Fig. 8 shows the temporal correlation coefficient between the off-diagonal tensor component of the Reynolds stress tensor (i.e. N = u v in Eq. (51 )) and of the tensor S W − W S (i.e. c = − 1 2 ζ˜ D in Eq. (15) ). For both of the coarse grainings shown ( n = 4 in Fig. 8 a and n = 8 in Fig. 8 b) there is overall a strong correlation between the tensor components, and as expected from Eq. (59) the correlation is negative. A similar strong correlation is obtained for the tensor on-diagonal components (i.e. M = 1 2 u u − 1 2 v v in Eq. (51) correlated against b = 1 2 ζ D in Eq. (15) ; not shown). The overall strength of the correlation degrades with increasing coarsegrain gridbox size, as expected due to the Taylor series becoming a poorer approximation as n increases. Somewhat weaker correlations than those seen in Fig. 8 were obtained in the direct comparison of the vorticity forcing ( Fig. 6 a), which is subject to numerical degradation since two orders of spatial differentiation are required to compute vorticity forcing from a stress tensor (see Eqs. (21) and ( 25 )). Fig. 8 gives further support to the use of T ( κ) as an approximation to the part of the Reynolds stress tensor that affects the vorticity budget (i.e. the second term in Eq. (59) ). Since only a single order of differentiation is required to calculate the forcing in models that prognostically evolve momentum (see Eqs. (17) and ( 24 )), less numerical degradation would be expected in such models. (51) ) and the tensor T = S W − W S (i.e. c = − 1 Since the Taylor series as defined by Eq. (57) approximates the spatial structure of the eddies by definition, a possible interpretation of Eq. (59) is that favourable correlations between the forcing by T ( κ) and the actual eddy forcing are due to T ( κ) being equivalent to the part of the approximated Reynolds stress tensor relevant to the vorticity budget. Any potential parametrization for use in eddy-permitting models will somehow mimic the structure of the eddies, i.e. the spatial structure of the advective term in the fluid equations, since it is the behaviour of the advective term that the parametrization is trying to approximate. It is therefore encouraging that the parametrization is consistent with the approximated spatial structure of the flow as represented by a Taylor series expansion, but it can also be viewed as a filtered version that isolates the vorticity budget contribution of the Reynolds stress tensor. The use of a Taylor series has some precedent in the turbulence literature, where it is referred to as the nonlinear gradient model ( Meneveau and Katz, 20 0 0;Nadiga, 20 08;Nadiga and Bouchet, 2011 ). Its use has been characterized as assuming a "self-similarity" of the eddies, whereby the spatial structure of the large-scale (i.e. coarse-grained) flow can be used to predict the spatial structure of the smaller-scale, unresolved flow. However, Sections 2.3 and 2.4 indicated that T ( κ) has distinct properties, which are not evident from expressing the coarse-grained Reynolds stress tensor as a simple Taylor series. These include energy conservation, enstrophy dissipation, and behaving qualitatively as an anisotropic viscosity that allows both strengthening and weakening of vorticity gradients. These properties suggest that T ( κ) can contribute to the upgradient fluxes required to sharpen and maintain realistic jet structures.

Conclusions
A parametrization of the horizontal Reynolds stresses due to mesoscale eddies for use in eddy-permitting ocean models has been proposed. It is applicable to eddy-permitting models because it depends on the flow deformation and therefore requires that some eddy activity be resolved by the model. Numerical simulations with an idealized eddy-resolving primitive equations model indicate that the parametrization can partially mimic the coarsegrained eddy vorticity forcing. Since the parametrization is ex-pressed in terms of a stress tensor, it is equally applicable to models that prognostically evolve momentum or vorticity (see Eqs. (24) or ( 25 ), respectively, for the explicit expressions giving the parametrized tendencies).
The parametrization has been developed within a twodimensional framework, by considering momentum fluxes and neglecting vertical motion and buoyancy effects that are important in the real ocean. Since the gridscale of eddy-permitting ocean models is near the deformation radius, the buoyancy term should not be large and a barotropic approach is reasonable; nevertheless, future work could explore possible 3D effects. Diagnostic analysis using more realistic ocean states may also be useful.
Key properties of the parametrization include small-scale enstrophy dissipation, conservation of energy (hence no exacerbation of spurious small-scale energy dissipation), and the possibility of upgradient momentum fluxes. Moreover, our model diagnostics indicate that the constant coefficient κ can be specified as having no spatial or temporal variation itself, and scales with the model resolution. To verify whether these benefits can be realized in practice, without other issues arising such as undesirable physical properties or numerical instability, requires implementation of the parametrization, which is being addressed in a separate study (Bachman et al., in preparation). Implementation is not expected to incur significant additional computational overhead in a model since no higher spatial derivatives than those used in conventional eddy viscosity are required.