## Abstract

Mesoscale eddies play a major role in the transport of tracers in the ocean. Focusing on a sector in the east Pacific, the authors present estimates of eddy diffusivities derived from kinematic tracer simulations using satellite-observed velocity fields. Meridional diffusivities are diagnosed, and how they are related to eddy properties through the mixing length formulation of Ferrari and Nikurashin, which accounts for the suppression of diffusivity due to eddy propagation relative to the mean flow, is shown. The uniqueness of this study is that, through systematically varying the zonal-mean flow, a hypothetical “unsuppressed” diffusivity is diagnosed. At a given latitude, the unsuppressed diffusivity occurs when the zonal-mean flow equals the eddy phase speed. This provides an independent estimate of eddy phase propagation, which agrees well with theoretical arguments. It is also shown that the unsuppressed diffusivity is predicted very well by classical mixing length theory, that is, that it is proportional to the rms eddy velocity times the observed eddy size, with a spatially constant mixing efficiency of 0.35. Then, the suppression factor is estimated and it is shown that it too can be understood quantitatively in terms of easily observed mean flow properties. The authors then extrapolate from these sector experiments to the global scale, making predictions for the global surface eddy diffusivity. Together with a prognostic equation for eddy kinetic energy and a theory explaining observed eddy sizes, these concepts could potentially be used in a closure for eddy diffusivities in coarse-resolution ocean climate models.

## 1. Introduction

Fluctuations on scales of roughly 20–300 km, which derive their energy primarily from the baroclinic instability of the large-scale density field, pervade the global ocean and contain a large fraction of the ocean’s energy (Gill et al. 1974). These fluctuations, known as mesoscale eddies, dominate the dispersion of particles and the mixing of tracers on large space and time scales (Lumpkin and Elipot 2010). Understanding the mixing induced by mesoscale eddies is a problem of both fundamental theoretical interest and practical importance.

The practical importance arises in coarse-resolution ocean climate models, which do not resolve mesoscales. Such models simulate the transport of tracers using a Reynolds-averaged formulation of the tracer conservation equation: where is the average concentration of the tracer in a grid box, is the resolved velocity field, **u**′ is the eddy velocity, and is an unresolved, subgrid-scale “eddy” flux, which must be parameterized. By far the most common approach is to use a closure of the form . To isolate the physical processes underlying the eddy flux, the tensor , which relates the three-dimensional eddy fluxes to the local background gradients, is usually expressed in terms of an advective component and a diffusive component, and the diffusive component is rotated into an isopycnal component (i.e., along density surfaces) and a diapycnal component (i.e., across density surfaces) (Redi 1982; Griffies et al. 1998). The diapycnal part of , related to breaking internal waves and other small-scale processes, is an active topic of research but is not the focus here. The isopycnal part of represents lateral mixing by mesoscale eddies.

Ocean climate models originally specified the components of as constants, but ample research based on both observations and eddy-resolving models has revealed strong spatial variability in mesoscale eddy mixing rates (Holloway 1986; Stammer 1998; Zhang et al. 2001; Zhurbas and Oh 2003, 2004; Ferreira et al. 2005; Marshall et al. 2006; Rypina et al. 2012; Fox-Kemper et al. 2013; Abernathey and Marshall 2013). Models that have incorporated spatial variation in mixing parameters have been able to reduce their biases, but there is still much room for improvement (Visbeck et al. 1997; Ferreira et al. 2005; Danabasoglu and Marshall 2007). Arguably, the two greatest obstacles limiting progress are 1) poor knowledge of the eddy diffusivities in the real ocean and 2) lack of understanding of how these diffusivities are related to the large-scale fields.

In this study, we demonstrate that eddy diffusivities derived from satellite observations are consistent with a simple formula, whose ingredients are the root-mean-square (rms) eddy velocity *u*_{rms}, which is proportional to the eddy kinetic energy (EKE), the eddy size, the eddy propagation speed, and the large-scale mean flow speed. Together with a prognostic equation for EKE (Eden and Greatbatch 2008; Marshall and Adcroft 2010) and a theory for the eddy sizes, these results can then be used to develop a complete eddy closure. This study is focused on relating observable eddy properties to diagnosed diffusivities; a related study (Bates et al. 2013, manuscript submitted to *J. Phys. Oceanogr.*) describes how to implement these insights in a complete eddy parameterization.

The central tool in our study is the 20-year record of sea surface height (SSH) observed by satellite altimetry. These data provide both time-mean statistical properties of eddies and a continuous record of surface geostrophic velocity fields. Past authors have used statistical information from the altimeter to estimate eddy diffusivity but lacked a “ground truth” based on an empirical observation of the eddy flux against which to compare their predictions (e.g., Holloway 1986; Stammer 1998). More recent studies have used the altimetric velocities to directly simulate the transport of passive tracers, leading to data-derived estimates of (Marshall et al. 2006; Ferrari and Nikurashin 2010; Abernathey and Marshall 2013). We follow this approach here. For simplicity we first focus on a sector in the east Pacific with symmetry in longitude but strong variations of eddy and mean flow properties with latitude (see Fig. 1). We calculate the meridional diffusivity *K*_{obs} as a function of latitude for a wide range of zonal-mean flows. The data and numerical experiments are detailed in section 2. We then use these experiments to validate theoretical predictions for the diffusivity based on mixing length theory. The theoretical framework is outlined in section 3. Section 4 contains a discussion of the difference between nonlinear eddies and linear Rossby waves, concluding that our framework is not applicable to the wavelike equatorial region. In section 5, we infer the phase speeds of eddies from our simulations, and in section 6 we examine the length scales associated with unsuppressed diffusivities. Section 7 calculates the suppression factor for a realistic zonal-mean flow. Finally, in section 8 we extrapolate our results to the global ocean and discuss predicted eddy diffusivities. A concluding discussion is undertaken in section 9.

