<figure><img src="image_1.png" /></figure>

HAL
open science

# Diffusivity, Kinetic Energy Dissipation, and Closure Theories for the Poleward Eddy Heat Flux

Guillaume Lapeyre, I. Held

► To cite this version:

Guillaume Lapeyre, I. Held. Diffusivity, Kinetic Energy Dissipation, and Closure Theories for the Poleward Eddy Heat Flux. Journal of the Atmospheric Sciences, 2003, 60 (23), pp.2907-2916. 10.1175/1520-0469(2003)0602.0.CO;2. hal-03358685

HAL Id: hal-03358685

https://hal.science/hal-03358685

Submitted on 9 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diffusivity, Kinetic Energy Dissipation, and Closure Theories for the Poleward Eddy Heat Flux

G. LAPEYRE

Atmospheric and Oceanic Sciences Program, Princeton University, Princeton, New Jersey

I. M. HELD

NOAA/Geophysical Fluid Dynamics Laboratory, Princeton, New Jersey

(Manuscript received 19 August 2002, in final form 16 June 2003)

## ABSTRACT

Diffusive eddy closure theory for estimating the poleward heat flux is reexamined and tested in the context of a two-layer homogeneous model. Consideration of the inverse energy cascade induced by baroclinic turbulence on the $\beta$ plane leads to an expression for diffusivity in terms of the kinetic energy dissipation and the $\beta$ effect. A key step in the closure is the identification of this diffusivity with that for potential vorticity in the lower of the two layers in the model. This assumption is then combined with an exact expression relating the diffusivity to the baroclinic energy generation and the mean vertical shear. The theory is closed by identifying the kinetic energy dissipation entering the inverse cascade argument with the baroclinic energy production. It is found that the first constraint in isolation based on inverse cascade arguments between the diffusivity of lower-layer potential vorticity and the kinetic energy dissipation is robust and accurate, whereas the final theory relating diffusivity to vertical shear remains useful but has somewhat degraded accuracy and is more sensitive to model parameters, such as numerical resolution and small-scale dissipation. In the limit of large supercriticality, this theory reduces to that of Held and Larichev. However, it is much more accurate in reproducing numerical results from a two-layer homogeneous model on a $\beta$ plane for the moderate supercriticalities that are typical of model atmospheres. The problems involved in generalizing this result to models with more layers on the vertical or with a continuous stratification are discussed.

## 1. Introduction

The theory for baroclinic eddy heat fluxes attributed to Held and Larichev (1996, hereafter HL) can be recast in a form in which one simultaneously satisfies two constraints between the eddy diffusivity and the rate of baroclinic energy production. The first of these constraints arises from the assumption that the energy-containing eddies stem from an inverse energy cascade that is halted by the $\beta$ effect, as discussed originally by Rhines (1975). From this assumption, one can relate the eddy length and velocity scales to $\epsilon_d$, the rate per unit of mass at which kinetic energy is flowing through the inverse cascade and being dissipated at large horizontal scales by surface friction. From these length and velocity scales, one can estimate a diffusivity from their product:

$$ \mathcal{D} \propto \epsilon_d^{3/5} \beta^{-4/5}. \quad (1) $$

This is the unique way of constructing a kinematic dif-

fusivity, given only $\epsilon_d$ and $\beta$, and this scaling underlies the theory of HL [see their Eqs. (26) and (32)]. Among others, Maltrud and Vallis (1991), Panetta (1993), Spall (2000), and Danilov and Guranie (2002) evaluate aspects of the predictions of eddy length and velocity scales following the Rhines argument. Smith et al. (2002) discuss scaling arguments of this type for the diffusion of a passive tracer by two-dimensional turbulence.

The second half of the theory of HL consists in relating $\epsilon_d$ with the production of available potential energy. In a quasigeostrophic system in which baroclinic production is the only significant eddy source, eddy energy is created through downgradient eddy heat fluxes and is typically dissipated at large scales through friction (vertical turbulent diffusion in surface boundary layers). In the Boussinesq approximation, for example, the baroclinic production is

$$ \epsilon_p = - \left[ \frac{\partial B}{\partial y} \frac{\partial B}{\partial z} \right]. \quad (2) $$

Here, the climate is assumed to be zonally symmetric, $B$ is the zonal mean buoyancy, $v'b'$ is the eddy horizontal buoyancy flux, and the brackets denote an av-

<footnote>Corresponding author address: G. Lapeyre, LPO, IFREMER, BP 70, Plouzané 29280, France.
E-mail: glapeyre@ifremer.fr</footnote>
erage over latitude and height. If small-scale dissipation is negligible, the production of available potential energy is equal to the energy dissipation at large scales.

To close the theory, one needs to relate the buoyancy flux to the diffusivity in (1). In the two-layer model of HL, one will obtain different answers, in general, if one chooses to diffuse buoyancy, lower-layer potential vorticity, or upper-layer potential vorticity with the value of the diffusivity set by (1). However, in the limit of strong supercriticality, in which the vertical shear is much larger than the minimum shear required for instability, the distinction between diffusing buoyancy and potential vorticity disappears in this two-layer model. We return to this point below as it is a key ingredient in the modification proposed here to the theory of HL. Putting this question aside for the moment, and setting $v'b' = -\mathcal{D}\partial B/\partial y$, we write

$$ \epsilon_p = \left[ \mathcal{D} \frac{(\partial B/\partial y)^2}{\partial B/\partial z} \right] = \left[ \frac{\mathcal{D}}{\mathcal{T}^2} \right], \quad (3) $$

