## Abstract

Given its importance in parametrizing eddies, we consider the orientation of eddy flux of potential vorticity (PV) in geostrophic turbulence. We take two different points of view, a classical ensemble- or time-average point of view and a second scale decomposition point of view. A net alignment of the eddy flux of PV with the appropriate mean gradient or the large-scale gradient of PV is required. However, we find this alignment to be very weak. A key finding of our study is that in the scale decomposition approach, there is a strong correlation between the eddy flux and a nonlinear combination of resolved gradients. This strong correlation is absent in the classical decomposition. This finding points to a new model to parametrize the effects of eddies in global ocean circulation.

## 1. Introduction

Ocean circulation is characterized by interactions over a vast range of spatial and temporal scales. Consequently, direct numerical simulation (DNS) of ocean circulation on climate time scales is unlikely in the foreseeable future. The problem of modelling ocean circulation is then one of how best to abstract important physics represented in the full governing equations at a lower cost. Note that while the Navier–Stokes equations are the governing equations for ocean circulation, the previous statement holds for further approximations of the system such as the primitive equations, the quasi-geostrophic (QG) equations and others.

One choice is to average the system over time or ensembles. Such an averaging is called Reynolds averaging (RA) in classical turbulence (Pope 2000). A model based on RA aims to solve for the mean aspects of the flow. Thus, while a model based on RA is capable of capturing many important aspects of a turbulent flow, it is unable to accurately predict spatio-temporal characteristics of the flow. One way to improve on this is to adopt unsteady RA. In this case, averages are considered over time intervals that are large with respect to the time scale of turbulence, but small compared with the variability time scales of interest. Clearly, unsteady RA-based models can be successful for flows with distinct separation of turbulent and variability time scales.

The implicit assumption of statistical stationarity in RA limits the applicability of that approach to modelling ocean circulation since ocean circulation is dominated by extremely rich variability on a wide range of time scales: the lifetime of mesoscale eddies are of the order of a few months, massive basin-wide changes like those due to El Niño occur on the annual to interannual time scale, the paths of important currents like the Kuroshio current change on the decadal time scale, abrupt changes in the thermohaline circulation (THC) with drastic implications for climate, as evidenced in palaeoclimate records, take place on the decadal to centennial time scales, whereas thermodynamic equilibration and gradual changes of the THC occur on the millennial time scale. The unsteady RA approach would, therefore, seem more applicable. Nevertheless, its applicability is limited to oceanic regimes where there is a separation between turbulence time scales and the variability time scales of interest. For example, the unsteady RA approach would be appropriate to parametrize the effects of baroclinic instability and mesoscale eddies, as through Gent–McWilliams parametrization, in simulations that do not resolve the (first) Rossby radius of deformation in coarse-resolution studies of THC variability on the centennial and millennial time scales. However, the unsteady RA approach would be inappropriate to parametrize the effects of subgrid scales in the mesoscale eddy-permitting and eddy-resolving simulations that are now becoming more commonplace.

State-of-the-art ocean general circulation models (OGCMs) are RA based. (For brevity, we will use RA to refer to unsteady RA as well, when there is no ambiguity.) In fact, turbulent fluxes are typically modelled using turbulent diffusivity/viscosity hypotheses. This includes the Gent–McWilliams parametrization. The parametrization of adiabatic flattening of isopycnals by baroclinic instability, as described by Gent & McWilliams (1990), is achieved effectively by turbulent diffusion of layer thickness (buoyancy), after noting that the bolus (advective) transport of scalars results from isopycnal thickness averaging of velocity. Turbulent diffusion of potential vorticity (PV) in the interior accompanied by buoyancy diffusion in surface and bottom layers has also been proposed (Treguier *et al*. 1997) as an alternative to thickness diffusion. As discussed above, this is justifiable where there is a reasonable separation of scales as in centennial time-scale climate change studies. However, with increasing computational resources, the same turbulent viscosity/diffusivity-based RA approaches are being used in mesoscale eddy-permitting and eddy-resolving simulations on the shorter annual to decadal time scales; settings where a separation of time scales is more difficult to justify.