## 2. Satellite data and diffusivity calculations

Obtaining direct estimates of mesoscale eddy diffusivity from observations is challenging. The most straightforward method is to construct a local diffusivity based on mooring measurements of the eddy flux of heat (Bryden and Heath 1985). This approach is compromised by the shortness of the available time series (Wunsch 1999) and by the presence of “rotational fluxes,” nondivergent components of the eddy flux vector field that can dominate local measurements (Marshall and Shutts 1981). Lagrangian observations, from surface drifters and subsurface floats, can also be used to estimate diffusivity [review by LaCasce (2008) and Rypina et al. (2012)]. This approach has been more successful but is still limited in accuracy and spatial resolution by the number and location of drifter trajectories (Klocker et al. 2012a,b). In recent years, an efficient, accurate method has been developed that uses real surface velocities (as observed by satellite) in combination with simulated passive tracers to infer eddy diffusivities (Marshall et al. 2006; Abernathey et al. 2010; Ferrari and Nikurashin 2010; Klocker et al. 2012b; Abernathey and Marshall 2013). The tracer-based approach offers global coverage for the entire satellite era, which began in 1992.

Satellite observations of sea surface height anomaly provide weekly geostrophic velocities at the sea surface and resolve the largest mesoscale eddies (Chelton et al. 2007). Eddy properties vary widely across the globe. To simplify our problem somewhat, we focus on a sector in the east Pacific between 180° and 130°W that is mostly free of land and in which eddy properties are relatively homogeneous with longitude. This allows us to focus on the variation of eddy properties and mixing rates as a function of latitude only. A snapshot of the sea surface height anomaly in this sector is shown in Fig. 1, revealing the meridional variations in eddy propagation speed and intensity. This sector was also analyzed by Abernathey and Marshall (2013), and our study builds on that work.

To calculate meridional eddy diffusivities, we simulate the advection of a passive tracer using satellite-derived velocity fields. The velocity dataset we employ is the surface geostrophic velocity anomaly from Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO). We use 17 years worth of observations, beginning with 6 January 1993. In addition to the standard geostrophic balance, this dataset employs the empirically validated “equatorial–geostrophic” approximation of Lagerloef et al. (1999) to calculate velocities near the equator (between 5°N and 5°S). The AVISO velocity anomalies are interpolated to a 1/10° grid, and a small correction is applied to remove divergence and to enforce periodicity in longitude. For further details of the data and processing, the reader is referred to Abernathey and Marshall (2013). In addition to the AVISO velocity anomalies, different zonal-mean flows are superimposed, as described in the following section.

Given a velocity field, we then solve the advection–diffusion equation for a passive tracer *C* using the Massachusetts Institute of Technology General Circulation Model (MITgcm) in offline mode (Marshall et al. 1997). The initial condition for the passive tracer is simply *C* = *ϕ*; that is, the tracer is proportional to latitude. Snapshots of the tracer field after 3 months of advection are shown in Fig. 2. The tracer is reset to the initial condition each year. The eddy stirring produces a meridional flux of tracer ; the overbar indicates an average taken zonally (eliminating rotational fluxes), in time over each year, and over an ensemble of 17 years. These calculations provide our empirical estimate of diffusivity in the meridional direction:

Extensive validation of this method has been performed in previous studies, confirming the insensitivity of *K*_{obs} to the averaging period (provided it exceeds a few months) and the numerical parameters (Marshall et al. 2006; Abernathey and Marshall 2013). It was also confirmed that that *K*_{obs} agrees very well with estimates based on Lagrangian trajectories (Klocker et al. 2012a,b) and with the “effective diffusivity” of Nakamura (1996). In what follows, we describe a series of experiments designed to investigate how *K*_{obs} is related to the properties of the flow. But first we review the necessary theoretical background.

## 3. Mixing length suppression

### Mixing length theory

We seek to relate Eulerian eddy statistics to the Eulerian flux gradient diffusivity in Eq. (1). (The relationship with the Lagrangian framework is discussed in appendix A.) A common approach in this context is to express the eddy diffusivity in terms of mixing length arguments (Taylor 1915; Prandtl 1925), which state that diffusivity is proportional to the rms eddy velocity, where , times a mixing length *L*_{mix}:

Here, Γ is a mixing efficiency usually assumed to be an order-one constant. The mixing length is conceptually analogous to the concept of the mean free path in thermodynamics: a fluid parcel will conserve its properties for a characteristic length *L*_{mix} before mixing with the surrounding fluid. From the mixing length perspective, a self-contained closure theory for eddy diffusivity should predict *u*_{rms} and *L*_{mix} based on the properties of the large-scale flow. This approach underlies many proposed theoretical closures for geostrophic turbulence and baroclinic equilibration, in both the ocean and atmosphere (Bretherton 1966; Green 1970; Stone 1972; Held and Larichev 1996; Eden and Greatbatch 2008; Marshall and Adcroft 2010; among many). It has also been used to estimate eddy diffusivities from observations (Holloway 1986; Stammer 1998). The observational problem is somewhat easier because *u*_{rms} is readily observed by satellite or drifters; most of the theoretical difficulty therefore lies in specifying *L*_{mix} and Γ. In the works cited above, *L*_{mix} is often assumed to be the length scale of the largest eddies. However, recent developments (and this study) argue that this is not necessarily the case; the mixing length can in fact be much smaller than the eddy scale, for reasons discussed in the next paragraph. Note that the mixing length expression in Eq. (2) is different from mixing length methods based on tracer variance, which have been used in some observational studies (Armi and Stommel 1983; Naveira Garabato et al. 2011) and that are briefly described in appendix B.