where $\mathcal{T}^{-1}$ is an inverse timescale proportional to the local Eady growth rate $fRi^{-1/2}$. If one can now replace $\mathcal{D}$ and $\mathcal{T}$ with appropriately weighted global means, equate the resulting value of $\mathcal{D}$ with that in (1), and equate $\epsilon_p$ and $\epsilon_d$, one obtains the result in HL:

$$ \mathcal{D} \propto \frac{1}{\beta^2 \mathcal{T}^3}. \quad (4) $$

Recently, Barry et al. (2002) have examined the scaling properties of the eddy heat transport for a set of general circulation experiments by varying external parameters, and their results support the connection between diffusivity and $\epsilon$ in the more realistic context of atmospheric general circulation models. However, they argue that their results are not well fit by (4), although it is implicit that (1) is well satisfied. There are several reasons why the theory of HL might break down when applied to comprehensive GCMs. Some of these can be addressed in the context of dry quasigeostrophic dynamics, and others relate especially to the extra complexities introduced by moisture. This paper is devoted to some of the dry quasigeostrophic issues. In particular we reexamine the dry two-layer, homogeneous model of HL and test the different parts of the theory that leads to (4). We provide a modified heuristic theory that fits the numerical results much more closely. This modified theory is of particular relevance when the supercriticality of the flow is relatively weak. In the process we make connections with some results of Pan van and Held (1996) for horizontally inhomogeneous quasigeostrophic flows. The impact of moist thermodynamics on these issues will be addressed in a further study.

Through the different sections of this paper, we discuss and assess the different steps leading to (4). We use the numerical setting of a two-layer model but believe that similar issues arise in other models. Indeed,

we briefly refer to the three-layer model of Held and O'Brien (1992) to help make this point. The numerical model is presented in the appendix. In section 2, we examine the arguments leading to (1) and test its accuracy. In section 3, we discuss the relation between the energy production and dissipation and also how to relate our present theory to the results of Barry et al. (2002).

## 2. Relation between diffusivity and cascading eddies

### a. Inverse cascade scalings

The relation (1) is tightly linked to the properties of the barotropic inverse cascade of energy on the $\beta$ plane. Given assumptions about this inverse cascade, several relations exist between the characteristic eddy velocity scale $V$ and length scale $k_0^{-1}$. The first of these assumes that there exists an inertial range with the universal $C\epsilon_d^{2/3}k^{-5/3}$ form. This assumption means that the cascade is localized in the spectral sense. From this, just as in the case of the inertial range in three-dimensional turbulence, one can relate the rate at which energy is flowing through the inverse cascade to $V$ and $k_0$:

$$ \epsilon_d \propto V^3 k_0. \quad (5) $$

Here, $\epsilon_d$ is the dissipation of barotropic energy by the large-scale drag. See, for example, Smith and Vallis (2002) for a discussion of this relation in the two-dimensional context.

We can use this relation to define a length scale $k_B^{-1}$ such that

$$ k_B = \epsilon_d/V^3. \quad (6) $$

The second assumption is that the eddies are characterized by the Rhines scale $k_R^{-1}$, that is,

$$ k_R = (\beta/V)^{1/2}, \quad (7) $$

where $V$ is the rms velocity. Rhines (1975; see also Panetta 1993) introduced this scale to predict the stopping scale of $\beta$-plane turbulence. However, as noted originally by Rhines and as discussed by Danilov and Gurarie (2002) and Smith et al. (2002), the $\beta$ effect alone cannot arrest the cascade completely; less efficient energy transfer to zonal flows of larger and larger meridional scale continues until arrested by frictional drag. If we want to identify $V$ and $k_R$ with scales relevant for meridional mixing, we need to assume that it is the point of transition from a more or less isotropic cascade to the anisotropic interaction with jets that is the determining characteristic of the flow, and not the larger jet scales that evolve if the damping is weak enough. See Smith et al. (2002) for support for this picture.

One can introduce another scale involving $\beta$ and $\epsilon_d$:

$$ k_B = \beta^{3/5} \epsilon_d^{-1/5}. \quad (8) $$

If $k_B$ and $k_R$ are proportional to each other, then they
are also proportional to $k_β$. Also, assuming that $ε_d$ and $β$ are the only relevant quantities and by a dimensional analysis, one finds that the three scales are proportional. Vallis and Maltrud (1993) argue that (8) has advantages over (7) in predicting the stopping scale for barotropic turbulence as it is a Galilean invariant. However, the advantage of (7) is that it gives information on the energy content of the flow. This can be seen in the context of 2D turbulence when one varies bottom drag keeping fixed $ε$ (see, e.g., Smith et al. 2002).

Assuming that the eddy scales (6), (7), and (8) are the same, one obtains immediately a characteristic diffusivity $V_{k_0^{-1}}$ defined by (1). The differences between these scales can be examined in Fig. 1. This figure shows results from our 256² and 1024² simulations (see the appendix) and results from the experiments of HL, as a function of supercriticality $ξ$ (shear divided by the critical shear required for inviscid instability). For clarity, the scales have been multiplied by appropriate constants so that they equal the Rhines scale $k_R^{-1}$ at $ξ = 4$. Also, $ε_d$ was evaluated as the dissipation of barotropic energy by Ekman friction. One can see that the different scales match each other for $ξ < 5$. For strongly supercritical flows ($ξ > 5$), the two Rhines scales $k_R$ and $k_β$ remain very similar, but $k_B$ strongly overestimates the other scales in the HL experiments and underestimates them in our 1024² runs. The actual relationship between these scales could involve the properties of the inverse cascade, such as the presence of inertial range and the effect of large-scale dissipation. The equality between $k_R$ and $k_β$ contrasts with the one-layer $β$-plane case studied by Vallis and Maltrud (1993), who found important differences between these scales. In the present model, we can understand the similarity between $k_R$ and $k_β$ as follows. If one assumes