In classical turbulence, the computational cost of DNS of large Reynolds number turbulent flows increases as the cube of the Reynolds number (Pope 2000) and is therefore prohibitive. Further, RA-based approaches have failed to provide detailed local spatio-temporal flow characteristics in a predictive physics-based fashion. On the other hand, it is almost always the case that in fully resolved computations, a disproportionately high fraction of the computational effort is expended on the smaller scales whereas energy is predominantly contained in the larger scales (Pope 2000; Geurts 2004). This has led to the technique of large eddy simulations (LESs) in which the large-scale unsteady motions that are driven by the specifics of the flow geometry and forcing and that are not universal are computed explicitly and the smaller subgrid motions (that are presumably more universal) are modelled (Pope 2000; Geurts 2004). Historically, however, LES had its origin in the modelling of geophysical flows—Smagorinsky model (Pope 2000).

LES is a turbulence modelling approach that is intermediate between DNS- and RA-based approaches. In a typical ocean simulation, the (time varying) eddy flux is calculated from the time-varying mean flow gradient by the parametrization derived from RA. As mentioned earlier, for coarse-resolution models, this might not be a problem because time variation of the resolved scales is (almost) absent. However, for the mesoscale eddy-permitting and eddy-resolving simulations that are becoming increasingly commonplace, such RA-based closures are more problematic since it involves modelling the effect of (fluctuations of) the full range of spatial scales on the appropriate time-mean circulation that is being resolved. On the other hand, in the LES approach, it is only the effects of the unresolved scales on the resolved scales, which have to be modelled, and hence, for these latter cases, LES is conceptually more appropriate as a modelling framework than RA. In this context, we note that it is unfortunate, but not uncommon, that in ocean modelling literature, some authors refer to *Reynolds* stress in the RA approach as *subgrid* stress.

Given the different functions of the model in the two approaches, the structure of these models of turbulence can be expected to be significantly different. To date, intuition about the nature of the turbulent flux of PV in geostrophic turbulence has been built exclusively using RA concepts. In this article, we take the first step towards investigating the nature of the equivalent object that needs to be modelled in an LES setting, the turbulent subgrid flux of PV. In so doing, we find that there are distinct advantages to the LES approach as compared with the traditional RA approach. However, we will see in §2 that it is straightforward to implement models based on scale decomposition in ocean models since the form of the governing equations in the two approaches is similar.

The simplest subgrid closure in LES is that of an eddy or subgrid viscosity. This viscosity may be based on mixing length ideas or the local magnitude of strain as in the classical Smagorinsky model (Pope 2000). Such eddy or subgrid-scale (SGS) viscosity is designed to represent the mean damping of the resolved scales by the unresolved scales. However, the two-way nature of interaction across the cut-off scale leads to a forcing of the large scales due to their nonlinear interactions with sub-filter-scale eddies (Leith 1990). In the framework of LES, the effects of such backscatter can be accounted for by a stochastic subgrid model in conjunction with eddy viscosity (e.g. Leith 1990). As we will see later, in the forward cascade of enstrophy in geostrophic turbulence, backscatter can be large, suggesting the importance of a stochastic component to a putative scalar eddy viscosity-based subgrid closure. On the other hand, we will also see that the eddy flux of PV aligns well with a nonlinear combination of resolved gradients, suggesting that the significant backscatter can be represented rather naturally in such a nonlinear closure.

In our discussions here, it should be clear that the underlying paradigm is one in which the effects of unresolved dynamical scales on resolved scales in large-scale simulations of geophysical turbulent flows are represented by suitable models. However, Majda *et al*. (2005) approach the problem of simulating the atmosphere–ocean system differently and have successfully developed stochastic mode reduction strategies, wherein they are able to describe essential large-scale dynamics using a very small number of large-scale modes.

In this article, we consider eddy-resolving simulations of the classical wind-driven ocean circulation and analyse them from both RA and LES perspectives to elucidate the differences in eddy–mean flow relations that result from different definitions of eddy and mean in these two frameworks. In the rest of the article, we first describe the governing equations that we work with. Next, we apply RA and scale decomposition to the governing equations and motivate the nonlinear gradient model. Finally, we present computational results and end with a discussion.

## 2. Governing equations and decompositions

Consider the equation for the evolution of QG PV in the layered form. Using standard notation,(2.1)where *q* is the PV; * u* is the velocity;

*F*is the forcing;

*D*is the dissipation; and subscript

*i*refers to the layer number. The non-divergent nature of the two-dimensional advecting geostrophic velocity allows for the introduction of a stream function

*ψ*such that

**=**