Recent work based on both observations and theoretical arguments has shown that the presence of eddy propagation relative to a background mean flow suppresses eddy diffusivities in the across-current direction (Marshall et al. 2006; Ferrari and Nikurashin 2010; Klocker et al. 2012a,b; Abernathey and Marshall 2013; Tulloch et al. 2013, manuscript submitted to *J. Phys. Oceanogr.*). This phenomenon can be understood heuristically as follows: if there is no mean flow, and if the eddies are stationary, the eddy stirring acts coherently on the same water masses for a long time, efficiently mixing their properties. But if an individual eddy moves relative to the underlying water, then it no longer acts on the same water masses, and consequently mixing is less efficient. Because *u*_{rms} is independent of mean flow or eddy propagation, this suppression effect must be interpreted as the result of a reduced mixing length.

An analytical model of this phenomenon was developed by Ferrari and Nikurashin (2010), who showed that, for a weakly nonlinear isotropic eddy field, the mixing length can be written as

where *γ*^{−1} is the eddy decorrelation time, *k* is the eddy wavenumber, *c*_{w} is the absolute eddy phase speed,^{1 }*U* is the background zonal velocity (depth dependent in general, but here always the value at the surface). The term *L* is the mixing length for isotropic turbulence, which we assume to be the dominant Eulerian length scale of the eddies, such that *L* = 2*π*/*k*. [It will be shown subsequently that this assumption works well. The relationship with the Lagrangian formulation of Ferrari and Nikurashin (2010) is discussed in appendix A.] If *c*_{w} = *U*, no suppression occurs, and Eq. (3) reduces to *L*_{mix} = *L*; that is, the mixing length depends directly on the eddy size, and the diffusivity is maximized. We call the corresponding “unsuppressed” diffusivity

But if *c*_{w} ≠ *U*, then *L*_{mix} < *L*, leading to a smaller *K*. Note that Eq. (3) relies on several assumptions. The two major assumptions are that (i) there is a scale separation between the eddies and mean flow, and (ii) the flow is a parallel shear flow. These assumptions likely fail in certain parts of the ocean where strong currents curve significantly, for example, in the western boundary currents or downstream of major topographic features in the Antarctic Circumpolar Current (ACC) (Naveira Garabato et al. 2011). The latter case will be discussed in detail in section 7. Meanwhile, we will focus on the Pacific sector, where no such issues arise, and, as we shall demonstrate, Eq. (3) is quite accurate. The symbols used in the paper are summarized in Table 1.

This new understanding of mixing length suppression represents a significant advancement in the theory of mesoscale mixing, but so far the issue has been examined only in the Southern Ocean. Furthermore, Eqs. (2) and (3) are only useful for coarse-resolution climate models if the parameters *u*_{rms}, *L*, *c*_{w}, and *γ* can be specified easily based on large-scale fields. Moreover, many simplifying assumptions underlie the theories explained above, and their validity must be tested experimentally. In this study, we test these ideas and their potential limitations in our Pacific sector (Fig. 1) by systematically varying an imposed zonal-mean flow and examining the resulting dependence of *K*_{obs} at each latitude. This is a purely kinematic exercise, because the flow fields are not dynamically consistent in general (Haynes et al. 2007).

We run a series of simulations where the tracer *C* is advected by the velocity field **u** = **u**′ + *U*_{0}**i**, where **u**′ is the time-dependent, AVISO-derived eddy velocity; *U*_{0} is a constant zonal-mean velocity independent of latitude; and **i** is the zonal unit vector. The set of 60 simulations employs values of *U*_{0} between −0.2 and 0.5 m s^{−1}. Meridional diffusivity *K*_{obs} [Eq. (1)] is calculated for each *U*_{0} experiment as a function of latitude. The resulting diffusivities are shown as gray shading in Fig. 3a, with examples for particular latitudes plotted in Fig. 3b. First, this shows unequivocally that the diffusivity is highly dependent on the mean flow. Further, as expected from Eqs. (2) and (3), for some value of mean flow (different at each latitude), *K*_{obs} reaches a maximum. We identify this maximum diffusivity with the unsuppressed eddy diffusivity . The diffusivities drop off for smaller and larger values of *U*_{0}. The skill of Eq. (3) is evident in the shape of the *K*_{obs}(*U*_{0}) curves, which match very closely with the Lorentzian functional form predicted by the equation. (The skill breaks down closer to the equator, for reasons discussed below.) Furthermore, this method not only identifies *K*^{0}, the unsuppressed diffusivity, but also effectively provides an independent estimate of *c*_{w}, the eddy propagation speed, as the value of *U* that maximizes the diffusivity at each latitude.

## 4. Linear Rossby waves versus nonlinear mesoscale eddies

Before interpreting our results, it is helpful to distinguish Rossby waves from nonlinear mesoscale eddies. Linear Rossby waves are zonally propagating oscillations that arise due to the presence of a planetary vorticity gradient. Mesoscale eddies are strongly nonlinear turbulent fluctuations that derive their energy primarily from baroclinic instability. In relation to tracer transport, the most important difference between Rossby waves and nonlinear mesoscale eddies is that eddies can cause irreversible mixing of tracers, whereas linear Rossby waves cannot; distinguishing regions dominated by linear Rossby waves from those of nonlinear mesoscale eddies (or geostrophic turbulence) is therefore important for the interpretation of the results presented in this paper. The analytical model for the eddy diffusivity discussed above [e.g., Eqs. (2) and (3)] applies only to nonlinear eddies. In principle waves and nonlinear eddies of different scales can coexist, given the broad spectrum of variability in the ocean (Wunsch 2010). However, it is evident that at low latitudes, eddy activity is suppressed and energy is transferred instead into waves, while at mid- and high latitudes, eddies are much more prevalent.

Several different approaches have recently been employed to distinguish regions dominated by linear Rossby waves from more nonlinear regions. One such approach is to use the nonlinearity parameter *u*_{rms}/*c* (Chelton et al. 2007, 2011; Early et al. 2011). (Here *c*, the intrinsic phase speed, differs from *c*_{w}, the absolute phase speed; *c* does not include the Doppler shift by the depth-averaged mean flow.) If *u*_{rms}/*c* > 1, that is, if the rotational eddy velocity *u*_{rms} exceeds its translational speed *c*, transforming the coordinates into a co-moving frame will lead to closed streamlines within the eddy. This leads to an eddy with an inner core, which traps fluid and advects it along its path and an outer ring, which continuously entrains and sheds fluid (Early et al. 2011). This outer ring mixes fluid over the eddy size *L* with the rotational velocity *u*_{rms}, leading to the mixing length arguments in Eq. (2). If *u*_{rms}/*c* < 1, contours are not closed and the anomalies are more wavelike. Theiss (2004) instead framed his analysis in terms of the Rhines scale *L*_{Rh} and the deformation radius *L*_{D}; he argued that when *L*_{Rh}/*L*_{D} > 1, waves are not able to transfer energy into alternating zonal flows, permitting isotropic eddies to exist. This argument implies a critical latitude at *L*_{Rh}/*L*_{D} = 1. A similar transition from waves to turbulence was demonstrated by Tulloch et al. (2009), who determined the wavelength required to optimally fit the phase speed predicted by linear theory to that observed by altimetry. They concluded that at low latitude there is an overlap between geostrophic turbulence and Rossby wave time scales due to an upscale energy transfer that produces waves, whereas at high latitudes there is no such overlap; that is, waves were not produced and eddies can exist.

Inspired by the studies cited above, we attempt to characterize wavelike versus turbulent regions using the ratio

that is, the ratio between the rotational eddy velocity *u*_{rms} and the eddy phase speed (or translational velocity) *c*. The Rhines scale is defined as , while *c* is related to *L*_{D} through the long-wave Rossby-wave dispersion relation . The second equality demonstrates the equivalence between the approach of Chelton et al. (2011) and that of Theiss (2004). For regions where *r* < 1, we therefore expect an upscale energy transfer to produce waves, whereas for regions with *r* > 1, the gap in time scales between Rossby waves and turbulence supports the existence of nonlinear mesoscale eddies. The critical latitude is at *r* = 1.

The ratio *r* for the Pacific patch used in this study is shown in Fig. 4a, demonstrating a transition from waves (*r* < 1) to eddies (*r* > 1) equatorward of the critical latitude at approximately 18°. Figure 4b shows the number of eddies found per square degree latitude by an eddy-tracking algorithm (Chelton et al. 2011), where an eddy is defined by a closed sea surface height contour. Equatorward of the critical latitude there are almost no eddies observed; presumably instabilities at these low latitudes transfer their energy into alternating zonal flows or waves, rather than producing closed eddies. Figure 4c shows the deformation scale *L*_{D} and the Rhines scale *L*_{Rh}, which are equal at the critical latitude. One consequence of this transition from waves to turbulence is that we expect different mixing behavior in the two regions. In particular, we do not expect the model of Ferrari and Nikurashin (2010) [Eq. (3)] to be valid in the more linear equatorial region.

## 5. Eddy phase speeds

Past studies have analyzed the propagation of sea surface height anomalies in the altimetric record (Chelton and Schlax 1996; Chelton et al. 2007, 2011). In recent years, there has been an effort to reconcile such observational estimates with theoretical predictions, and a debate has arisen over what dynamics are observed in altimetric SSH anomaly signals—linear Rossby waves or nonlinear mesoscale eddies (Chelton and Schlax 1996; Killworth et al. 1997, 2004; Chelton et al. 2007; Tulloch et al. 2009; Chelton et al. 2011; O'Brien et al. 2013)? Our experiments can shed light on this debate by providing an independent measurement of the phase speed through an inversion of Eqs. (2) and (3). In the experiments with a range of constant *U*_{0} (Fig. 3), we found a maximum diffusivity at each latitude. Equation (3) says that the value of *U*_{0} that maximizes *K*_{obs} is equal to *c*_{w}, the Doppler-shifted eddy phase speed; this maximum value is indicated by the blue line in Fig. 3a. Furthermore, because linear waves do not cause mixing, we know that the value of *c*_{w} estimated in this way is associated with nonlinear eddies.

Nonlinear mesoscale eddies are expected to propagate at the speed predicted by the linear baroclinic Rossby wave theory (McWilliams and Flierl 1979; Killworth 1986; Cushman-Roisin et al. 1990). In the long Rossby wave limit (|**k**| → 0), and taking into account Doppler shifting by the depth-averaged mean flow (Klocker and Marshall 2013, manuscript submitted to *Geophys. Res. Lett.*), this speed can be written as

To evaluate this equation for the Pacific region, we borrowed an estimate of *L*_{D} from Tulloch et al. (2009), who calculated it from the mean stratification of an ocean state estimate (Forget 2010). (The term *L*_{D} itself is plotted in Fig. 4c.) The depth-averaged mean flow component is calculated using the depth-averaged mean flow of the same hydrography. In Fig. 3a, the phase speed estimated using the inversion of eddy diffusivity estimates (blue line) is compared to the phase speed from Eq. (6) (red line). These two values agree very well poleward of the critical latitude. This good agreement therefore confirms that nonlinear mesoscale eddies, which we expect to dominate the ocean poleward of the critical latitude, do travel at the linear Rossby wave speed Doppler shifted by the depth-averaged mean flow. This is encouraging from the perspective of parameterizing *K*_{obs}, because *c*_{w} is straightforward to calculate given the large-scale density field. In section 7, we will use Eq. (6) in Eq. (3) to quantitatively estimate the mixing length suppression effect.

## 6. Eddy length scales

We now turn to understanding the eddy length scales relevant for mixing. To achieve this without the distorting effect of eddy propagation, we use estimates of the unsuppressed eddy diffusivities , for which Eqs. (2) and (3) reduce to . Figure 4 shows the empirically calculated unsuppressed mixing length *L*, the observed eddy size *L*_{obs}, and the Rossby radius of deformation *L*_{D}. For these calculations we use a constant mixing efficiency parameter of Γ = 0.35. This gives the best fit between *L* and *L*_{obs} and also has a physical rationale explained below. The observed eddy size *L*_{obs} is estimated using an eddy-tracking algorithm and defined as the radius of a circle with an area equal to that enclosed by the contour of SSH within the eddy around which the circum-average speed is maximum (Chelton et al. 2011). For regions of nonlinear mesoscale eddies (*r* > 1), the observed eddy size is slightly larger than the Rossby radius of deformation. This is consistent with a weak, upscale transfer of kinetic energy from a source with scales near the Rossby radius of deformation; such an inverse energy cascade is expected from geostrophic turbulence theory (Rhines 1975; Scott and Wang 2005).

The agreement between the observed eddy size *L*_{obs} and the empirically calculated eddy size *L* (using Γ = 0.35) can be seen even more clearly in Fig. 5. This figure shows how eddy length scales depend on the ratio *r*. Figure 5a plots the ratio *l*/*L*_{D} versus *r*, where *l* is either *L*_{obs} (blue dots) or *L* (red dots). Figure 5b plots *L*_{obs} and *L* directly against *r*. Both Figs. 5a and 5b unequivocally show that *L* agrees very well with *L*_{obs} in regions of *r* > 1, with more complicated behavior in regions of *r* < 1. In the past, the relationship between *L* and observed eddy size was more obscure, because previous studies tried to relate suppressed mixing lengths *L*_{mix} to the eddy length scales (e.g., Lumpkin et al. 2002; Eden and Greatbatch 2008; Marshall and Adcroft 2010). Only if the suppression effects can be removed, does the mixing length reduce to the length scale of the largest eddies, as has been assumed many times before (Bretherton 1966; Green 1970; Stone 1972; Held and Larichev 1996). This agreement between the observed eddy size and mixing length means we can accurately estimate the unsuppressed eddy diffusivities using . This approximation is shown in Fig. 6 (red dashed line) and agrees extremely well with the value of derived directly from the tracer fluxes.