$$k_R/k_B = \xi^x,$$

one finds

$$k_R/k_β = \xi^{-x/5}.$$

Thus, even if one finds a large difference between $k_R$ and $k_B$, the difference between the two Rhines scales $k_R$ and $k_β$ will remain small.

To assess further if a theory for eddy scales can be expressed only in terms of $ε_d$ and $β$, we can examine how the velocity scale $V$ scales as a function of $ε_d$ and $β$. From $k_R = k_β$, one would have

$$V \propto \beta^{-1/5} \varepsilon_d^{2/5}. \qquad (9)$$

As can be seen in Fig. 1d, this works indeed very well for all of the experiments, even when $k_B$ departs from $k_β$. Assuming as before that $k_R/k_B = \xi^x$, we find

$$V = \varepsilon_d^{2/5} \beta^{-1/5} \xi^{2x/5}.$$

This relation shows that the velocity scaling is only weakly sensitive to the departure of $k_B$ from $k_β$. Such a departure is expected to occur if the model does not

produce the Kolmogorov spectral shape, that is, if it does not have a well-defined inverse energy cascade.

The barotropic energy spectra from our two-layer model are presented on Fig. 2. The spectra are steeper than $k^{-5/3}$. This is presumably related to the fact that the barotropic energy production (i.e., the baroclinic to barotropic energy transfer) is not localized in spectral space but is distributed over a wide range of scales, as was shown by HL. It could also result from nonlocality of energy transfers, and/or large-scale anisotropic structures (e.g., zonal jets). This helps explain why $k_B$ behaves differently from the other scales.

These spectra help one to understand some important differences between the different sets of simulations. In our 256² case, the friction is not large enough to prevent the inverse cascade from reaching the largest scale when the supercriticality is large. So the barotropic kinetic energy peaks at the domain scale (see Fig. 2a). Because of this piling up of energy at the domain scale, the eddies are more energetic than predicted: $V/(β^{-1/5}ε_d^{2/5})$ increases as $ξ$ increases (see inset in Fig. 1d). On the other hand, in the 1024² case, bottom friction is large enough that the stopping scale is not set anymore by $β$ alone: as can be seen in Fig. 2b, the eddy length scale is no longer set by the Rhines scale. Both $β$ and drag are important in this case. Presumably, the velocity scale is no longer set by $ε_d$ and $β$ and should also depend on friction. This behavior is consistent with the physical fields as we find that the eddies are almost isotropic in this case (implying that $β$ is not the essential parameter). As a result, the energy-containing eddies have less energy for a given energy generation rate than predicted: $V/(β^{-1/5}ε_d^{2/5})$ decreases as $ξ$ increases (inset in Fig. 1d). These results show that, for a given barotropic energy injection rate, the energy of the eddies may vary substantially, depending on the efficiency of the drag in arresting the cascade. Danilov and Gurarie (2001) show that for barotropic f-plane turbulence, large-scale drag can be responsible for non-universal properties of the inverse cascade. Here, another source of nonuniversality is the broad distribution of the barotropic energy generation. In consequence, the precise equality of all three scales defined above may be difficult to realize.

### b. Diffusivity scaling

Introducing a mixing scale $k_d^{-1}$ such that

$$k_d = V/\mathcal{D}, \qquad (10)$$

and equating it with $k_β$ and $k_R$, we obtain

$$\mathcal{D} = c\beta^{-4/5}\varepsilon_d^{3/5}. \qquad (11)$$

We could also derive this relation from $\mathcal{D} = V_{k_0^{-1}}$, using the length scale $k_0 = k_β$ and the velocity scale defined by (9).

One now needs to equate the diffusivity $\mathcal{D}$ either with the diffusivity of buoyancy $\mathcal{D}_b = -[v'b']/\partial_y B$ or of
<figure><img src="image_3.png" /><figcaption>FIG. 1. The length scales ($k_\beta\lambda)^{-1}$ (circles), ($k_R\lambda)^{-1}$ (plus symbols), ($k_\beta\lambda)^{-1}$ (crosses), as defined in (6)–(8), plotted as a function of supercriticality $\xi$ for the three experiments (a) 256² runs, (b) 256² HL runs, and (c) 1024² runs. The scales have been multiplied by appropriate constants so that they equal the Rhines scale at $\xi = 4$. (d) Barotropic rms velocity as a function of prediction $\varepsilon_d^{2/5}\beta^{-1/5}$. Inset corresponds to the ratio of barotropic velocity and prediction as a function of $\xi$.</figcaption></figure>

upper-/lower-layer potential vorticity (PV; i.e., $\mathcal{D}_1 = -[v'q']_1/\partial_y Q_1$ and $\mathcal{D}_2 = -[v'q']_2/\partial_y Q_2$). In the two-layer model, the potential vorticity flux integrates to zero over the domain:

$$ \sum_k H_k [v'q']_k = 0. $$

Here, $k = 1$ refers to the upper layer, $k = 2$ to the lower layer, $H_k$ is the mean depth of the layer, and $[v'q']_k$ is the poleward eddy PV flux. Brackets now refer to a horizontal average. The uniform PV flux in the upper layer is equal and opposite to the uniform PV flux in the lower layer. The mean PV gradients are $\partial_y Q_1 = \beta(1 + \xi)$ in the upper layer and $\partial_y Q_2 = \beta(1 - \xi)$ in the lower layer. In this model the ratio of the diffusivities in the two layers is constrained by the ratio of the mean PV gradients (Vallis 1988):

$$ \frac{\mathcal{D}_1}{\mathcal{D}_2} = \frac{\xi - 1}{\xi + 1}. $$