*u**k*×∇

*ψ*. In (2.1), PV, for example, for a two-layer system, is given aswith index

*i*=1 corresponding to the top layer and 2 to the bottom layer. For a multilayer system with more than two layers, the stretching term in the definition of PV above, has contributions from both the top and bottom interfaces for layers other than the top and bottom layers. The effect of (shallow) topography is felt only by the bottom (

*N*th) layer and enters the definition of PV of the bottom layer as another term

*f*

_{0}/

*H*

_{N}

*Z*

_{b}, where

*Z*

_{b}is the bottom elevation.

For convenience, 1 may be written as(2.2)

### (a) Reynolds averaging

RA proceeds by assuming a decomposition of the formwhere is the ensemble or time mean. With such a decomposition, and . The evolution of mean PV is then given by(2.3)where(2.4)Evolution of mean potential enstrophy is obtained by multiplying the above equation with as(2.5)where and(2.6)Further, since PV is an advected scalar in an inviscid and unforced setting, the evolution equation of its variance (eddy potential enstrophy) is revealing and is given by(2.7)where and is eddy potential enstrophy.

From (2.5) and (2.7), it is clear that is a transfer of potential enstrophy from the mean to the eddies. With steady wind stress, variations of (upper) layer thickness do not change PV forcing in quasi-geostrophy (, where is the undisturbed depth and so *F*′=0). Thus, in a statistically stationary turbulent flow driven by steady forcing, since dissipation is negative definite, and the second term on the left only serves to redistribute, it is clear that has to have a net (i.e. in the domain-integrated sense) positive value.

The requirement of the net positive value for has motivated the commonly used local downgradient closure,(2.8)Clearly, for this approximation to be valid locally, inspection of (2.7) (assuming statistical stationarity and steady forcing again) advection of perturbation potential enstrophy would have to be negligible. As Rhines & Holland (1979) point out, this is the case when the lateral scale of the mean PV field far exceeds the displacement of fluid particles over a few eddy periods. In effect, this is a restatement, from a Lagrangian point of view, of the requirement of a scale separation between the turbulence that is being parametrized and the scales of the flow that are being studied. Unfortunately, this is not the case in eddy-permitting simulations of ocean circulation with their intense western boundary currents and strongly curved flows. Thus, we will see later that while has a net positive value, the above local turbulent viscosity hypothesis is generally not valid.

### (b) Scale decomposition

In LES, the resolution of energy-containing eddies that dominate flow physics is made computationally feasible by introducing a formal scale separation (Pope 2000). The scale separation is achieved by applying a low-pass filter to the original equations. To this end, consider filtering at scale *l* so that fields ** u**,

*q*, etc. can be split into large- and small-scale components aswherethe filter function

*G*is normalized so thatand where the integrations are over the full domain

*D*. In contrast to Reynolds decomposition, however, generally,

*u*_{ll}≠

*u*_{l}and

*u*_{sl}≠0.

Applying such a filter to (2.2) leads to an equation for the evolution of the large-scale component of PV that is the primary object of interest in LES:(2.9)where(2.10)is the turbulent sub-filter PV flux that may in turn be written in terms of the Leonard stress, cross stress and Reynolds stress (Pope 2000) as(2.11)However, while *σ* itself is Galilean invariant, the above Leonard and cross stresses are not Galilean invariant. Thus, when these component stresses are considered individually, the following decomposition, originally due to Germano (1986), is preferable:(2.12)

The filtered equations, which are the object of simulation on a grid with a resolution commensurate with the filter in LES, are then closed by modelling SGS stresses to account for the effect of the unresolved small-scale eddies. In this case, (2.9) will be closed on modelling the turbulent subgrid PV flux *σ*.

We next consider the transfer of enstrophy from large to small scales at every location in physical space. To this end, multiplying (2.9) by *q*_{l} leads to an equation for the evolution of large-scale potential enstrophy (2.13)and where(2.14)The flux on the l.h.s. can only spatially redistribute large-scale potential enstrophy, where as the term represents the flux of potential enstrophy out of the large scales and into the small scales. To wit, the term appears with the opposite sign in the equation governing the evolution of the potential enstrophy involving small scales,(2.15)and where