We have employed a mixing efficiency of Γ = 0.35, but what sets this value? One attempt at explaining this parameter was the work of Taylor (1915), which is also the first publication to use a mixing length approach (to explain vertical eddy transport of heat in the atmosphere) known to the authors. Taylor (1915) wrote the vertical eddy diffusivity as , where was the typical vertical velocity within an eddy, and *d* was the average height through which an eddy moves in the vertical before mixing with its surroundings. He explains the factor of 0.5 by saying “the air at any given point is equally likely to be in any portion of the path of the eddy so that the average value of *z* − *z*_{0} [the Lagrangian displacement] should be approximately equal to *d*/2.” Our *L* is analogous to Taylor’s *d*. If we assume that our *u*_{rms} represents the rotational velocity of a circular eddy, this would be related to by a factor of *π*/2. Correcting our mixing efficiency by this factor, we obtain , which is very close to Taylor’s factor of , especially if one considers the uncertainties in defining the observed eddy size *L*_{obs}. This provides some physical rationale for the use of Γ = 0.35, beyond the simple fact that it produces the best agreement with *K*^{0}.

## 7. The suppression factor

So far we have focused on *K*^{0}, the unsuppressed diffusivity, using it to understand the relationship between the eddy kinetic energy, size, and mixing in the absence of any suppression effects. But *K*^{0} is a hypothetical construction, not the diffusivity experienced in the real ocean. We will now build on these results to understand eddy diffusivities in the presence of a realistic mean flow. To do so, we empirically calculate the meridional diffusivity using the passive tracer with velocity field **u** = **u**′ + *U*(*y*), where *U*(*y*) is the zonally averaged surface zonal flow [here taken from the ocean state estimate of Forget (2010)]. This diffusivity, which we call *K*_{obs}, is shown as a black solid line in Fig. 6; it is clearly less than in most of the ocean. Now, we test to what extent this reduction in diffusivity can be explained quantitatively by the mixing suppression model discussed in section 2.