Held and Larichev (1996) considered the limit of $\xi \to \infty$, in which case the PV gradients in either layer are equal to the horizontal buoyancy gradient, to within a constant, and we do not need to choose between diffusing PV or buoyancy. The question is how best to proceed for general values of $\xi$. Our proposal is not based on a fully coherent theory, but given the very good fit that we obtain with numerical results, we believe that the assumptions resulting in this theory deserve close scrutiny. Here our choice is to use the lower-layer diffusivity as our best guess; that is,

$$ \mathcal{D} \propto \mathcal{D}_2 = \frac{[v'q']_2}{\beta(\xi - 1)}. $$
<figure><img src="image_4.png" /><figcaption>FIG. 2. Spectra of barotropic kinetic energy for different values of supercriticality (a) for our 256² simulations and (b) for our 1024² simulations. The upper curves correspond to larger supercriticality. The stars correspond to the Rhines scale $k_R = (\beta/2V)^{1/2}$, where V is the barotropic velocity rms. The figure insets are the compensated energy spectra (i.e., multiplied by $k^{-5/3}$). The circles still correspond to $k_R$.</figcaption></figure>

As discussed in numerous studies (see, e.g., Tréguier et al. 1997; Marshall et al. 1999), it seems more appropriate to diffuse a quantity conserved following the geostrophic flow (when adiabatic and inviscid), such as PV rather than buoyancy. The choice of the vertical level to use for diffusing PV is less clear cut. In the two-layer model, in the presence of the $\beta$ effect, typical eddy phase speeds are such that $U-c$ is smaller in the lower than the upper layer, so the lower layer more nearly serves as the steering level for these eddies. As a result, meridional particle displacements are larger for a given streamfunction perturbation in the lower layer, and the eddies can be expected to mix PV more efficiently. The upper layer is characterized by coherent zonal jets in which large values of $U-c$ and PV gradients tend to restrict the large-scale meridional transport. The lower layer flow creates mixing that is more plausibly described by simple diffusion.

Using the lower-layer PV diffusivity for $\mathcal{D}$ in (1) does work quite well for all the simulations, as can be seen on Fig. 3. Actual and predicted diffusivities have been nondimensioned by $U\lambda$, where $\lambda$ is the deformation radius, and we multiplied each set of predictions by a constant chosen to best fit with the numerical results. The value of the constant is provided in the caption. As evident in the figure, there is a clear range of values where this scaling for the lower-layer diffusivity works quite well, and this agreement is improved for the 1024² simulations. The fact that this theory works so well may be partly due to the fact that, if there is some departure from the scaling expected for a pure inertial range, that is, $k_R/k_B = \xi^{\alpha}$, this gives a diffusivity $\mathcal{D} \approx \beta^{-4/5}\epsilon_d^{3/5}\xi^{3\alpha/5}$. Thus, there is a rather weak dependence on discrepancies of this type in the length scales.

We have also observed that the result is not very

sensitive to the exact definition of $\epsilon$. In the above we have used the energy lost to surface friction from the barotropic mode, but an equally good fit is obtained if one uses the total energy dissipation by surface friction (in the baroclinic as well as the barotropic mode). Replacing $\epsilon_d$ with the energy production $\epsilon_p$, which balances energy lost through both large-scale drag and small-scale mixing, does not alter the fit either (not shown). The 3/5 power law in the expression for $\mathcal{D}$ makes the result rather insensitive to these types of differences.

If one instead plots the diffusivity for the upper-layer PV against $\beta^{-4/5}\epsilon_d^{3/5}$, the resulting fit is clearly not as good, as shown in Fig. 3b. It does appear to be appropriate, in this two-layer model, to associate the lower-layer diffusivity of PV, and not the upper-layer diffusivity, with the energy throughflow in this way.

Concerning the lower-layer diffusivity scaling, we observe some departure at both ends of the range of supercriticality. For small $\xi$, the diffusivity increases with decreasing $\xi$ despite the fact that energy decreases (not shown). We relate this behavior to the effect of bottom friction and frictional instability that seem to increase the PV fluxes, as we find that a passive tracer advected by the bottom layer has a smaller diffusivity than does PV. For large $\xi$, the simulations present different behaviors: the 256² simulations show that the diffusivity is larger than its predicted value, whereas the opposite is true for the 1024² case. The causes of this discrepancy relate to Ekman friction and numerical resolution: in the 256² case, we have seen that friction is inefficient in arresting the inverse cascade, so that eddies feel the size of the domain at large $\xi$ and are more energetic. This leads to a diffusivity larger than predicted. On the other hand, in the 1024² case, bottom friction is quite efficient
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Nondimensional diffusivity $D/U\lambda$ as a function of $c\epsilon^{3/5}\beta^{-4/5}/U\lambda$. The crosses are the data from HL varying $\beta$, whereas the circles are our 256² dataset varying $\lambda$. The constant $c$ was adjusted for each case to give the best fit for the scaling (15); that is, $c = 1$ for the crosses and $c = 1.44$ for the plus symbols. The solid line corresponds to the perfect fit. (b) Nondimensional diffusivity as a function of $c\epsilon^{3/5}\beta^{-4/5}/U\lambda$ for the 1024² simulations. The stars correspond to the lower-layer diffusivity and the plus signs to the upper layer; $c$ equals 1.25 in this case.</figcaption></figure>

in arresting the cascade and the eddies are less energetic; the diffusivity is also smaller than predicted by focusing on the effect of $\beta$ alone.