The irreversible forward cascade of potential enstrophy requires that domain integral of be positive. On integrating (2.13) over the domain, since the dissipation of enstrophy at large scales is small, in a statistically stationary state, enstrophy input by forcing is balanced by a flux of enstrophy out of the large scales and into the small scales. For this to happen, the sub-filter PV flux would have to have a net alignment down the gradient of the large-scale PV in the domain-integrated sense.

In the following, we show that this required alignment between the turbulent sub-filter flux of PV and the gradient of large-scale PV is weak. On the other hand, we find that the turbulent sub-filter flux of PV is remarkably well aligned with a nonlinear combination of large-scale gradients.

In either of these approaches, it is clear that it is only the divergent component of the eddy flux that directly drives mean circulation. However, given the non-unique (arbitrary) nature of a decomposition of the eddy flux into rotational and divergent components, we prefer to analyse the full flux. The exception is the decomposition of Marshall & Shutts (1981) but which we find does not help in general.

## 3. Orientation of the eddy flux of PV

We consider the classic configuration used in Holland & Rhines (1980) to study eddy-induced ocean circulation. It consists of a two-layer ocean basin on a midlatitude beta plane. A characteristic non-uniform ocean stratification is considered in which the undisturbed top layer depth *H*_{1} is 1 km and the bottom layer depth *H*_{2} is 4 km. The lateral geometry consists of a square domain with a latitudinal and a longitudinal extent of 2560 km. The midlatitude beta plane is such that the Rossby deformation radius is approximately 40 km and the grid spacing is 10 km. The forcing consists of a steady double gyre wind stress with a peak amplitude of 1 dyne cm^{−2}. The dissipation consists of a combination of bottom drag and small-scale mixing of relative vorticity that is either Laplacian or biharmonic but with a (Munk) scale defined as (*ν*_{2}/*β*)^{1/3} or (*ν*_{4}/*β*)^{1/5} of approximately 10 km.

An enstrophy and energy-conserving form of spatial discretization is used for the Jacobian and time stepping is carried out using an adaptive fifth-order embedded Runge–Kutta Cash–Karp scheme (Press *et al*. 1992) to minimize time truncation errors. An advantage of this choice of numerical discretization is that the resultant code is nonlinearly stable and allows for simulations at any level of dissipation including no dissipation at all.

The long time-mean circulation and PV in the two layers are shown in figure 1. The phenomenology of this flow is classic and often referred to in connection with the PV homogenization theory. The reader is referred to Holland & Rhines (1980) for details. While the mean circulation in the lower layer is entirely eddy driven, the mean circulation in the upper layer is modified at *O*(1) by the eddies. The ratio of the eddy to mean kinetic energy is approximately 3 in the upper layer and approximately 7 in the lower layer. We have carried out similar computations in three-layer configurations in order to have a layer that is shielded from both direct wind stress curl forcing and bottom friction, and the qualitative nature of the flows remains unchanged.

That the flow is resolved should be clear from figure 2 where the horizontal spectral distributions of the total energy (solid line) and the barotropic (dashed line) and baroclinic (dot-dashed line) kinetic energies are shown. The bulk of the total energy resides in the large-scale sloping of the isopycnals (layer interface) and constitutes the available potential energy of the system, as seen by the difference between the total energy (solid line) and the sum of the baroclinic kinetic energy (dot-dashed line) and the barotropic kinetic energy (dashed line). And the energy flow in the system from available potential energy to baroclinic kinetic energy to barotropic energy follows the classic phenomenological picture of geostrophic turbulence (e.g. Salmon 1998). Considering the non-uniform stratification, the distribution of kinetic energy in the top layer looks more similar to the distribution of kinetic energy in the baroclinic mode and likewise with the lower layer barotropic pair. This is consistent with present interpretation of altimetric data on large scales (Smith & Vallis 2001).

### (a) Inclination to mean or large-scale PV gradient

#### (i) Reynolds decomposition

Figure 3 shows the distribution of the angle between the eddy flux of PV and the time-mean PV gradient using Reynolds decomposition in the two layers. In this case and all other cases to follow, we obtain the angle between two vectors in the usual manner of dividing their dot product by their magnitudes. If either of the magnitudes is of the order of double-precision round-off, the angle is undefined and not included in the analysis. Distribution of an angle is obtained by considering the angle at all interior grid points (unless stated otherwise) at a few different times (typically approx. 10) separated by the eddy-averaging time (e.g. at various *T*, where *τ* is the averaging time).

The required alignment of the eddy flux down the gradient of mean PV is verified in the mean angle in the two plots being greater than *π*/2 (approx. 1.6 and 1.7 in the two layers, respectively). What is surprising, however, is how weak this alignment is, particularly in the top layer. It is clear from figure 3 that the local turbulent viscosity hypothesis (2.8) is generally not valid.

Time averaging of the eddy flux is performed over a period of approximately 3.4 gyre turnaround times where the lateral scale is the domain size and velocity scale is the Sverdrup velocity. Varying the averaging time from 0 (instantaneous) to 34 gyre turnaround times improves the mean inclination in the lower layer from close to *π*/2 to approximately 2.2, with almost no change to the mean inclination in the upper layer (1.6). The averaging over long times such as 34 eddy turnover times would be appropriate for a RA-based steady simulation; for that a local downgradient eddy flux approximation is reasonable for the lower layer only. Averaging over approximately 3.4 eddy turnover times seems to be appropriate for an unsteady RA-based simulation, and, for this reason, we fix averaging time at approximately 3.4 eddy turnover times for all RA-based analysis. For this case, a local downgradient eddy flux approximation is becoming inappropriate for either of the layers.

It is often assumed that wind forcing in the upper layer complicates an examination of the downgradient nature of eddy flux of PV and that only subsurface layers should be investigated in this respect. Equation (2.7) shows that this is not the case, as pointed out by Drijfhout & Hazeleger (2001). Only the variable part of the forcing, and then only that part that is correlated with the eddy variability itself, affects the downgradient eddy flux of PV. Furthermore, forcing in our case is steady. Also, note that the variations of upper layer thickness do not change PV forcing in quasi-geostrophy. Thus, while upper layer wind forcing allows flow there to cross PV contours, the forcing does not directly affect the orientation of the eddy flux of PV with respect to the mean gradient.

It is only the divergent component of the eddy flux of PV that directly drives mean circulation, and so, it would be just as well if only the divergent component of the eddy flux is directed downgradient. The decomposition of the eddy flux into rotational and divergent components has proceeded along two distinct approaches (i) using a Helmholtz decomposition and (ii) considering components along and across the mean gradient in conjunction with a rotational gauge. The first of these approaches leads to a non-unique decomposition owing to the arbitrariness of the boundary conditions that have been specified independently on the rotational and divergent components (e.g. Fox-Kemper *et al*. 2003). This renders the downgradient nature of any such redefined eddy flux questionable and we do not attempt such a decomposition. The non-uniqueness in the second approach is related to the gauge invariance, but has the advantage of allowing for a purely local decomposition. While the reader is referred to Eden *et al*. (2007) for a comprehensive presentation of the various forms this approach has taken, and for additional physical requirements that are used to constrain the indeterminacy, we consider here only the suggestion of Marshall & Shutts (1981).

Marshall & Shutts (1981) suggested that if the mean circulation contours do not deviate much from the mean PV contours, then a two-way balance is possible in (2.7). That the mean advection of perturbation enstrophy could be balanced by a rotational PV flux aligned along contours of perturbation potential enstrophy, leaving the rest of eddy–PV flux to be downgradient after neglecting triple correlations. Thus, from the eddy enstrophy equation with RA (2.7), this amounts to a two-way balance(3.1)(3.2)again after neglecting triple correlations.

Figure 1 displays a significant degree of co-parallelism between contours of time-mean stream function and PV. We, therefore, check the possibility of a better alignment down the gradient of a component of the eddy flux as suggested by Marshall & Shutts (1981). However, since a scatter plot of mean stream function and mean PV displays large scatter, we use (3.1) as a definition, and look at the distribution of(3.3)The two-signed nature of this quantity in figure 4 shows that the alignment to the mean gradient using such a decomposition is not much better either. Note that since we do not have the actual divergent component of the eddy–PV flux (in this decomposition), we cannot determine the inclination of the divergent component to the mean gradient. So, we rely on the fact that negative/positive values of the above quantity correspond to down/upgradient fluxes. In figure 4, only the western central quarter of the domain is considered since that is where the eddy activity primarily is and the two-way balance is expected. The same figure using the full domain does not show any improvement in the downgradient alignment.

#### (ii) Scale decomposition

For the scale decomposition analysis, we choose a filter that corresponds to an inversion of the Helmholtz operator(3.4)where *L*_{f} is the filter width and *q*_{l}=*q* on the boundaries and so also for any of the other variables. We choose a filter width of 63 km, approximately 1.5 times the Rossby deformation radius of 40 km. There are only minor differences on varying the filter width from 40 to 80 km and these differences are not considered further. The nature of the spectrum in figure 2 suggests that our results are not likely to be sensitively dependent on the width of the filter over a wider range of filter widths.

Similarly, figure 5 shows the distribution of the angle between the (sub-filter) eddy flux of PV and the large-scale gradient of PV using the LES formalism in the two layers. The required alignment of the eddy flux down the gradient of mean PV is verified in the mean angle in the two plots being greater than *π*/2 (approx. 1.6 and 1.7 in the two layers, respectively). Again the alignment is weak, and the distribution of angles very similar to the distribution with Reynolds decomposition. With the eddy flux, most often pointed perpendicular to the large-scale gradient of PV, a local downgradient closure is again not justifiable.

### (b) Inclination to a nonlinear combination of gradients

In eddy-permitting simulations, some of the range of scales of turbulence is explicitly resolved. Therefore, information about the structure of turbulence at these scales is readily available. In LES formalism, there is a class of models that attempts to model the smaller unresolved scales of turbulence based on the assumption that the structure of the turbulent velocity field at scales below the filter scale is the same as the structure of the turbulent velocity field at scales just above the filter scale (Meneveau & Katz 2000).

Further expansion of the velocity field in a Taylor series and performing filtering analytically results in(3.5)a quadratic nonlinear combination of resolved gradients for the subgrid model. The interested reader is referred to Meneveau & Katz (2000) for a comprehensive review of the nonlinear gradient model.

Equivalently, the expansion of *u*_{l} and *q*_{l} in the Galilean invariant form of the Leonard stress component of the sub-filter eddy flux of PV (2.12) in a Taylor seriesproduces at the first order(3.6)again a quadratic nonlinear combination of resolved gradients and where *C*_{2l} is the second moment of the filter used. In the two-dimensional context, this model has been derived by Eyink (2001) without the self-similarity assumption, but rather by assuming scale locality of contributions to *σ* at scales smaller than the filter scale, and its use has been investigated by Bouchet (2003) and Chen *et al*. (2003).

The distribution of the angle between sub-filter PV flux vector and the nonlinear vector combination ∇*u*_{l}.∇*q*_{l} is shown for the two layers in figure 6. The strong alignment of the eddy flux of PV in the direction of the nonlinear combination is evident. We also note a small random component not aligned with the nonlinear combination. We anticipate synthesizing an analysis of this stochastic component together with the deterministic component based on the nonlinear combination of gradients to develop a parametrization in the future.

In analogy with the scale decomposition approach, we check for the alignment of Reynolds eddy flux of PV with in figure 7. We note that Olbers *et al*. (2000), in the context of a two-layer wind-forced zonally homogeneous channel and using Reynolds decomposition, consider a form that is somewhat similar (their form would involve a cubic nonlinearity as opposed to a quadratic nonlinearity here). In comparison with the alignment with the gradient in figure 3, this alignment is qualitatively different with a slight preference for the up- and downgradient directions, but the overall alignment is no better and thus does not prove to be an obvious candidate to consider for parametrizations. In this RA context, it would also be interesting to check if with the TRM-G suggestion of Eden *et al*. (2007) for a rotational component of the eddy flux, the alignment between and will improve.

Figure 8 shows instantaneous spatial structure of the alignments (in radians) with ∇*q*_{l} and ∇*u*_{l}*.*∇*q*_{l} in the scale decomposition formalism. In figure 8*b*, poor alignment with ∇*q*_{l} is indicated by an almost equal distribution of red–yellows and green–blues. On the other hand, in figure 8*c*, strong parallel alignment with ∇*u*_{l}*.*∇*q*_{l} is indicated by a predominance of blue–blacks over red–yellows. A similar plot of the time-averaged PV, and alignments of Reynolds eddy flux of PV with and , shows poor alignments in both cases in figure 9.

## 4. Discussion

Given its multiscale nature, DNS of ocean circulation on time scales of interest are unlikely. Modelling the effects of a substantial range of unresolved scale on larger scale circulation is therefore necessary. While RA-based modelling is appropriate for coarse-resolution simulations with resolutions of the order of a 100 km, an LES approach seems to be more appropriate for the eddy-permitting and eddy-resolving simulations that are presently possible for the interannual to decadal time scale.

Using eddy-resolving simulations, we examined the nature of eddy flux of PV in the RA and scale decomposition approaches in an inhomogeneous (basin) setting. In either of these approaches, a net alignment of the eddy flux with the large-scale gradient of PV is required and verified. However, this alignment is found to be weak in both approaches.

To our knowledge, previous analyses of the eddy flux of PV have only employed the RA approach. We do not attempt a comprehensive discussion of that literature here since it is far too extensive. However, the required net alignment in the RA approach and extensive research by the atmospheric community, using two-layer QG dynamics to study homogeneous *β*-plane turbulence in the presence of prescribed vertical shear, has been used as a justification for modelling the eddy flux as *locally* pointing down the mean gradient. Our finding of weak downgradient alignment is in-line with previous such findings (e.g. Holland & Rhines 1980; Olbers *et al*. 2000; Drijfhout & Hazeleger 2001) that when the turbulence is not homogeneous the local downgradient closure is not verified. This is because advection of eddy enstrophy plays a very significant role, particularly so in a basin configuration. Thus, a local downgradient closure may be seen as increasingly inappropriate in going from fully homogeneous settings to zonally homogeneous settings to basin configurations.

On the other hand, the eddy flux of PV in the scale decomposition approach aligns well with a nonlinear combination of large-scale gradients motivated by the notion of similarity of the structure of turbulence below and above the scale of interest. The eddy flux of PV in the Reynolds decomposition approach does not align well with a similar nonlinear combination of mean gradients. This reiterates the appropriateness of the LES approach to mesoscale eddy-permitting and eddy-resolving regimes.

A transformed Eulerian mean (TEM) or temporal residual mean (TRM; e.g. Eden *et al*. 2007 for a discussion) formulation of the mean PV budget relates the eddy PV flux directed perpendicular to the mean PV gradient to an eddy-induced advection velocity. Combining scalar diffusivity with such an eddy-induced advection also leads to a tensorial structure for the diffusivity. The poor alignment of with that we find in our preliminary tests suggests that the tensor diffusivity is not strongly related to the mean velocity deformation tensor. This finding, however, needs to be examined further.

The complications posed by rotational eddy fluxes have long been recognized and various attempts have been made to define a dynamically relevant rotational eddy flux. A comprehensive discussion of a local decomposition can be found in Eden *et al*. (2007). We considered one simple instance of such a decomposition, as suggested by Marshall & Shutts (1981) and found that a balance between the advection of eddy enstrophy and a rotational component of the eddy flux does not improve the alignment between the remaining component and mean PV gradient in the downgradient sense. In a more generalized sense, we could have checked the alignment between this remaining component and , or for that matter also use the TRM-G prescription for the rotational component and check for alignments of the remaining component. These aspects of analysing are, however, of tangential interest to this paper and we did not pursue them.

The good correlation that we obtain between the sub-filter eddy flux of PV and the nonlinear model in the scale decomposition approach is not entirely surprising. Similar good correlations using the nonlinear model have been seen both in three- and two-dimensional turbulent flows (e.g. Meneveau & Katz 2000; Bouchet 2003; Chen *et al*. 2003). In some of these other contexts, the nonlinear model has been recognized to be insufficiently dissipative. This was actually noted in the *a posteriori* sense in the barotropic double gyre ocean circulation context in Holm & Nadiga (2003) as well.

How well a model performs is best assessed by comparing results from simulations that use and do not use the model against available data, some of which may come from appropriately resolved simulations (*a posteriori* testing of LES). However, because results of simulations contain integrated effects of numerical discretization and other artefacts, besides the effects of the model, *a posteriori* testing is not very insightful (e.g. Ghosal 1996; Meneveau & Katz 2000; Holm & Nadiga 2003; Nadiga & Livescu 2007). We have instead, conducted preliminary *a priori* testing, wherein sub-filter turbulent fluxes are checked for correlations against larger scale quantities in a resolved flow. While we find good correlations to a nonlinear combination of gradients, this correlation is neither necessary nor sufficient for LES using this model to perform well. It, however, provides a good physical basis to design parametrizations.

The findings on the orientation of the eddy flux of PV that we present here are of significance to both deterministic and stochastic representation of subgrid processes in geostrophic turbulence. First, from the deterministic modelling point of view, we expect that parametrizations based on the nonlinear model will be better than those based on downgradient closures.

In the scale decomposition approach, the poor alignment between the eddy flux of PV and the large-scale gradient and the peaking of the angle between them at *π*/2 implies the importance of the two-way communication between the sub-filter scales and the larger scales. There is simultaneous (i) drain of enstrophy from large to small scales and (ii) backscatter from small to large scales. These two processes are of comparable importance and there is only a slight difference between the two processes that result in the net downgradient nature of the alignment. Eddy (subgrid) viscosity approaches, on the other hand, end up representing only the small net downgradient nature of the sub-filter resolved-scale interactions. With the situation being similar in the forward energy cascade regime of three-dimensional turbulence, improved subgrid models have been obtained by representing the two processes distinctly (Leith 1990; Chasnov 1991). This forms the basis of stochastic LES in the context of three-dimensional turbulence.

One may similarly anticipate better subgrid models of geostrophic turbulence when both the drain and backscatter processes are represented rather than just the net effect. In fact, Jung *et al*. (2005) have successfully employed such a strategy in an operational atmospheric circulation model to reduce certain model biases and improve the simulation of certain weather regimes. In more idealized oceanic set-ups, Berloff (2005), Nadiga *et al*. (2005), Duan & Nadiga (2007) and Nadiga & Livescu (2007) have performed statistical analyses of sub-filter terms and reported preliminary results from such stochastic sub-filter closures.

Thus, if the subgrid modelling problem is approached from the point of view of (scalar) eddy viscosity, then, a stochastic representation of backscatter has shown promise. However, the good alignment of eddy flux of PV with ∇*u*_{l}.∇*q*_{l} that we find suggests that a parametrization of the eddy flux based on ∇*u*_{l}.∇*q*_{l} has the potential to model the forward and backscatter rather naturally and deterministically.

From the point of view of an eddy viscosity closure, the good correlation between the eddy flux of PV and the large-scale nonlinear gradients implies a tensorial rather than a scalar form for the eddy viscosity. This is understandably so, since in LES there is not a good separation of scales that is required for a scalar eddy viscosity closure to hold (Kraichnan 1967). In analogy with the TEM and TRM approaches, the tensorial form of eddy viscosity implies an eddy-induced advection in addition to diffusion (and anti-diffusion) of PV in the scale decomposition or LES approach.

That the resolved component of velocity deformation determines the structure of the eddy viscosity tensor, implies owing to the incompressible nature of the velocity field, that certain directions will be subject to negative viscosity—an issue that we will address and report on elsewhere. However, we note that the angle between the sub-filter flux and the nonlinear gradient itself is a random variable with a distribution peaked at 0. This implies that there is a small component of the sub-filter eddy flux that is not well represented by the nonlinear gradient. We are studying the nature of this component and anticipate the model for this component to play a significant role in developing a parametrization based on the nonlinear combination of gradients.

State-of-the-art OGCMs are RA based. This is, however, not a handicap for implementing models based on scale decomposition. This is evident in the form of the governing equations (2.3) and (2.9) in the present context in the two approaches. Thus, given an ocean model based on RA, it is straight forward to implement a model based on scale decomposition ideas. It should be noted, however, that while in the RA approach, *Σ* is a model for the fluctuating component of the full spectrum of spatial scales, in the LES approach, *σ* is a model only for the spatial scales that are not resolved. Further, when such a scale decomposition-based model has been implemented in a conventional ocean model, it is important to interpret the variables in the scale decomposition sense rather than in the RA sense.

Finally, we expect that, although the findings of this article are not based on simulations of extensive ranges of parameters, the qualitative nature of these results will hold-up in flow situations where eddies play an important role in shaping the large-scale and/or mean circulation.

## Acknowledgments

I would like to thank Greg Eyink for extensive discussions and Andy Majda for encouraging me to pursue and write-up some early work that I presented at a workshop in Banff. Constructive criticism by Carsten Eden and an anonymous referee of the original manuscript has led to significant improvement in the presentation. B.T.N. was supported in part by the Climate Change Prediction Program of DOE and the LDRD Program at LANL (20030038DR).

## Footnotes

One contribution of 12 to a Theme Issue ‘Stochastic physics and climate modelling’.

- © 2008 The Royal Society