As can be seen from Eqs. (2) and (3), and using a constant mixing efficiency Γ, the suppressed eddy diffusivity can be written as

The variables that determine the mean flow suppression are the mean flow speed relative to the eddy phase speed (*c*_{w} − *U*), the eddy wavenumber *k*, and the eddy decorrelation time scale *γ*^{−1}. The wavenumber is simply *k* = 2*π*/*L*_{obs}.

The crucial parameter in determining the strength of the suppression effect is the inverse time scale *γ*. In the stochastic models of Ferrari and Nikurashin (2010) and Klocker et al. (2012b), *γ* arises as the linear damping of potential vorticity anomalies, meant to represent the damping aspect of nonlinear wave–wave interactions. In these models, it is equivalent to the eddy velocity decorrelation time scale, meaning that we can estimate it using the Lagrangian formula

where *K*^{0} is the unsuppressed diffusivity discussed previously. (The subscript _{L} indicates that this is the Lagrangian time scale.) Together with the values of *L*_{obs} and *c*_{w}, which were also discussed in preceding sections, this allows us to test the agreement of Eq. (7) against *K*_{obs}. The comparison is shown in Fig. 6, where the red solid line shows *K*_{mix} as estimated using Eq. (7).

From Fig. 6, it is clear that Eq. (7) captures the suppression effect qualitatively. The strong meridional variation in *K*_{obs} is captured particularly well, with the peaks and valleys located at the right latitudes. However, it is also clear that significant quantitative disagreement in magnitude exists in certain regions. In general, Eq. (7) underestimates *K*_{obs} at most latitudes using *γ*_{L}. There are two possibilities to explain this remaining misfit within the context of Eq. (7): either the mixing efficiency Γ now changes its value with latitude or else our choice of *γ*_{L} was incorrect. It is not possible to distinguish these possibilities through our experiments. Because we have no physical model with which to explain a nonconstant Γ, we prefer to treat it as constant and instead focus on *γ*.

Rather than using *γ*_{L} in Eq. (7), we can find the value that gives the minimum misfit between *K*_{obs} and Eq. (7), keeping the other parameters fixed. This leads to an alternate estimate of *γ*, which we call *γ*_{fit}. In Fig. 7, we plot these two estimates of *γ* side by side. (Note that the figure plots *γ*^{−1}, in order to give units of time.) While *γ*_{L} varies between (4 days)^{−1} and (8 days)^{−1}, *γ*_{fit} is much more constant at approximately (4 days)^{−1}. From this point of view, the fact that *γ*_{L} is in general smaller than *γ*_{fit} is responsible for the overestimated suppression effect in Fig. 6.

The mismatch between *γ*_{L} and *γ*_{fit} suggests a potential shortcoming in the stochastic models of Ferrari and Nikurashin (2010) and Klocker et al. (2012b); these models predict the two quantities should be the same. It is noteworthy that the region of best agreement is the ACC region, where the aforementioned studies were focused. This is one possible reason why the issue was not noticed in those earlier studies. It is well known that the nature of the relationship between Eulerian and Lagrangian scales varies considerably throughout the global ocean (Lumpkin et al. 2002). By examining a wide range of latitudes, our results suggest that some aspects of the model could be refined. Keeping in mind that *γ* was introduced through a linearization of turbulent nonlinear wave–wave interactions it seems unsurprising that this parameter is a source of uncertainty.

Despite the overestimation of the mixing suppression, we emphasize that, given the wide range of latitudes and *K*_{obs} values over this sector, it is evident that Eq. (7) is skillful at reproducing the observed diffusivities and can therefore provide a robust basis for eddy parameterizations. Even employing *γ*_{L} in Eq. (7), the agreement with *K*_{obs} is decent. Even better agreement could be achieved by simply using a constant *γ* = (4 days)^{−1}, although this option is somewhat less satisfying without a deeper explanation for how its value arises.

## 8. Global extrapolation

Motivated by the skill of Eq. (7), we now attempt to extrapolate to the full global ocean using observed eddy kinetic energy and size. We begin with an isotropic unsuppressed diffusivity *K*^{0}, given by the product of *u*_{rms} and *L*_{obs} globally. This is shown in Fig. 8a. We then calculate different suppression factors based on the mean velocities and long-wave linear Rossby wave speed in each direction:

where *c*_{x} is now the zonal phase speed from Eq. (6), *U* is the zonal-mean flow as before, *c*_{y} is the meridional phase speed that we set equal to the depth-averaged meridional-mean flow,^{2} and *V* is the mean meridional velocity. The suppressed eddy diffusivity shown in Fig. 8b is then the minimum between the zonal and meridional diffusivities, *K*^{min} = min(*K*_{x}, *K*_{y}), because we are interested in the cross-stream eddy diffusivity, that is, the direction along which diffusivities are smallest. The resulting eddy diffusivities, shown in Fig. 8b, reveal an extremely high degree of spatial variability over more than an order of magnitude. Low diffusivities are found in regions of strong mean flow, such as the ACC and western boundary currents, with high diffusivities on their flanks. Both the magnitude and spatial variability of eddy diffusivities are consistent with recent studies based on observations (Marshall et al. 2006; Rypina et al. 2012; Naveira Garabato et al. 2011; Abernathey and Marshall 2013).

One aspect of these global diffusivity estimates merits closer consideration. Recent work on eddy mixing in the ACC has highlighted the importance of “hotspots” for cross-frontal exchange in the lee of topography (Naveira Garabato et al. 2011; Thompson and Sallée 2012). Naveira Garabato et al. (2011) hypothesized that these hotspots are due to the mixing length arguments in Eq. (3) breaking down in these regions because these topographic features lead to a nonzonal shear flow, hence violating the assumptions behind the mixing length model in Eq. (3). To understand these hotspots in the ACC in more detail we show *K*^{min}, EKE, and the inverse of the suppression factor in Fig. 9. As shown before (e.g., Naveira Garabato et al. 2011; Thompson and Sallée 2012), the ACC displays strong zonal inhomogeneities, with regions of particularly strong EKE found downstream of major topographic features such as the Kerguelen Plateau, Campbell Plateau, and Drake Passage. Together with the distribution of the mean flow, the mixing length arguments predict strong suppression of eddy mixing over topography with very weak suppression in the lee of these features, leading to eddy diffusivities that are enhanced by more than an order of magnitude in these mixing hotspots. Hence, even though we cannot prove here that the assumptions behind the mixing length model in Eq. (3) are not violated in these regions, we can show that the mixing length approach is capable of predicting very high eddy diffusivities in the lee of topography, which are in the same locations and of the same magnitude as described elsewhere (Sallée et al. 2011; Thompson and Sallée 2012; Abernathey and Marshall 2013).

## 9. Conclusions

Our goal in this study has been to attempt to explain observed patterns of global surface eddy diffusivity using the mixing length approach. In this case, our “observed” diffusivities were generated from kinematic tracer advection experiments using observed surface geostrophic velocities. Provided one properly accounts for the suppression effect of wave propagation relative to the mean flow, we have demonstrated that knowledge of eddy size, EKE, depth-averaged mean flow, and the Rossby deformation radius (needed to calculate the long-wave Rossby wave phase speed) permits the accurate estimation of spatially varying eddy diffusivities using a simple formula. The resulting eddy diffusivities vary by more than an order of magnitude. This exercise can been seen as an updated version of the studies by Holloway (1986) and Stammer (1998), who also attempted to make estimates of eddy diffusivity from satellite observations. These earlier studies lacked two key elements that we have included here: an updated mixing length theory accounting for the suppression effect and a kinematic tracer simulation to provide a ground truth against which to test the theory’s skill.

In the course of making our estimate and validating the mixing length formulas, we gained two additional insights along the way: 1) nonlinear eddies responsible for mixing travel at the Doppler-shifted long-wave Rossby wave phase speed, and 2) there is a direct relationship between observed eddy size and mixing length, but only for the hypothetical unsuppressed diffusivity. As Eq. (3) makes clear, eddy size is crucial for setting the actual mixing length, but the suppression factor also plays an important role globally (and not just in the Southern Ocean). Another conclusion we reached [confirming work by Theiss (2004) and Tulloch et al. (2009)] is that there exists a critical latitude of approximately ±18° that divides an equatorial regime, dominated by wavelike behavior, from the rest of the ocean, where nonlinear eddies prevail. The mixing length theory we employed applies well to mixing in the nonlinear eddy regime, but not the low-latitude regime. Future work will have to develop a new approach to understand what determines eddy diffusivity in this low-latitude region.

It is important to understand what our study does not do: it does not make a complete closure for eddy diffusivity in terms of the mean state, in the vein of Green (1970), Stone (1972), or Held and Larichev (1996). To make such a closure using our approach, one must go farther and predict the EKE and eddy size based on the mean state. Stammer (1998) illustrates of how this might be done using linear baroclinic instability theory. One limitation of linear theory is its inability to account for the inverse cascade, which causes eddies to grow larger than the deformation radius (Smith 2007). As Fig. 5a illustrates, the strength of the inverse cascade depends on the nonlinearity parameter *r*. Another limitation of linear stability analysis that treats each profile independently (as in Stammer 1998; Smith 2007; Smith and Marshall 2009; Vollmer and Eden 2013) is that it cannot allow for nonlocal generation and dissipation of EKE. Grooms et al. (2013) recently showed that the eddy energy cycle can in fact be extremely nonlocal, complicating such closure attempts.

Two final ingredients necessary to turn these ideas into a full-fledged closure scheme in a coarse-resolution ocean model are 1) a theory for the vertical structure of the diffusivity and 2) a treatment of the advective part of eddy transport [related to the Gent and McWilliams (1990) parameterization]. Ferrari and Nikurashin (2010) and Klocker et al. (2012a) showed that it is relatively straightforward to extend the mixing length approach in the vertical, provided that the vertical structure of EKE is known. In the interior of the ocean, the mixing length approach would provide an isopycnal diffusivity. Complicating matters, Smith and Marshall (2009) and Abernathey and Marshall (2013) have demonstrated that isopycnal diffusivities have a nontrivial relationship with the Gent–McWilliams coefficient, especially where there is strong vertical structure in these coefficients. A related study (Bates et al. 2013, manuscript submitted to *J. Phys. Oceanogr.*) addresses this issue through bypassing the Gent–McWilliams approach completely and calculating the eddy-induced velocities in terms of potential vorticity mixing (Marshall 1981).

It is our hope that the ideas presented in this paper can contribute to parameterization efforts by providing simple relations between diffusivity and eddy properties. Beyond this practical goal, we consider these results an incremental step toward better understanding tracer transport by turbulent flows, a problem of fundamental theoretical importance in the ocean and atmosphere. While our results, together with recent studies by Ferrari and Nikurashin (2010) and Klocker et al. (2012a,b), suggest significant progress for extratropical oceanic flows with high zonal symmetry, much remains to be done. The mixing length ideas presented herein clearly fail in the equatorial regime, characterized by values of *r* < 1. Further research is required to understand how to better describe this regime. The issue is not simply of theoretical interest; the meridional ocean eddy heat transport is quite significant at low latitudes (Jayne and Marotzke 2002). Furthermore, although we have gone to great lengths to empirically validate the mixing length formulas in the Pacific sector, we have not done such a detailed comparison for the full ocean surface flow. Instead, we have simply extrapolated the results from the Pacific sector to a global scale. Although the resulting map is generally consistent with observational estimates (e.g., Rypina et al. 2012; Abernathey and Marshall 2013), understanding the detailed structure of the diffusivity tensor in regions of strong anisotropy, such as western boundary currents, should be a high priority for future research. Additionally, more research needs to be done to understand the variations of mixing efficiency, a topic that has seen little attention so far. Finally, because we have focused on cross-stream diffusivity, we have not attempted to include the effects of shear dispersion in our estimates; this process is expected to enhance diffusivity in the along-jet direction in the presence of mean flow shear (Young and Jones 1991).

## Acknowledgments

The authors thank M. Bates, S. Keating, D. Munday, S. Groeskamp, R. Ferrari, M. Nikurashin, M. Spydell, and two anonymous reviewers for very useful comments on the manuscript.

### APPENDIX A

#### Lagrangian Perspective

We framed our discussion of diffusivity and mixing length in section 3 in terms of Eulerian arguments. Here, we briefly review the Lagrangian interpretation. A Lagrangian approach to eddy diffusion was introduced by Taylor (1921). Assuming a homogenous and isotropic eddy field, he showed that one can write *K* as being equal to the EKE times an integral time scale *T*_{L}, with the integral time scale being equal to the time integral of the Lagrangian velocity autocorrelation *R*_{L}:

This approach has subsequently been refined to allow for spatially varying diffusivities and mean flow and to be applicable to ocean float data (Davis 1987, 1991). Numerous applications followed, with one recent example being the Diapycnal and Isopycnal Mixing Experiment in the Southern Ocean (DIMES) (Gille et al. 2012), which used floats to estimate eddy diffusivities in the ACC. Klocker et al. (2012a) demonstrated how wave propagation relative to the mean flow can also reduce the Lagrangian time scale, analogously to the suppression of mixing lengths. Specifically, for a weakly nonlinear eddy field, they showed that the integral time scale can be written as

where is the decorrelation time scale. [This expression rests on the same assumptions as the expression for the mixing length in Eq. (3).] Therefore, *T*_{L} is composed of an exponential decay term and an oscillatory term. The exponential decay term depends only on eddy properties (giving the unsuppressed integral time scale) and the oscillatory term is due to across-current suppression by the mean flow, leading to a shorter integral time scale and therefore a smaller eddy diffusivity. Both the mixing length arguments leading to *K*_{mix} and the Lagrangian approach leading to *K*_{L} give identical results when using the relation *T*_{L} = *γ*^{−1} = *K*^{0}/EKE because the effect on the mean flow on shortening the mixing length *L* is equivalent to the effect on the shortening of the integral time scale *T*_{L}.

### APPENDIX B

#### Tracer-Based Mixing Length

The theoretical expression for the mixing length *L*_{mix} [Eq. (3)] depends only on eddy properties and the mean flow. An alternative, empirical approach to define mixing length, now dependent on a tracer, is (Armi and Stommel 1983; Thompson and Young 2007; Naveira Garabato et al. 2011)

where is the rms tracer fluctuation and is the large-scale tracer gradient. The simple idea behind this definition is that tracer anomalies arise from turbulent mixing across a mean gradient, which means that the amplitude of the tracer anomalies must be related to the distance over which the eddies mix tracer. An eddy diffusivity is then defined as

Note the different mixing efficiency *α* in contrast to Γ in Eq. (3). This mixing length definition also rests on two major assumptions (Naveira Garabato et al. 2011): (i) tracer fluctuations are generated by local stirring of the large-scale tracer gradient, with the advection of tracer variance from upstream regions being negligible; and (ii) the gradient of the tracer varies slowly over the distance *L*_{tracer}.

We calculated *K*_{tracer} from our experiment with (*c*_{w} − *U*) = 0, that is, the velocity field that produces unsuppressed diffusivities, to examine the relationship between and . As shown as dashed lines in Fig. B1, these diffusivities agree very well with a mixing efficiency of *α* = 0.15, which is very close to the mixing efficiency of *α* = 0.16 estimated by Wunsch (1999). This suggests that, in the unsuppressed mixing regime, eddy diffusivities can in principle be estimated equivalently using either the mixing length formation based on the Eulerian velocity field, based on tracers, or from the Lagrangian time scale.

We now repeat the exercise for suppressed diffusivities using realistic *U*(*y*), with both *K*_{obs} and *K*_{tracer} shown as solid lines in Fig. B1. From this figure it can be seen that the eddy diffusivity estimated from Eq. (B2) is approximately a factor of 2 larger than *K*_{obs}. As in Eq. (7), the mixing efficiency itself is quite dependent on the properties of the flow, a conclusion also reached by Thompson and Young (2007) in idealized simulation of baroclinic turbulence. In contrast to the mixing efficiency Γ, we cannot see an obvious relationship of *α* with the mean flow.

## REFERENCES

*Geophys. Res. Lett.,*

**34,**L15606,

*J. Phys. Oceanogr.,*

**20,**758–768.

*Ocean Circulation and Climate: A 21st Century Perspective,*G. Siedler et al., Eds., Elsevier, 185–209.

*CLIVAR Exchanges,*No. 17, International CLIVAR Project Office, Southampton, United Kingdom, 46–48.

**70,**569–602.

*J. Geophys. Res.,*

**115,**C12017, doi:10.1029/2010JC006338.

*J. Geophys. Res.,*

**116,**C09019,

*Remote Sens. Environ.,*

**129,**

*J. Geophys. Res.,*

**114,**C02005,

**72,**198–209.

*J. Geophys. Res.,*

**106**(C7), 13 817–13 836.

*J. Geophys. Res.,*

**109,**C05015, doi:10.1029/2003JC002241.

## Footnotes

^{1}

By absolute eddy phase speed *c*_{w}, we mean the phase speed that would be measured by a stationary observer. For baroclinic Rossby waves, this phase speed is determined by a dispersion relation that includes a Doppler shift by the depth-averaged zonal-mean flow , which is not necessarily the same as *U* at the surface.

^{2}

Whereas the depth-averaged zonal-mean flow very accurately explains the Doppler shift of mesoscale eddies, this is not the case for the meridional phase speed of eddies; see Klocker and Marshall (2013, manuscript submitted to *Geophys. Res. Lett.*) for details. Nevertheless, describing the meridional phase speed with the depth-averaged meridional-mean flow has negligible consequences for the estimates of eddy mixing.