A last question concerning relation (11) is how to choose the appropriate constant $c$. This constant might have some dependence on large-scale friction because the cascade cannot be fully stopped by the $\beta$ effect alone [see discussion in Smith and Vallis (2002); Danilov and Gurarie (2002)]. It may also depend on small-scale dissipation if this dissipation is large enough to affect the energy content of the eddies and the precise shape of the energy spectra. Actually, the spectral shape of the energy could be relevant for the constant $c$ since the eddy velocity $V$ and the cascade rate $\epsilon$ are precisely linked through this shape: the universal form $C\epsilon^{2/3}k^{-5/3}$ would lead to (5). Thus, the departure from this form (due to $\beta$ or the absence of inertial range) could affect the constant $c$ as is observed. However, we have not been able to make a quantitative connection between the two.

In the results described above, the two layers in the model have equal mean thicknesses. We have also examined a two-layer model with layers of unequal thicknesses. The results are equally good: we still observe relation (1) to be valid, using the lower-layer PV diffusivity (not shown).

The question is how to proceed in a multilayer model, or with a continuous stratification, within this horizontally homogeneous framework. Is there a particular layer in which the diffusivity for PV can be predicted with the same scaling that works for the lower layer in the two-layer model? Two hypotheses are suggested: use

the eddy steering level, or use the lowest model level. As discussed above, eddies are known to be efficiently mixed at the steering level where Rossby waves have the same speed as the zonal flow [see Haynes et al. (2001) for an illustration]. In this case, the PV diffusivity at the steering level would intervene in (1). For the continuously stratified case in which the source of the instability is the surface temperature gradient, it may be best to think of temperature at the surface as diffusing according to (1). Indeed surface temperature in continuously stratified baroclinic instability often plays a similar role as does the lower-layer PV in the two-layer model.

To begin testing this scaling in a more general setting, we use the model of Held and O’Brien (1992), which generalizes the present model to the case of three layers of equal depth and is forced by a mean zonal wind shear whose distribution is changed in each simulation. From the data provided in Held and O’Brien (1992), the PV diffusivity can be computed in each layer and compared with prediction (1). The total generation of eddy energy was used instead of the energy dissipation in this computation. Figure 4 shows that (1) fits quite well for the lower-layer diffusivity. This confirms that, in this case as well, there exists a particular layer where we are best able to equate $D$ with the PV diffusivity. Only further study can elucidate if it is an interior layer closest to the typical steering level or the surface layer.

Even if we can predict the diffusivity of a particular layer, one still needs the vertical distribution of the diffusivity for a full theory of the eddy transport, even in horizontally homogeneous models such as these (see
<figure><img src="image_6.png" /><figcaption>FIG. 4. Potential vorticity diffusivities as a function of prediction (11) using the three-layer model of Held and O'Brien (1992). The PV diffusivities were shifted to separate each layer.</figcaption></figure>

Tréguier et al. 1997; Smith and Vallis 2002). We can hope to write

$$ \mathcal{D}(z) = \gamma(z)\mathcal{D}_0, \qquad (12) $$

where $\gamma$ is the distribution of the diffusivity in the vertical, $\gamma(z_0) = 1$ for some $z_0$ at which the scaling described above works best, and $\mathcal{D}_0$ satisfies

$$ \mathcal{D}_0 = c\beta^{-4/5}\varepsilon_d^{3/5}. $$

Because of zero total PV flux, one has the extra constraint

$$ [v'q'] = \mathcal{D}_0 \int \gamma(z) \partial_y Q(z) dz = 0, $$

but it is only in the two-layer model that this constraint, together with the expression above for the diffusivity in a particular layer, fully determines the PV fluxes [without using (12)].

### 3. Energy production

#### a. Relation between diffusivity and supercriticality

To close the theory in the two-layer model, one has to find another relation between $\varepsilon_d$ and $\mathcal{D}$. The generation of eddy kinetic energy plus available potential energy per unit mass can be written as

$$ \varepsilon_p = -\frac{1}{H_1 + H_2} \sum_{k=1}^{2} H_k U_k [v'q']_k. \qquad (13) $$

This expression reduces to the two-layer version of (2) in the horizontally homogeneous case (in which there is no transfer of eddy kinetic to mean kinetic energy). In a statistically steady state, the production of baro-

clinic energy, $\varepsilon_p$ is exactly compensated by the dissipation of baroclinic and barotropic energy by surface friction and small-scale mixing. Given the insensitivity of our tests of the relation (1) to the precise definition of $\varepsilon_d$, we proceed by identifying $\varepsilon_p$ with $\varepsilon_d$ in (1). In this picture, we assume that all the energy produced baroclinically cascades barotropically to large scales. Then we have

$$ \varepsilon_d \approx \varepsilon_p = \mathcal{U}[v'q']_2 = \mathcal{U}\beta\mathcal{D}_2(\xi - 1), \qquad (14) $$

where $\mathcal{U} \equiv (U_1 - U_2)/2$ so that $\xi = \mathcal{U}/(\beta\lambda^2)$. This results in the following expression for the diffusivity:

$$ \mathcal{D}_2 \propto \mathcal{U}\lambda\xi^2(1 - \xi^{-1})^{3/2} = \beta^{-2}\mathcal{U}^3\lambda^{-3}(1 - \xi^{-1})^{3/2}, \qquad (15) $$

which we can write in the form (4) if we set

$$ \mathcal{T} \equiv (1 - \xi^{-1})^{-1/2} \lambda / \mathcal{U}. $$

The identification of $\mathcal{D}$ with $\mathcal{D}_2$ has the immediate advantage that the resulting expression for $\mathcal{D}$ tends to zero as $\xi \to 1$. The expression of HL results by taking the limit $\xi \to \infty$.

In a comparison of homogeneous and channel quasigeostrophic models, Pavan and Held (1996) fit the lower-layer PV fluxes in the homogeneous model with a diffusive closure and then use this diffusivity to obtain good fits to numerical integrations for inhomogeneous baroclinically unstable jets. The expression that they obtained is precisely (15), although they were unaware of the argument of the preceding paragraph.

Given the good agreement with (1) whether we use $\varepsilon_d$ or $\varepsilon_p$, and that (14) is exact for our two-layer model, we might expect the Pavan and Held formula (15) to work equally well. However, it is informative, in this regard, to consider the implications of assuming a small departure in the $\mathcal{D}-\varepsilon$ scaling. More precisely, suppose that, for fixed $\beta$,

$$ \mathcal{D} \propto \varepsilon^{3/5+\alpha}, \qquad (16) $$

where $\alpha$ is a small number. Using (11), we find

$$ \mathcal{D} \propto \xi^{\mu}(\xi - 1)^{\mu}; \qquad \mu = \frac{3/5 + \alpha}{2/5 - \alpha}. \qquad (17) $$

Assuming $\alpha$ to be small, the scaling exponent is $3/2 + 25\alpha/4$. Thus, a small imprecision in $\alpha$ is amplified nearly sixfold in the expression for $\mathcal{D}(\xi)$. The implication is that it is much more demanding to obtain an accurate theory for heat flux, or diffusivity, as a function of the mean gradients than a theory for diffusivity as a function of $\varepsilon$.

The numerical results are displayed in Fig. 5a for our new simulations and for those tabulated in HL. For the 256² simulations, the scaling works only marginally and the differences between the various versions of the model are larger than in Fig. 3a, as anticipated from the argument of the preceding paragraph. The 1024² simulations show a better fit with the scaling provided by (15), especially if we discard the larger values of $\xi$ for which we know that the surface drag is stopping the cascade before $\beta$ can play its role. Additionally, the fit
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Nondimensional diffusivity $D_2/U\lambda$ as a function of $c^{5/2}\xi^2(1 - \xi^{-1})^{3/2}$. The circles are our 256² simulations, the crosses are HL simulations, and the stars are the 1024² simulations. The solid line is the perfect fit. (b) Nondimensional diffusivity $D_2/U\lambda$ as a function of $c^{5/2}\xi^2(1 - \xi^{-1})^{3/2}$ (stars) and $c^{5/2}\xi^2$ (plus signs) for the 1024² simulations. The value of c is the same as in Fig. 3.</figcaption></figure>

is much better than the fit with the theory of HL that uses only the large supercriticality limit of (15), as can be seen in Fig. 5b. In light of the various reasons to expect some discrepancies, we are encouraged by the closeness of the fit, especially in the 1024² simulations.

Keeping in mind how one might generalize our results to a vertically continuous model, one can start by considering a diffusivity of the form $D(z) = \gamma(z)D_0$, with $D_0$ determined as above at some level $z_0$, where $\gamma(z_0) = 1$. In a quasigeostrophic model with mean vertical shear $U(z)$, the energy production is set by the appropriate generalization to (13). One obtains

$$D_0 = \beta^{-2} \left[ \frac{1}{H} \int U(z) \gamma(z) \partial_y Q(z) dz \right]^{3/2}$$

as a generalization of (15). We thus have a generalized timescale

$$T = \left[ \frac{1}{H} \int U(z) \gamma(z) \partial_y Q(z) dz \right]^{-1/2}.$$

The resulting timescale depends on how the vertical shear projects on the meridional PV gradient. (If the rhs of this expression is negative, the implication would be that one cannot diffuse PV downgradient with this vertical structure in the diffusivity and simultaneously generate eddy energy, so the diffusivity should be set to zero.) Given the distribution of $\gamma(z)$, we can determine $D_0$ and, therefore, $D$ at each level.

b. Energy production and entropy

Rather than using the generation of available potential energy or the destruction of barotropic energy by large-scale friction, Barry et al. (2002) estimate $\varepsilon$ by computing the entropy destruction by the mean radiative

and boundary fluxes and assuming that a fixed fraction $\eta$ of this entropy destruction is balanced by the entropy generated when kinetic energy is dissipated:

$$\frac{\varepsilon}{T_0} = \eta \left[ \frac{Q}{T} \right]. \qquad (18)$$

Here, $T_0$ is the average temperature at which kinetic energy is dissipated, we interpret $Q$ as the divergence of the large-scale eddy flux of energy, and brackets denote a global mean. One then has approximately

$$\varepsilon \propto Q_0 \frac{\delta T}{T_0}, \qquad (19)$$

where $\delta T$ is a characteristic temperature difference between the heated and cooled regions and $Q_0$ the magnitude of the corresponding characteristic heating and cooling rates.

There is a close relationship between estimates of kinetic energy dissipation based on the entropy and available potential energy balances. Indeed, Lorenz's original motivation for introducing the concept of available potential energy was to provide a more robust way of estimating $\varepsilon$ from observations than that provided by expressions along the lines of (18) (see Lorenz 1967).

In our simple Boussinesq two-layer model, there is no distinction between entropy production and the generation of available potential energy. Consider a two-layer flow forced by an imposed mass exchange between the two layers, with lower-layer mass converted to upper-layer mass in low latitudes to represent heating and an equal amount of upper layer converted to lower layer in high latitudes to represent cooling. The rate of generation of available potential energy per unit mass of the fluid by this process in the quasigeostrophic limit is
$$ \varepsilon = \frac{g^*}{2H_0} [Qh], \qquad (20) $$

where $Q$ is the mass source per unit time in the lower layer, $h$ is the lower-layer depth, $g^*$ is the reduced gravity, and $H_0$ is the mean depth of either layer. This expression is analogous to (18) with $h$ playing the role of $1/T$. If there were significant thermal damping resulting in the loss of available potential energy, an effect absent in our model and ignored in our analysis above, this would result in a reduction in “efficiency,” analogous to the factor $\eta$ in (18).

Assume for simplicity that eddy momentum fluxes are negligible in the time mean. Then, in the quasigeostrophic approximation, there is no Eulerian mean meridional flow and $Q$ must be balanced by the eddy thickness flux:

$$ Q = \frac{\partial}{\partial y} \overline{v' h'}. \qquad (21) $$

Integrating (20) by parts one obtains

$$ \varepsilon = - \frac{g^*}{2H_0} \left[ \frac{\partial h}{\partial y} \overline{v' h'} \right]. \qquad (22) $$

Using the thermal wind equation to relate $g^* \partial h / \partial y$ to the mean vertical shear, and setting the eddy PV flux equal to $(f/H_0) \overline{v' q'}$, one recovers (14). The constraint between $\varepsilon$ and $\mathcal{D}$ provided by the global mean entropy or available potential energy budgets looks like it is fundamentally nonlocal and inhomogeneous. But this impression is created by the integration by parts. To the extent that eddies are produced and dissipated at more or less the same latitude, (22) or (13) can be satisfied at each latitude, the homogeneous model should capture the essence of the eddy maintenance, and we can think of an approximate diffusivity that is a function of the local environment.

## 4. Discussion and conclusions

The relation (1) between eddy diffusivity and energy generation rate is potentially of considerable significance for our understanding of atmospheric circulation, as is implicit in Barry et al. (2002). We have tried to clarify the role of this constraint in the theory of Held and Larichev (1996). The latter is a closure theory for the eddy buoyancy and PV fluxes in a horizontally homogeneous, baroclinically unstable, two-layer model. We first verify that this expression is indeed satisfied by numerical solutions to this homogeneous model in the weakly unstable limit $\xi < 5$. A three-layer model also confirms our findings. The reason why (1) works reasonably well is that it is not very sensitive to the exact shape of the barotropic energy spectrum: the scaling works much better in our case than the $k^{-5/3}$ power-law prediction. We cannot rule out the possibility that there are other ways of motivating (1) without referring

to the inverse energy cascade at all. Our theory relies only on the existence of an upscale cascade for the barotropic energy associated with a significant part of energy dissipation at large scales, and on the role played by $\beta$ in the inverse cascade.

The key ingredient in the closure of HL is the identification of the diffusivity in (1) with the diffusivity computed from an appropriate flux in the homogeneous model and the associated mean gradient. In the large supercriticality limit considered in HL, the diffusivities for buoyancy, upper-layer PV, and lower-layer PV are all equal, so one does not have to choose among these, but the resulting theory does not fit results for moderate supercriticality very well. Here, we demonstrate that the simple choice of identifying the diffusivity in (1) with the diffusivity for lower-layer PV results in a good fit with numerical simulations, especially the 1024² simulations. In all cases, the fit is superior to the large supercriticality limit of HL. The resulting theory has, in fact, already been suggested by Pavan and Held (1996), based not on a “derivation” of the kind described here, but rather on a simple fit to numerical results.

We do not have a clear understanding of why selecting lower-layer diffusivity for the theory is appropriate. One possibility discussed above is that this theory is most relevant for levels close to the typical steering level of the dominant waves. An alternative heuristic viewpoint is as follows. An aspect of the theory of HL, following Rhines (1979) and Salmon (1980), is that a substantial inverse energy cascade results in the energy-containing eddies being dominated by the barotropic mode. However, if the supercriticality of the flow is not very large, there is not enough of an inverse cascade to create a fully barotropic flow. Instead, the eddy kinetic energy tends to be concentrated in the upper layer, analogous to the upper-tropospheric maximum in eddy kinetic energy observed in the atmosphere. As long as the horizontal scales are comparable in the two layers, the eddy potential entrophy will also be larger in the upper than in the lower layer. One can always think of the flow in a two-layer quasigeostrophic model as the linear superposition of a part induced by the upper-layer PV distribution and a part induced by the lower-layer PV distribution. If the ratio of potential entrophy in the upper to the lower layer is large enough, then the flow in the lower layer can be predominantly induced by the PV in the upper layer. As a consequence, the lower-layer PV will behave like a passive tracer. In the upper layer, the eddies are advected by a flow that they themselves induce to a great extent, so arguments based on passive advection are less directly relevant. This upper/lower-layer asymmetry is ultimately due to the drag in the bottom layer and the weaker PV gradient there.

The two-layer homogeneous model has the distinctive feature that the upper-layer PV flux is known once the lower-layer flux is known—they are simply equal in magnitude and opposite in sign. Since the generation of eddy energy is determined by the PV fluxes, one has
another relation between $\epsilon$, the lower-layer PV diffu-
sivity and the mean PV gradients that one can couple
with (1) to determine the energy throughput and the
diffusivity separately. In a model with more than two
levels in the vertical, this argument is no longer closed,
as one needs a theory for the vertical distribution of the
PV flux. See Held and O'Brien (1992), Smith and Vallis
(2002), and Tréguier et al. (1997) for a discussion of
this issue, which we consider an unsolved problem.

*Acknowledgments.* The authors would like to thank Shafer Smith, Paul Kushner, and Geoff Vallis for helpful comments about this work. All simulations were performed on the supercomputers at NOAA's Geophysical Fluid Dynamics Laboratory, using a quasigeostrophic code originally written by S. Smith and G. Vallis and parallelized by GL and S. Smith. GL was funded under a grant/cooperative agreement with the National Oce- anic and Atmospheric Administration (NA07RG0002).

APPENDIX

Numerical Model

We have reexamined integrations similar to those in
HL, partly as a point of departure for the moist model
described in Lapeyre and Held (2003, unpublished man-
uscript). We take the standard horizontally homoge-
neous model on the $\beta$ plane with a zonally symmetric
climate and layers of equal mean depths, as in HL. The
mean zonal velocity in each layer is $U_1 = U$ and $U_2 = -U$ where the subscripts 1 and 2 refer respectively to
the top and bottom layers of fluid. Each layer has equal
mean depth $H_1 = H_2 = H$. The supercriticality $\xi \equiv (U_1 - U_2)/(\beta\lambda^2)$ (where $\lambda$ is the radius of deformation) must
be larger than unity for instability.

We examined different sets of simulations, using both
256² resolution as in HL, but also the higher resolution
of 1024². The current model uses an exponential filter
to dissipate small scales [see Smith et al. (2002) for a
description of the scheme], rather than the $\nabla^8$ damping
used in HL. This modifies the energetics of the system
as we now find typically that less than 3% of the energy
input rate $\epsilon_p$ is dissipated by small-scale diffusion, leav-
ing 97% to be destroyed by Ekman friction (the anal-
ogous figures were only 30% and 70% in HL). As a
consequence our 256² and 1024² simulations are sub-
stantially more energetic than in HL, holding other pa-
rameters fixed. Additionally, we vary the supercriticality
$\xi$ in our 256² simulations by varying $\lambda$, rather than $\beta$
as in HL, motivated by the desire to test if the depen-
dence of the dry model results on static stability can be
used to help understand our companion moist simula-
tions. The largest wavelength in the domain ranges from
$20\pi\lambda$ for the smaller supercriticality to $34\pi\lambda$ for the
largest. This leaves more room for the direct cascade of
enstrophy compared to the simulations of HL where the
largest scale was chosen to be $100\pi\lambda$. For our 1024²

simulations, we chose 100$\pi\lambda$ for the largest scale as in
HL, but we used a larger bottom friction $r = 0.32 U\lambda^{-1}$
instead of $r = 0.16 U\lambda^{-1}$ for the 256² simulations here
and in HL. In summary, our new simulations are less
viscous than those in HL, and, in conjunction with HL,
provide some information on the robustness of the scal-
ing to the size of the domain and to the strength of
surface friction. The range of supercriticality $\xi$ we use
in these different simulations is between 1.1 and 14 in
order to push the limits of the theory to the weakly
unstable limit.

REFERENCES

Barry, L., G. C. Craig, and J. Thuburn, 2002: Poleward heat transport by the atmospheric heat engine. *Nature*, **415**, 774–777.

Danibu, S., and D. Gurarie, 2001: Nonuniversal features of forced two-dimensional turbulence in the energy range. *Phys. Rev.*, **63E**, R020203.

—, and —, 2002: Rhines scale and spectra of the beta-plane turbulence with bottom drag. *Phys. Rev.*, **65E**, R067301.

Haynes, P. H., J. F. Scinocca, and M. D. Greenslade, 2001: Formation and maintenance of the extratropical tropopause by baroclinic eddies. *Geophys. Res. Lett.*, **28**, 4179–4182.

Held, I. M., and E. O'Brien, 1992: Quasigeostrophic turbulence in a three-layer model: Effects of vertical structure in the mean shear. *J. Atmos. Sci.*, **49**, 1861–1870.

—, and V. D. Larichev, 1996: A scaling theory for horizontally homogeneous, baroclinically unstable flow on a beta plane. *J. Atmos. Sci.*, **53**, 946–952.

Lorenz, E. N., 1967: *The Nature and Theory of the General Circulation of the Atmosphere*. WMO, 161 pp.

Maltrud, M. E., and G. K. Vallis, 1991: Energy spectra and coherent structures in forced two-dimensional and beta-plane turbulence. *J. Fluid Mech.*, **228**, 321–342.

Marshall, D. P., R. G. Williams, and M.-M. Lee, 1999: The relation between eddy-induced transport and isopycnic gradients of potential vorticity. *J. Phys. Oceanogr.*, **29**, 1571–1578.

Panetta, R. L., 1993: Zonal jets in wide baroclinically unstable regions: Persistence and scale selection. *J. Atmos. Sci.*, **50**, 2073–2106.

Pavan, V., and I. M. Held, 1996: The diffusive approximation for eddy fluxes in baroclinically unstable jets. *J. Atmos. Sci.*, **53**, 1262–1272.

Rhines, P. B., 1975: Waves and turbulence on a $\beta$-plane. *J. Fluid Mech.*, **69**, 417–443.

—, 1979: Geostrophic turbulence. *Annu. Rev. Fluid Mech.*, **11**, 404–441.

Salmon, R. S., 1980: Baroclinic instability and geostrophic turbulence. *Geophys. Astrophys. Fluid Dyn.*, **15**, 167–211.

Smith, K. S., and G. K. Vallis, 2002: The scales and equilibration of midocean eddies: Forced-dissipative flows. *J. Phys. Oceanogr.*, **32**, 1699–1721.

—, G. Boccaletti, C. C. Henning, I. N. Marinov, C. Y. Tam, I. M. Held, and G. K. Vallis, 2002: Turbulent diffusion in the geostrophic inverse cascade. *J. Fluid Mech.*, **469**, 14–47.

Spall, M. A., 2000: Generation of strong mesoscale eddies by weak ocean gyres. *J. Mar. Res.*, **58**, 97–116.

Tréguier, A. M., I. M. Held, and V. D. Larichev, 1997: Parameterization of quasigeostrophic eddies in primitive equation ocean models. *J. Phys. Oceanogr.*, **27**, 567–580.

Vallis, G. K., 1988: Numerical studies of eddy transport properties in eddy-resolving and parameterized models. *Quart. J. Roy. Meteor. Soc.*, **114**, 183–204.

—, and M. E. Maltrud, 1993: Generation of mean flows and jets on a beta plane and over topography. *J. Phys. Oceanogr.*, **23**, 1346–1362.