## Abstract

Temperature–salinity profiles from the region studied in the North Atlantic Tracer Release Experiment (NATRE) show large isopycnal excursions at depths just below the thermocline. It is proposed here that these thermohaline filaments result from the mesoscale stirring of large-scale temperature and salinity gradients by geostrophic turbulence, resulting in a direct cascade of thermohaline variance to small scales. This hypothesis is investigated as follows: Measurements from NATRE are used to generate mean temperature, salinity, and shear profiles. The mean stratification and shear are used as the background state in a high-resolution horizontally homogeneous quasigeostrophic model. The mean state is baroclinically unstable, and the model produces a vigorous eddy field. Temperature and salinity are stirred laterally in each density layer by the geostrophic velocity and vertical advection is by the ageostrophic velocity. The simulated temperature–salinity diagram exhibits fluctuations at depths just below the thermocline of similar magnitude to those found in the NATRE data. It is shown that vertical diffusion is sufficient to absorb the laterally driven cascade of tracer variance through an amplification of filamentary slopes by small-scale shear. These results suggest that there is a strong coupling between vertical mixing and horizontal stirring in the ocean at scales below the deformation radius.

## 1. Introduction

The temperature–salinity (*T*–*S*) relationship of the global ocean is often characterized by thermohaline fluctuations with little signature in density. Vertical profiles display *T*–*S* intrusions of tens to hundreds of meters that tend to cancel their joint effect on density (e.g., Schmitt 1994; Lilly et al. 1999; Ferrari and Polzin 2005). Horizontal maps of *T*–*S* properties show abundant examples of thermohaline fronts of a few kilometers with strong compensation between the *T* and *S* gradients (e.g., Roden 1977; Rudnick and Ferrari 1999; Ferrari and Polzin 2005). The processes that generate this *T*–*S* variability have concerned oceanographers for decades, because these processes provide clues to the dynamics that set the *T*–*S* relationship of the global ocean.

At least two explanations have been proposed for the generation of thermohaline structures with weak density signature in the ocean interior below the surface mixed layer. The first one invokes double-diffusive instabilities at molecular scales as the source of the observed variability. The idea is that *T*–*S* differences in water masses with stable density structure can be unstable because of the different molecular diffusivities of heat and salt. Vertical *T*–*S* differences generate double-diffusive layers—that is, well-mixed laminated structures up to 100-m thick and of large lateral extent (Schmitt 1994). Lateral *T*–*S* gradients drive double-diffusive intrusions because of the lateral buoyancy flux convergences due to differential heat and salt fluxes (Ruddick and Gargett 2003). In both cases the resulting filaments are characterized by strong *T*–*S* compensation, because double diffusion is arrested as soon as *T* and *S* are not opposing in their effect on density. Even though double diffusion is clearly responsible for the observed *T*–*S* variability in some regions of the ocean, such as the Sargasso Sea (Ruddick and Gargett 2003), it is not clear how ubiquitous this process really is.

A second explanation of the observed *T*–*S* compensation is that the large scale *T*–*S* water mass anomalies are broken up by the mesoscale eddies, yielding strong small-scale thermohaline fronts. Mesoscale motions in the ocean interior are mostly along-isopycnal surfaces and are therefore not effective at increasing horizontal density gradients but are efficient at stirring properties whose isopleths are inclined to isopycnals (MacVean and Woods 1980).

The consequence is that mesoscale eddies can generate compensated *T*–*S* gradients but not density gradients. Ferrari and Polzin (2005) show convincing evidence that this second mechanism is at work in the generation of *T*–*S* variability in the eastern North Atlantic at the water mass contrast between the salty Mediterranean Outflow Water and the fresher North Atlantic Central Waters. Using numerical simulations of a fully turbulent, time-evolving quasigeostrophic (QG) flow, Klein et al. (1998) show that this mechanism occurs generically in regions where baroclinic eddies coincide with mean lateral gradients of temperature or salinity.

The two mechanisms imply different pathways of *T*–*S* variance through the ocean. In the first case, *T*–*S* gradients are generated as a result of instabilities at molecular scales, and the transfer of variance is from small to large scales. In the second case, mesoscale eddies generate small-scale gradients by twisting and folding the large scale *T*–*S* anomalies. There are numerous observations of vertical *T*–*S* gradients at scales below a few hundred meters and of thermal dissipation rates at molecular scales (e.g., Polzin 1996). This has allowed a careful examination of double-diffusive processes in the ocean. Comparatively less is known about the generation of small-scale gradients by eddy stirring, because observations on lateral scales of 1–100 km are very rare. As a result, there has been a tendency to focus exclusively on double-diffusive processes in interpreting the *T*–*S* variability in the ocean. The goal of this paper is to quantify the efficiency of eddy stirring at creating small-scale gradients and to test the hypothesis that the transfer of variance from large to small scales is the primary mechanism for the generation of compensated *T*–*S* variability in the ocean interior.

We investigate the stirring of thermohaline anomalies in a fully turbulent, stratified flow using the large-scale *T*–*S* patterns found in the Eastern North Atlantic at the water mass contrast between Mediterranean and North Atlantic Central Waters as a background state. The key assumption is that the interior ocean dynamics on scales of 1–500 km can be described by the quasigeostrophic approximation. This allows us to explain the structure of the thermohaline anomalies in terms of well-known properties of forced QG turbulence, as suggested by Klein et al. (1998). The simulated flow stirs the observed large-scale *T* and *S* gradients, producing a forward cascade of *T*–*S* variance. The simulated *T*–*S* structure is quantitatively consistent with the North Atlantic Tracer Release Experiment (NATRE) observations, producing 0.2-psu salinity anomalies at the depth of the Mediterranean salinity tongue. Furthermore, the simulated generation of temperature variance by geostrophic stirring at the depth of the Mediterranean salinity tongue is near to but somewhat *larger* than microstructure measurements of thermal variance dissipation, as found by Ferrari and Polzin (2005). Given that lateral mesoscale generation is typically neglected, the fact that a moderately energetic eddy field can generate sufficient variance to account for the observed dissipation is significant. Moreover, the approximate agreement between the simulated generation and measured dissipation rates is quite remarkable because it links processes spanning five orders of magnitude in spatial scales.

If lateral stirring by the mesoscale field generates a forward cascade of *T*–*S* variance at the submesoscales, what halts this cascade? One possibility is that the formation of isopycnal *T*–*S* gradients is arrested in the horizontal by three-dimensional turbulence as a result of internal wave breaking. However, this scenario would lead to quite small *T*–*S* features. Ledwell et al. (1998) find that filaments of a tracer released in the main thermocline, as part of NATRE, had widths of about 2 km, once the initial tracer patch had been stirred by the mesoscale field for about six months. Given the observed stretching rate of about 3 × 10^{−7} s^{−1} and the tracer-based estimate of an isotropic turbulent diffusivity *κ* ≃ 1 × 10^{−5} m^{2} s^{−1}, a steady-state balance between stretching and turbulent mixing would give filaments of approximately 6 m. Some enhancement of the tracer filament width might be expected by internal wave shear dispersion (Young et al. 1982), which couples turbulent diffusion and internal wave shear to produce an enhanced horizontal diffusivity. Ledwell et al. (1998) estimate that shear dispersion would give filament widths of 300 m, still an order of magnitude below the observed value.

Clues into what arrests the cascade of the *T*–*S* cascade come from a careful examination of the generation of the *T*–*S* filaments, which in many ways resemble the filaments observed in stratospheric tracers. Stratospheric ozone, for example, is filamentary in the horizontal, with typical widths of order hundreds of kilometers, and laminated in the vertical, with typical thicknesses of order 1 km (Orsolini 1995; Appenzeller and Holton 1997). The well-observed and straight-forwardly modeled nature of the three-dimensional structure of stratospheric tracers has enabled progress on the problem of what sets their horizontal and vertical scales. Haynes and Anglade (1997) argued that the observed stratospheric tracer structure can be understood as resulting from a coupling between horizontal strain, vertical shear, and turbulent diffusion: tracer filaments are continuously strained in the horizontal and sheared in the vertical, so that the horizontal and vertical scales of the filaments shrink in unison. By considering the evolution of a tracer blob in the presence of a constant horizontal strain, vertical shear, and laplacian diffusion, they show that the ratio of horizontal–to–vertical scales of the tracer filaments (defined as *α*) is proportional to the ratio of shear over strain. Thus, when the shear is much larger than the strain, the vertical scale of the tracer will be much smaller than its horizontal scale. Diffusion will therefore first act to limit the minimum vertical scale, but, by its linkage to the horizontal strain, it will also effectively limit the minimum horizontal scale, yielding an effective horizontal diffusion that is much larger than the bare background isotropic diffusivity. This model appears to match well with observations and simulations of stratospheric tracer distributions (Legras et al. 2003; Haynes and Vanneste 2004; Legras et al. 2005).

Here we show that this relationship holds in a fully turbulent, geostrophic flow and, further, that *α* is then proportional to *N*/*f*, where *N* is the buoyancy frequency and *f* is the Coriolis parameter. Geostrophic stirring drives a cascade of tracer variance both in the horizontal and in the vertical, because the vertical shear continuously tilts the filaments generated by lateral stirring. The tracer filaments are much thinner in the vertical than in the horizontal, and the cascade of tracer variance is arrested once the thickness, not the width, of the filaments becomes of the same order as the scale of the isotropic turbulence that characterizes internal wave breaking. The picture is somewhat different than the atmospheric case, for here the strong vertical shears arise from the forward cascade of potential enstrophy and are both larger in magnitude and more vertically structured than the background mean shear. This conjecture is supported by simulations in which vertical diffusion alone acts to absorb the tracer variance cascade. The coupling of horizontal and vertical gradients in the forward tracer variance cascade provides a route to dissipation of the laterally generated temperature variance. We further show that the tracer filament widths span a wide range, because the distribution of tracer slopes is quite broad around the mean value of *N*/*f*. As a result, the range of filament widths is quite broad, even if the vertical mixing scale is sharp. For realistic values of turbulent diffusivity, the widths of tracer filaments span the range of ∼100 m to ∼10 km, which includes Ledwell et al. (1998)’s observations.

The paper is organized as follows: Observational data provided by NATRE is reviewed in section 2, setting the bar for our numerical investigation. In section 3 we rationalize the model framework, discuss the method by which salt and temperature with arbitrary mean gradients can be advected by a quasigeostrophic flow, and present a theory for the production of compensated *T*–*S* fine-structure by mesoscale stirring. In particular, we consider how vertical *T*–*S* gradients are generated by horizontal straining and shear in the velocity field. The setup for and main results of our central simulation are presented in section 4. Details of the generation of vertical gradients by horizontal straining are discussed in section 5. In section 6 we demonstrate that the coupling of horizontal and vertical gradients in tracer filaments allows background vertical diffusion to effectively limit the forward cascade of tracer variance; thus, the horizontal filament widths observed in NATRE can be understood as being set by vertical mixing alone. The paper concludes with section 7.

## 2. A review of observations

The eastern subtropical North Atlantic was the subject of an intensive series of field programs in 1991–93 as part of the North Atlantic Tracer Release Experiment (Ledwell et al. 1998) and the Subduction Experiment (Joyce et al. 1998). These observations provide a detailed picture of the mean flow patterns and the thermohaline variability in the region and are used as ground truth for the study of eddy stirring of thermohaline anomalies in the ocean.

The NATRE campaign was conducted southwest of the Canary Islands (26°N, 29°W) during April 1992. The first part of the campaign was devoted to collecting vertical profiles with a high resolution profiler (HRP) that measured *T*, *S*, and velocity shear every meter. Most profiles reached 2000 m, whereas a few reached 4000 m. One hundred profiles were regularly spaced on a 400 km × 400 km grid, and 44 profiles were clustered near the center of the grid. Ferrari and Polzin (2005) provide a thorough review of the observational program. For present purposes, we focus on the *T*–*S* relationship from all the profiles (Fig. 1). The *T–S* curves are very tight in the thermocline waters, characterized by weak isopycnal gradients. Immediately beneath the thermocline, contrasts between saline Mediterranean Water (MW) and fresher Antarctic Intermediate Water and Labrador Sea Water result in strong *T*–*S* gradients along-isopycnal surfaces. The *T*–*S* variability disappears at even greater depths.

The HRP survey in NATRE was near the center of a moored array deployed as part of the Subduction Experiment (Weller et al. 2004). Mooring C at 25.5°N and 29°W was equipped with vector averaging current meters at depths of 200, 300, 310, 1500, and 3500 m. Here 2-yr time series (summer 1991–spring 1993) of velocity and temperature are available at most depths. These data are used to test the skill of our numerical simulations in reproducing the observed geostrophic eddy field.

At the end of the NATRE campaign, an anthropogenic tracer (SF_{6}) was released on isopycnal *σ _{θ}* = 26.75 kg m

^{−3}, at a depth of about 313 m (Ledwell et al. 1998). The tracer was injected over the period of 5–13 May 1992 and sampled immediately thereafter over the period of 14–31 May, an average time difference of two weeks (1.2 × 10

^{6}s). The tracer was injected as a set of nine streaks of length

*L*= 10 km. These nine streaks defined an initial patch of about 45 km × 45 km extent. At two weeks’ time, the tracer streaks had formed a nearly continuous patch that had been stretched zonally and compressed meridionally by a factor of 2 or so.

_{i}Subsequent sampling of the tracer in October and November 1992 showed that the tracer was teased into filaments by the mesoscale eddy field (Fig. 2). There was a clear distinction between the tracer containing filaments and the surrounding water. Ledwell et al. (1998) inferred a total filament length of about 1800 km and, from the October tows, an rms streak width of 3 km (and a range of 2–5 km). The along-streak growth was assumed to have resulted from a rate of stretching (*σ*) along each streak element. This stretching rate was estimated assuming an exponential rate of increase (Garrett 1983) for the streak length *L*: *σ*_{6 months} = log(*L*/*L _{i}*)/6 months = [3 ± 0.5] × 10

^{−7}s

^{−1}. Assuming a balance between stretching and diffusion in the across-streak direction, Ledwell et al. (1998) estimated a horizontal diffusivity of

*κ*= 1–3 m

_{h}^{2}s

^{−1}.

How are the distribution of the tracer and the *T*–*S* measurements related? A corollary to the hypothesis we explore in this paper is that the distribution of the injected tracer should, at some stage, mimic the distribution of *T* and *S*. Here the *T* and *S* variance are continually forced by their large-scale along-isopycnal gradients, set by water mass contrasts between the salty Mediterranean Waters and the fresher North Atlantic Central Waters. The injected tracer forms a freely dispersing cloud, as described. When large-scale straining (at the radius of the filament cloud) has teased the blob into a series of filaments, the filaments presumably achieve widths that represent a balance between large-scale straining and submesoscale mixing (at a few kilometers and below). At this stage, the processes controlling the distribution of the injected tracer and the *T*–*S* fronts should be equivalent.

## 3. The model framework

Large- and mesoscale motions in the ocean interior are directed primarily along-isopycnal surfaces and are therefore not effective at increasing horizontal density gradients, but they are efficient at stirring tracers whose isopleths are inclined to mean isopycnals (MacVean and Woods 1980). “Spice,” the combination of *T* and *S* that does not contribute to density, is effectively a passive tracer that can be stirred by along-isopycnal velocities (this will be made mathematically explicit below). The consequence is that mesoscale eddies can generate sharp spice gradients—that is, *T* and *S* gradients that are compensated in their effect on density.

We seek here to model the isopycnal stirring of *T* and *S* at sufficient resolution and detail to draw quantitative conclusions about the validity of the mesoscale stirring hypothesis. The problem requires a realistic eddy field that spans scales from above the deformation scale to about 1 km. A standard numerical approach to this problem would employ a primitive equation or nonhydrostatic model in a regional domain configuration, run at relatively high resolution. However, even at state-of-the-art resolutions, such models are not numerically converged (Capet et al. 2008), because unbalanced instabilities (physical and artifactual) are excited at small scales and such computations are still limited by computational power. Moreover, it is difficult in such models to disentangle numerical from physical diffusion, isopycnal from diapycnal mixing, and balanced from unbalanced motions. Here we avoid these complications and stir the passive components of temperature and salinity with a QG model in which eddies are generated by the baroclinic instability of a prescribed background shear and buoyancy gradient. Because the deformation-scale dynamics that generate the eddies precisely characterize the asymptotic regime of validity for quasigeostrophy (QG) this is a very accurate approximation at the mesoscale. At smaller scales we ask the question, how much of the observed finescale structure can be explained by balanced mesoscale stirring?

### a. Advecting temperature and salt with quasigeostrophic flow

In QG models, the stratification is fixed in time, but it is nevertheless possible to advect spice, the passive component of the *T*–*S* field, as shown originally by Klein et al. (1998). Here we briefly outline a generalization of their method and relegate details to appendix A. Assuming a linear equation of state, the nondimensional potential density *ρ* and spice *γ* are

The dimensional potential density and spice are *ρ*_{0}(1 + *ρ*) and *ρ*_{0}(1 + *γ*), and the constants *ρ*_{0}, *T*_{0}, and *S*_{0} are reference values of density, temperature and salinity, respectively. Subtracting and adding the inviscid conservation equations for salinity and temperature (and neglecting vertical advection of eddy components, consistent with QG scaling), respectively, one can write the conservation equations for density and spice:

where *D*/*Dt* = ∂* _{t}* +

*u*∂

*+*

_{x}*υ*∂

*= ∂*

_{y}*+*

_{t}**u**·

**∇**is the horizontal advective derivative,

*D*and

_{ρ}*D*are dissipation terms, and the asterisk superscripted variables [

_{γ}*ρ** =

*ρ**(

*z*) and

*γ** =

*γ**(

*z*)] are related to the prescribed mean stratification, namely

*N*

^{2}= −

*g dρ**/

*dz*and

*N*

_{γ}^{2}= −

*g dγ**/

*dz*. The operator

**∇**= (∂

*, ∂*

_{x}*) is the horizontal gradient, unless stated otherwise. Eliminating*

_{y}*w*between the two equations leads to a purely horizontal advection equation for the tracer

*C*:

and *D _{C}* is the appropriate linear dissipation operator.

^{1}The equations for

*C*and density in Eqs. (3.2) and (3.3) are used to fully determine the evolution of spice,

*T*, and

*S*. Note that in the QG framework, the tracer

*C*is exactly conserved along isopycnals, unless an explicit diffusion or mixing mechanism is coded into the model. This is another advantage of a QG model versus primitive equation models that inevitably suffer from spurious diapycnal diffusion.

To include the observed mean *horizontal* gradients of temperature and salinity, we construct mean horizontal gradients of buoyancy and spice, and from these the mean horizontal gradients of the conserved tracer **Γ**_{C} = **∇***C*. The conservation equation for the eddy component of *C* (overbars denote horizontal and time average; primed quantities are deviations from the average) is then

where *N* ^{2}, *N _{γ}*

^{2},

**Γ**

*, and*

_{C}**u**= (

*u*, ) are prescribed functions of

*z*. The velocity and density perturbations

**u**′ and

*ρ*′ are obtained by solving the quasigeostrophic potential vorticity equation forced with the observed stratification

*N*

^{2}and mean velocities

**u**as described below. Given the solution to the passive tracer equation in (3.4) and the definition of

*C*in (3.3),

*T*and

*S*can be reconstructed from their relationships to density and spice in (3.1). Details of the method are provided in appendix A.

### b. Potential vorticity dynamics and the three-dimensional potential enstrophy cascade

The velocity field that stirs the tracer *C*′ is generated from the QG potential vorticity (QGPV) conservation equation

where *q*′ = ∇^{2 }*ψ*′ + ∂* _{z}*(

*f*

^{2}/

*N*

^{2}∂

*′) is the eddy QGPV,*

_{z}ψ*ψ*′ is the horizontal eddy streamfunction, and

**Γ**

_{q}=

**∇**

*q*is the mean horizontal QGPV gradient, given by

Note that Eq. (3.5) is isomorphic to the equation for *C*′ in (3.4). Both the tracer and QGPV can be considered as Lagrangian tags of water parcels advected by the flow. The only difference between the two is, of course, that the velocity is dynamically linked to the QGPV, whereas the tracer *C*′ has no feedback on the velocity. Nevertheless, if one knows the structure of QGPV that develops as eddies stir its mean gradient, one can infer the structure that the tracer will develop.

Charney (1971) argued that for small enough scales, one can treat *N* as constant, rescale the vertical coordinate as *z̃* = *Nz*/*f*, and write the QGPV as *q*′ = ∇̃^{2}*ψ*′, where ≡ (∂_{x}, ∂_{y}, ∂_{z̃}). In the vertically stretched coordinates, the QG potential vorticity advection problem is isomorphic to the stirring of vorticity in two dimensions. One can therefore apply the arguments of Kraichnan (1967) for the downscale cascade of enstrophy in two dimensions to the downscale cascade of potential enstrophy, /2, in three dimensions. The fundamental difference is that potential enstrophy is cascaded isotropically in the stretched coordinates; that is, *enstrophy is cascaded at a ratio N*/*f* in the vertical and horizontal directions.

Following Kraichnan, the spectrum 𝒫_{q} of potential enstrophy (QGPV variance) is predicted to be

where *μ*^{2} = *k*^{2} + ℓ^{2} + *λ*^{2} is the magnitude of the three-dimensional wavenumber in the stretched coordinate space, and the spectrum is defined such that

with 〈·〉 denoting a vertical integration over the depth of the fluid. Multiplying (3.5) by *q*′ and integrating over the domain, one finds that *η*, the rate of downscale spectral transfer of potential enstrophy, is −〈 · **Γ**_{q}〉.

The theory can be extended to nonconstant *N* by expanding *q*′ and *ψ*′ onto Fourier modes in the horizontal and eigenfunctions *F _{λ}*(

*z*) in the vertical:

where the eigenfunctions *F _{λ}*(

*z*) are solutions of the Liouville eigenvalue problem

and the eigenvalues *λ* are the internal deformation wavenumbers. If *N* is constant, then *λ _{n}* =

*nπf*/

*NH*(for mode

*n*), which makes clear the stretching of the vertical scale of motion by

*N*/

*f*. If the stratification

*N*is not constant, as in the present problem, then Charney’s argument still applies; however, it requires that the scale of significant variations in

*N*be larger than the vertical scale of the turbulence. In the spectral domain, the QGPV and streamfunction are related as

Here the symmetry of the inversion operator in (*k*, ℓ, *λ*) is obvious; however, when *N* is not constant, the spectral transfers may be nonlocal at large scale (Hua and Haidvogel 1986). The result is that for transfers among small vertical wavenumbers, there is an enhancement of first baroclinic energy and an inhibition of transfers to the barotropic mode (Fu and Flierl 1980; Smith and Vallis 2001); however, transfers involving larger vertical wavenumber remain isotropic. In the present problem, we are concerned with stirring and mixing at horizontal scales an order of magnitude smaller than the deformation scale; therefore, the vertical scales are sufficiently small that one can expect isotropy in (*k*, ℓ, *λ*).

### c. Energy and the three-dimensional cascade of density variance

Charney’s argument for an isotropic cascade of potential enstrophy in the stretched coordinates implies an equipartition of kinetic, 〈 + 〉/2, and available potential energies, *g*^{2}〈/*N* ^{2}〉/2. The total energy spectrum according to (3.7) is *μ*^{2}|*ψ̂*_{kℓλ}|^{2} ∝ *μ*^{−3}, and hence the available potential energy (APE) and density spectra might be expected to be ∝ *μ*^{−3}. The argument, however, fails because, for nonconstant *N*, the bulk of the kinetic and available potential energies are trapped in the first baroclinic mode (Fu and Flierl 1980; Smith and Vallis 2001), whereas the equipartition between kinetic and available potential energies holds at small horizontal and vertical scales. The density variance spectrum (which we define as 𝒫_{ρ}) can be easily computed from the potential energy spectrum, if one accounts for the fact that the bulk of the potential energy is in baroclinic mode 1. The result is

where we used the fact that |*ψ̂*_{kℓλ1}|^{2} ∝ *K*^{−5} in the potential enstrophy cascade regime.

The steep spectral roll-off in 𝒫_{ρ} confirms that, in QG models, eddy stirring generates little variability in *ρ*. Density variance remains trapped at large scales in the first baroclinic mode.

### d. The three-dimensional cascade of compensated temperature and salinity variance

The analogy between the QGPV and tracer equations suggests that Charney’s arguments for the spectrum of potential enstrophy should also apply to the tracer variance 〈/2〉. Defining 𝒫_{C}(*μ*) as the spectrum of the tracer variance, with the same notation as in (3.8), one expects

where *χ* = −〈 · **Γ**_{C}〉 is the rate of generation of tracer variance [obtained by multiplying (3.4) by *C*′ and integrating over the domain]. The forward cascade of tracer variance is therefore expected to be isotropic in (*k*, ℓ, *λ*); therefore, the tracer, like the QGPV, should exhibit horizontal and vertical scales of motion that maintain a ratio of *N*/*f* throughout the inertial range.

The passive tracer spectrum is predicted to have a shallow *μ*^{−1} slope, whereas the density is expected to have all its variance concentrated in baroclinic mode 1 at large scales. Because Γ, *T*, and *S* are all linear combinations of *C* and *ρ*, their high-wavenumber spectrum is also dominated by a *μ*^{−1} spectral roll-off. This in turn implies that *T* and *S* are expected to be rich with submesoscale structure but compensating in their effect on density. Klein et al. (1998) show convincing numerical evidence for this idea in a stratified quasigeostrophic simulation for an idealized flow regime. In the present paper, we apply this model framework to the observed mean structure in the NATRE region. One goal is to assess whether this idea is sufficient to *quantitatively* account for the observed *T*–*S* structure.

## 4. Central simulation of the NATRE region

The central calculation is performed on a domain the size of 1500 km × 1500 km × 5050 m domain with 1.5-km resolution in the horizontal and a variable vertical resolution, ranging from 35 m near the top to 120 m near the bottom for a total of 80 vertical levels. The NATRE HRP data used to initialize the simulations were also interpolated on the same vertical grid. The calculation is horizontally doubly periodic and uses a fully dealiased fast Fourier transform method to compute the nonlinear terms. This representation also allows the use of a very scale-selective exponential cutoff filter to absorb potential enstrophy and tracer variance and is therefore almost completely inviscid down to the smallest resolved horizontal scale. (In some explicit cases discussed later, the tracer variance filter is replaced by either explicit horizontal or vertical diffusion.) The model uses a finite-difference representation in the vertical and is bounded above and below by flat, rigid surfaces. There is no explicit representation of surface buoyancy here; however, 1) these effects are included approximately by the finite-difference representation (see Tulloch and Smith 2009 for details) and 2) the focus in this paper is on the interior, far below the depths potentially affected by surface dynamics. The model is essentially the same as that used by Smith and Vallis (2002). The Coriolis gradient *β* is set to zero, and eddy energy is removed from the flow via quadratic drag at the bottom. These latter two choices were made to 1) allow straightforward control of the eddy energy, 2) produce a horizontally homogeneous field of eddies, and 3) avoid the spurious emergence of jets at low eddy energy levels typical of simulations with flat bottom topography. We emphasize that our focus here is on the forward cascade, and these choices are intentionally made to produce as realistic eddy field as possible while avoiding many of the unresolved issues regarding the equilibration scale and energy of mesoscale eddies.

The high-resolution profiles of *T* and *S* from NATRE are used to generate the mean state for the model as follows: First, the United Nations Educational, Scientific and Cultural Organization (UNESCO) equation of state is used to generate the background density *ρ**(*z*) and the mean linear gradient of density **Γ**_{ρ} (see appendix B). A linear equation of state is defined from the data, with *T*_{0} = 6 C, *S*_{0} = 35.2 psu, *β _{T}* = 1.25 × 10

^{−4}K

^{−1},

*β*= 7.69 × 10

_{S}^{−4}psu

^{−1}, and

*ρ*

_{0}= 1028 kg m

^{−3}—chosen to best match the potential density measured at the depth of the Mediterranean salinity tongue. The coefficients of the linear equation of state for these reference parameters are used to compute the background spice,

*γ**(

*z*), and the mean linear gradient of spice

**Γ**

_{γ}according to Eq. (3.1). The background stratification profiles

*N*

^{2}and

*N*

_{γ}^{2}are computed from vertical derivatives of

*ρ**(

*z*) and

*γ**(

*z*), respectively. The first deformation radius derived from the mean stratification is 41 km. The model background temperature and salinity fields are then reconstructed from density and spice using the linear equation of state (3.1). Thus the background temperature and salinity—

*T**(

*z*) and

*S**(

*z*), respectively—differ slightly from the observed profiles as a result of the requirement that a linear equation of state connect all variables in the model. Finally, the mean velocity field

**u**is computed from

**Γ**

_{ρ}via thermal wind balance in (A.1). The resulting profiles of

*N*

^{2},

*N*

_{γ}^{2}, and

**u**(

*z*) as well as the horizontal gradients of

*T*and

*S*(

**Γ**

*and*

_{T}**Γ**

*, respectively) are shown in Fig. 3.*

_{S}The steady-state and horizontally and time-averaged eddy velocity is plotted as a function of depth in Fig. 4. An observational estimate, computed from the moored current array in NATRE, is also plotted for comparison. The model eddy kinetic energy (KE) is significantly larger than observed near the surface. We speculate that this is due to the lack of any surface energy dissipation mechanism, such as thermal drag, but we do not pursue this issue further in this paper. The model energy is close to the observed value near the depth of the salinity tongue, the focus depth of the present study. Moreover, despite the excess eddy kinetic energy near the surface, the vertical shear of the horizontal eddy velocity is very close to the shear observed in NATRE—both are plotted in the left panel of Fig. 5. (Shear from NATRE data is computed through thermal wind from 100 density profiles on the 400 km × 400 km NATRE square grid.) Note also that the mean shear contributes only about 10% of the total shear. The efficient production of so much vertical structure in the eddy field by large-scale, low-mode baroclinic instability is a major focus of this paper and will be discussed further in the next section.

A similar comparison of the horizontal gradients of the velocity field (the strain and vorticity) would be equally useful to consider, but the data grid is too coarse to compute horizontal derivatives of velocity. Nevertheless, Ledwell et al. (1998) estimate the mesoscale strain as *σ*_{6 months} = [3 ± 0.5] × 10^{−7} s^{−1}, based on the stretching rate of a tracer filament at the level of the tracer release (see section 2). The profile of strain from the model, plotted over the same depth range as the shear data, is shown in the right panel of Fig. 5. The simulated values are up to an order of magnitude larger than the data-based estimate. However, the actual *stretching* rate of a tracer can be significantly less than the strain,^{2} and the energy of the simulated flow is higher than the observed flow at 300-m depth.

Snapshots of the temperature, salinity, and density fields at a depth of *z* = −1100 m (near the level of the MW outflow) in a subregion of the domain are shown in Fig. 6. Temperature and salinity exhibit sharp, finescale fronts that are completely compensating on density below horizontal scales of about 100 km. Power density spectra of spice and density, shown in Fig. 7, exhibit −1 and −5 slopes, respectively, consistent with the physical space picture and theoretical expectations. A spectral roll-off of −1 is also reported from the analysis of along-isopycnal *T*–*S* fluctuations in the NATRE campaign (Ferrari and Polzin 2005).

The large *T*–*S* excursions that are apparent in the NATRE *T*–*S* relationship (Fig. 1) can be explored directly in the simulation. Figure 8b shows a *T*–*S* diagram constructed from subsurface (below 500 m) eddy temperature and salinity in a 400-km subregion (the size of the large NATRE grid). One finds that the width of the cloud is quantitatively similar to that observed over the same depths in NATRE, shown in Fig. 8a.

If lateral straining by mesoscale eddies generates the observed *T*–*S* variance, then some physical mechanism must dissipate the variance in steady state. Locally, the strain field stretches compensated filaments of *T* and *S*, whereas diffusion acts to spread them (the width of the filaments results from a balance between these two actions). In the numerical model, tracer variance is removed by a scale-selective filter; therefore, the width of filaments and level of variance in a vertical profile are set by a balance between mesoscale straining and filter dissipation. The diffusive power of the filter is adjusted to ensure that this dissipation scale is resolved, and as long as this is the case, one can take the variance generation as independent and determined solely by the large-scale flow. In the main simulation, we have used sufficient resolution to allow for 2-km filaments, the width observed in NATRE (Ledwell et al. 1998). The fact that the *T*–*S* cloud in the simulation also matches that observed in NATRE (see Fig. 8 again), with no other mechanism available to generate that variance, is consistent with the mesoscale stirring hypothesis. We address the physical dissipation mechanism in sections 5 and 6.

An advantage of the simulation is that it can be used to assess eddy statistics that are difficult to measure directly from the data. In particular, the mesoscale variance generation

can be computed without the need to introduce a poorly constrained mesoscale eddy diffusivity, as was done in Fig. 7 of Ferrari and Polzin (2005). For the central simulation, the mesoscale variance generation is computed and shown in Fig. 9 (light gray bars). Microstructure measurements of the dissipation of thermal variance at molecular scales (dark gray bars) and variance generation by turbulent motions acting on the mean diapycnal gradient (medium gray bars), both computed from NATRE data, are shown in the same figure. The mesoscale generation is small in the upper thermocline, where turbulent generation is sufficient to account for dissipation. At the Mediterranean salinity tongue level, however, dissipation exceeds turbulent generation, and the mesoscale generation more than makes up for it. Ideally, the sum of mesoscale and turbulent generation would match dissipation at every level, but it is too much to expect exact consistency between the simulated and measured values. We know, for example, that there is too much kinetic energy in the simulation and likely too much strain as well. Moreover, the averaging period and the area in the simulation are very different from that over which the measured values are computed. Finally, of course, this is a quasigeostrophic model with finite resolution. The widths of the bars on the simulated values represent the largest and smallest values of area-averaged generation over the averaging period and therefore do not take any of these factors into account. Nevertheless, it is remarkable that the simulation can come close to matching the eddy generation of thermal variance at hundreds of kilometers with in situ measurements of thermal dissipation at scales of a few centimeters.

## 5. Generation of vertical gradients by horizontal straining

The traditional view of finescale temperature variance, due to Osborn and Cox (1972), assumes that turbulence is produced by internal wave breaking, resulting in a turbulent vertical velocity field that stirs the large-scale mean vertical temperature gradient, thus generating temperature variance at the scale of breaking waves. The temperature variance generated is assumed to be dissipated isotropically at the molecular dissipation scale, yielding a relationship between the eddy temperature diffusivity and the molecular diffusion coefficient and thus providing a direct estimate of the former from microstructure measurements of the latter. Ferrari and Polzin (2005) review this picture in great detail and demonstrate that it neglects the generation of thermal variance by mesoscale geostrophic stirring along-isopycnals. Here we investigate the detailed mechanism by which the mesoscale eddy field generates the gradients that ultimately contribute to dissipation at middepth in the NATRE region.

### a. Tracer aspect ratio in the three-dimensional tracer cascade

As discussed in sections 3b–3d, the forward cascade of tracer variance is three-dimensional and therefore small vertical scales are readily produced. Moreover, because the cascade is expected to be isotropic at submesoscales in the wavenumber space (*k*, ℓ, *λ*), tracer filaments should maintain a scale-independent aspect ratio of *N*/*f* throughout the potential enstrophy inertial range. This is demonstrated in Fig. 10, which shows a horizontal and vertical slice of the salinity field in a subregion of the central simulation. Overlaid on the vertical slice are contours of the total (mean + eddy) density field, and a steeper white line with slope *f*/*N* (with *N* evaluated at the appropriate depth), oriented next to a filamentary structure. The filament is clearly much steeper than the isopycnals, and its slope (inverse aspect ratio) is consistent with the expected *f*/*N* scaling.

More general consistency with the theory is demonstrated as follows: An idealized tracer *ϕ* is driven by a unitary, zonal mean gradient **Γ**_{ϕ} = **x̂** and advected by the velocity field of the central simulation (this is actually the tracer called *A* in appendix A). The vertical structure of the ideal tracer is projected onto the vertical modes defined in (3.9) and onto horizontal Fourier modes. Figure 11 shows the time-averaged variance spectrum of this tracer as a function of the horizontal wavenumber *K* = *k*^{2} + ℓ^{2} and vertical wavenumber *λ* (left panel) and as a function of the three-dimensional wavenumber *μ* = *k*^{2} + ℓ^{2} + *λ*^{2} (right panel). The wavenumbers are all nondimensionalized by 2*π*/*L*_{0}, where *L*_{0} is the horizontal domain size 1500 km; therefore, wavenumber 1 just fills the domain. The white curves in the left panel are isolines of *μ*, and the fact that the isolines of the variance are close to these curves suggests consistency with the prediction of an isotropic spectrum in the stretched space. The one-dimensional spectrum in the right panel follows a −1 power law (the thin solid line has a −1 slope for comparison), also consistent with isotropy. Notice that the spectra become very noisy at the largest vertical wavenumbers. This happens because there is no dissipation on vertical motions and because the finite-difference computation of the high-vertical modes for the nonconstant stratification are inaccurate.

The same two diagnostics used in Fig. 11 for the tracer spectrum are shown for the total energy *μ*^{2}|*ψ̂*_{kℓλ}|^{2} (nondimensionalized by the square of the maximum zonal speed *U*_{0} = 0.012 m s^{−1}) in Fig. 12. The solid line in the right panel has a slope of −3, showing that the total energy has the slope predicted by Charney’s (1971) theory for the forward cascade of potential enstrophy (see also McWilliams 1989). The spectrum peaks strongly at low horizontal and vertical wavenumbers, especially at the first baroclinic mode, consistent with geostrophic turbulence in the presence of surface-intensified stratification (Fu and Flierl 1980; Smith and Vallis 2001). The peak defines the scale and energy of the large, baroclinic eddies that stir the tracer. The low modes, however, dominate only in the spectra of energy; the potential enstrophy and tracer, in contrast, have a great deal of variance at high wavenumbers and exhibit three-dimensional isotropy in the vertically stretched coordinate space (the tracer shown in Fig. 11 and the potential enstrophy spectrum follow immediately from multiplication of the energy by *μ*^{2}).

### b. The effect of eddy velocity structure on the tracer filaments

A second view of the three-dimensionally stirred tracer arises from the analysis of Haynes and Anglade (1997), who considered the evolution of a tracer under the action of kinematic horizontal strain and vertical shear, with applications to stratospheric mixing. Haynes and Anglade explored the effects of both a steady, linear strain and shear, and a random-straining model. The core result follows from the evolution of a tracer blob under the action of a steady velocity field with strain and shear—that is, **u** = (*σx*, −*σy* + Λ*z*, 0). The development of tracer gradients in such a flow can be easily computed, as illustrated in appendix C. At long time scales, the squared aspect ratio (inverse slope) of the tracer obeying this solution is

In other words, the aspect ratio of tracer filaments follows the aspect ratio of the velocity field. Haynes and Anglade (1997) make additional, more general arguments, based on a more complex range of assumptions about the velocity field (i.e., a random-straining model with assumptions about the decorrelation times of the strain and shear), but the core prediction does not change much.

Haynes and Anglade’s (1997) scaling for *α* should apply to any velocity field with horizontal strain and vertical shear (so long as they do not decorrelate in time too quickly), whereas our argument above is specific to a flow generated from a forward potential enstrophy cascade. For the two arguments to be consistent, the velocity field in the simulations presented here must have an *N*/*f* aspect ratio. Figure 5, discussed previously, demonstrates the ample production of shear by the geostrophic cascade. In Fig. 13 we compare *N* ^{2}/*f* ^{2} (solid line), the aspect ratio of the tracer (dashed–dotted line), and the aspect ratio of the velocity field /*σ*^{2} (dashed line). [Here *σ*^{2} = is the horizontal strain.] To a strong degree, the tracer and velocity aspect ratios are close to *N*/*f* at all depths.

The aspect ratios averaged over the domain are a bulk measure of the flow statistics, but they do not give any information about how horizontal and vertical gradients of velocity and tracer vary with each another. Does the aspect ratio scaling hold at both strong and weak gradients, or is it dominated by a few strong gradients? In Fig. 14 we show the 2D probability density functions (PDFs) of tracer gradient (top panel) and velocity gradient (bottom panel) at a depth of 1100 m. Also plotted is the line with slope *N*(1100 m)/*f*. The results indicate that both the velocity field and tracer maintain a nearly constant aspect ratio, over a wide range of magnitudes of their respective vertical and horizontal gradients. Notice, however, that there is substantial scatter about the mean *N*/*f* value (discussed further in the next section).

It is apparently the case that low modes dominate only in the spectra of energy; shear and strain, like the potential enstrophy and tracer variance, must have a great deal of variance at high wavenumbers. To verify this claim, vertical wavenumber spectra of horizontally and time-averaged vertical shear and 100 × horizontal strain are shown in Fig. 15. A reference line with slope −1/2 is shown for reference, indicating that both spectra are very flat. The noisiness at high vertical wavenumbers is even more pronounced for shear and strain than for the tracer because, apart from the lack of direct dissipation of vertical motions, the velocity field is forced at high vertical mode (the mean shear drives baroclinic instability at a range of vertical modes), and there are extra derivatives involved. Therefore, the spectra are plotted only up to *λ*_{40} = 220 (where *λ _{n}* is defined in section 3b). Regardless, it is clear that the geostrophic cascade can produce ample small-scale gradients in the velocity field.

The structure of the tracer at submesoscales is only moderately affected by the energy and scale of the mesoscale eddies and instead tends toward a universal aspect ratio *N*/*f*. This implies that the shear/strain ratio in QG is dominated by small-scale velocity gradients rather than by the large-scale ratio of APE to KE, or in other words, by the aspect ratio of the energy-containing eddies. We have partially investigated whether this is the case as follows: we ran a series of tracer advection calculations in which we retained only the first few vertical modes of the velocity field from the control simulation. As the number of vertical modes was decreased, the strain/shear ratio approached the large-scale KE/APE ratio and so did the tracer aspect ratio, but they were both much larger than *N*/*f*. This strongly suggests that the shear/strain ratio and the associated tracer slope are set by the submesoscale gradients and not by gradients at the energy-containing scales. A more thorough test of our claim would be to run a series of simulations in which the large eddy scale is varied by changing the bottom drag. This analysis is ongoing and will be reported elsewhere.

### c. The effect of mean gradients on the tracer slope

The idealized passive tracer *ϕ* used to produce the spectrum in Fig. 11 is stirred by a mean horizontal gradient with no vertical structure: **Γ**_{ϕ} = **x̂**. Because **u**′ is dominated by kinetic energy at large horizontal and vertical scales in QG models, the eddy stirring of the large-scale tracer gradient, **u**′ · **Γ*** _{ϕ}* in the equivalent of (3.4) with

*C*′ replaced by

*ϕ*′, generates only tracer structure at large scales. Small-scale

*ϕ*fluctuations are primarily the result of the turbulent cascade of variance toward smaller scales, as we assumed throughout our analysis.

Adding vertical structure to the background gradient of a tracer will introduce a new source of tracer variance at small scales, and in this case, the vertical cascade of tracer variance should not display an inertial range. Consider, for example, a mean gradient for *ϕ̃* that has vertical structure = Γ(*z*)**x̂**. Because the advection equation is linear, we may write it as

where we have ignored dissipation. Defining a rescaled tracer *ϕ* = *ϕ̃*/Γ, the variance of *ϕ* is again forced only by the large scales of *u*′, and *ϕ* filaments approach an *N*/*f* aspect ratio. The squared aspect ratio of the rescaled tracer can then be compared to that of the original tracer, *ϕ̃*, and the result is

The aspect ratios of filaments of *ϕ̃* are thus modified by a factor proportional to |1 + *h*_{filament}/*h*_{background}|, where *h*_{filament} is the thickness of *ϕ* filament and *h*_{background} is the characteristic scale over which the background gradient Γ varies. For a tracer with a smooth background gradient |1 + *h*_{filament}/*h*_{background}| ≈ 1, and the tracer aspect ratio approaches *N*/*f*. However, if *h*_{filament} ≥ *h*_{background}, then the aspect ratios of filaments of *ϕ̃* are larger than those of the ideal (rescaled) tracer *ϕ*, because the vertical derivative of the tracer is augmented by the structure of the background.

In our simulations, *T* and *S* are examples of tracers forced with **Γ** = **Γ**(*z*); lateral tracer gradients are larger at the Mediterranean salinity tongue than above or below. The aspect ratio of temperature filaments, shown in Fig. 16, is indeed larger than that of *ϕ*, especially in the thermocline where the vertical variations in **Γ*** _{T}* are largest (a similar result holds for the salinity). The quasigeostrophic PV is another example of a tracer that is forced with a depth-dependent gradient,

**Γ**

*. Hence, we do not expect QGPV filaments to have an aspect ratio of*

_{q}*N*/

*f*everywhere. In regions where

**Γ**

*has a rich vertical structure, it is difficult to predict the slope of filaments, because the structure of the background gradient affects the dynamics, in addition to augmenting the small-scale QGPV gradients. A full theory for the generation of QGPV gradients in environments with variable*

_{q}**Γ**

*is beyond the scope of this paper.*

_{q}## 6. Vertical mixing and the end of the tracer cascade

The simultaneous enhancement of vertical and horizontal gradients in the tracer field (Fig. 11 and Fig. 14) suggests multiple paths to the dissipation of variance. In the mesoscale stirring scenario, the steady-state temperature variance level is determined by a balance between mesoscale production and small-scale dissipation:

where *κ _{H}* and

*κ*are horizontal and vertical turbulent diffusivities, respectively. Depending on the ratio of

_{V}*κ*and

_{H}*κ*, one can envision situations in which the tracer variance is dissipated by horizontal or vertical mixing.

_{V}The traditional assumption is that small-scale turbulent mixing is associated with breaking internal gravity waves at scales of meters (Gregg 1989). Wave breaking is believed to occur in isotropic patches, hence *κ* = *κ _{H}* =

*κ*. If, as suggested in the previous section, temperature gradients maintain the ratio =

_{V}*α*

^{2}, then the variance budget is

as originally posited by Haynes and Anglade (1997) for stratospheric tracer concentrations. Vertical diffusion clearly dominates if the small-scale turbulent diffusivity is isotropic, because *α* ≫ 1. One therefore expects that the tracer cascade can be arrested even in the total absence of horizontal diffusion. Here we explore the effectiveness of vertical diffusion in halting the laterally driven variance cascade.

To investigate the effects of vertical and horizontal diffusion on a tracer, we ran a small suite of additional simulations at 512^{2} horizontal resolution, using the same stratification as for the central simulation, but with a simpler velocity structure obtained by setting the mean meridional component to 0, and replacing the zonal component *u* with its projection on the first baroclinic mode. These runs were also computed in a larger domain, of size 3000 km × 3000 km (these runs were originally purposed to investigate the effects of eddy scale on the tracer slope, which is the only reason the domain is larger, and this has no bearing on the results discussed here). Using an ideal tracer *ϕ* with unit mean gradient in the *x* direction (as in the last section), we replace the exponential cutoff filter with either a horizontal laplacian diffusion or a vertical laplacian diffusion. Figure 17 shows horizontal and vertical slices of the tracer in a simulation with only horizontal diffusion (left panels) and in a simulation with significant vertical diffusion (right panels). One finds that vertical diffusion is not only effective in damping vertical variance but also in smoothing the horizontal field. A closer inspection of the horizontal fields, however, reveals that there is a wider range of filament widths in the case with vertical diffusion than in that with only horizontal diffusion. Likewise, although horizontal diffusion apparently halts the vertical cascade of tracer variance, there is a wider range of vertical scales in that case than in the case with vertical diffusion.

The horizontal structure of the two cases is investigated as follows. In Fig. 18 we plot the horizontal spectra of tracer variance, at 1100-m depth, for four simulations: one with only a horizontal filter, one with a horizontal laplacian diffusion, and two with only vertical diffusion (see figure caption and inset for details). A few features are apparent. First, moving from a filter to laplacian diffusion in the horizontal primarily changes the scale of the roll-off; moreover, although broader for the laplacian diffusion, a roll-off scale can be readily identified in both cases. The roll-off scale is the typical filament width, given by the Batchelor scale for the effective horizontal diffusion *κ _{H}*/Σ, where Σ is the stretching rate experienced by the filament. Although vertical diffusion is effective at absorbing laterally created variance, it apparently does so in a far less scale-selective manner than horizontal diffusion of any type; a roll-off scale cannot be clearly defined.

The lack of a clear horizontal filament width in the vertical diffusion cases can be understood as follows. Considering again the right panel of Fig. 14, one sees that the typical aspect ratio of a tracer filament is *N*/*f*, but in fact filaments take on a wide range of ratios *α* ± *δα*. Hence, diffusion sets the vertical scale of the filaments to be *κ*/Σ, whereas the horizontal Batchelor scale takes a wide range of values given by

Those filaments with very shallow slopes will experience very little horizontal diffusion, whereas those with steeper slopes experience quite enhanced horizontal diffusion. The result is that the arrest of the cascade spans a wide range of scales, and thus no one filament width can be identified. Haynes and Vanneste (2004) reached a similar conclusion in their study: the lack of a clear roll-off scale in simulated stratospheric tracer concentrations was attributed to variability in the vertical shear, hence leading to variability in *α* and, ultimately, to the cutoff wavenumber .

Estimating the filament width from (6.2) is quite delicate, because one must calculate the mean stretching rate Σ experienced by filaments. Sundermeyer and Price (1998) and Hu and Pierrehumbert (2001) show that the stretching rate is typically much smaller than the rms strain rate *σ*, because *σ* is dominated by regions of large strain, whereas filaments tend to be expelled rapidly from high-strain regions. Techniques have been developed to estimate Σ (e.g., Pierrehumbert 2000), but for our purposes it is far simpler to estimate the filament widths from the tracer spectra. We then use (6.2) to extrapolate our results, obtained using two large values of *κ _{V}* to the smaller

*κ*values found in the ocean thermocline.

_{V}Consider the spectra in Fig. 18 for the two simulations where the cascade of tracer variance is arrested by a vertical diffusivity. We can estimate the largest possible filament width as the inverse of the wavenumber at which the spectrum deviates from that of the run with a horizontal filter. For the simulation with *κ _{V}* = 1.7 × 10

^{−2}m

^{2}s

^{−1}, this is about

*K*= 10, or a filament width of 3000 km/(2

*π*)/10 ∼ 50 km. For the simulation with

*κ*= 1.7 × 10

_{V}^{−3}m

^{2}s

^{−1}, this is about

*K*= 30, or a filament width of ∼16 km. The decrease in scale is proportional to the square root of the ratio of the diffusivities, as expected from the Batchelor scaling in (6.2). Extrapolating the scaling to a diffusivity of order

*κ*= 10

_{V}^{−5}m

^{2}s

^{−1}gives filament widths as large as about 3 km. Ledwell et al. (1998) also estimated tracer filament widths of a few kilometers from a tracer release experiment done as part of NATRE. However, the comparison between simulations and data should not be taken too seriously, because the set of runs used for this calculation was not designed to match the conditions encountered in NATRE; the simulations were run with only first baroclinic shear, a larger domain size, and a lower resolution than the central simulation.

Finally, one should compare the dissipation of thermal variance by geostrophic stirring versus that induced by internal wave shear. Young et al. (1982) show that the coupling between a wave with vertical shear and turbulent mixing results in an enhancement of diffusion in the horizontal, a result known as shear dispersion. They find that for the Garrett–Munk internal wave spectrum (Garrett and Munk 1979), which describes the internal wave continuum in the ocean away from boundaries, the enhancement is given by *κ _{H}* =

*κ*(

*N*/

*f*)

^{2}.

This results in a dissipation of thermal variance proportional to *χ*_{IW} = *κ*(*N*/*f* )^{2}*a* and filament widths ℓ_{filament} = (*κ*/Σ)^{1/2}*N*/*f*. At first glance, these scalings appear to be the same as those derived for the generation and dissipation of tracer filaments by geostrophic stirring because *α* ≈ *N*/*f*. However, there is a fundamental difference in the two problems. The enhancements of lateral diffusion by internal wave shear shows little departure from the scaling *κ _{H}* =

*κ*(

*N*/

*f*)

^{2}, because the Garrett–Munk spectrum is quite stable throughout the ocean. The scatter about

*κ*=

_{H}*κ*(

*N*/

*f*)

^{2}is instead very large in the geostrophic stirring problem because the PDF of tracer slopes is broad. Ledwell et al. (1998) refuted the hypothesis that internal wave shear dispersion set the width of the tracer filaments observed in NATRE because ℓ

_{filament}= (

*κ*/Σ)

^{1/2 }

*N*/

*f*≈ 400 m was nearly an order of magnitude smaller than observed. (In NATRE at the tracer release level,

*κ*≈ 1 × 10

^{−5}m

^{2}s

^{−1}, Σ ≈ 3 × 10

^{−7}s

^{−1}, and

*N*/

*f*≈ 60.) However, our estimates of tracer widths from the simulations show that in the geostrophic stirring problem, the largest filaments can be nearly an order of magnitude larger than the mean value. A strong observational test of the scenario described in this paper would be to compute the distribution of filament widths in a tracer release experiment on a wide range of scales from

*O*(1 km) and

*O*(100 m). If the distribution is wide, our results would be vindicated.

## 7. Summary and discussion

The main result of this study is that the thermal variance dissipated at molecular scales in the NATRE region below the main thermocline is generated through lateral eddy stirring of the large-scale thermal gradient. This was demonstrated by running a QG model forced with the observed large-scale shear, vertical stratification, and mean lateral *T*–*S* gradients. A vigorous geostrophic eddy field developed through baroclinic instability of the mean flow. The eddies stretched and twisted the climatological *T*–*S* gradients and produced sharp *T*–*S* filaments along-isopycnals. The eddy generation of *T*–*S* variance was particularly effective at the level where salty and warm Mediterranean Waters overflow into the fresher and colder Intermediate Atlantic Waters, resulting in large thermohaline contrasts. The thermal variance generated through eddy stirring was close to in situ measurements of thermal dissipation at molecular scales. Our work suggests that in regions characterized by large watermass contrasts, the generation of compensated *T*–*S* gradients by lateral eddy stirring can overcome the generation of *T*–*S* variance by vertical turbulence.

The QG model was used to study the cascade of *T*–*S* variance from climatological to dissipation scales. Following previous work (Haynes and Anglade 1997; Klein et al. 1998), we showed that geostrophic eddy stirring generates both horizontal and vertical gradients, that is, the *T*–*S* variance is cascaded both to small horizontal *and vertical* scales. The horizontal cascade is driven by lateral strain, which compresses the tracer isopleths. The vertical cascade is due to the tilting of the tracer filaments by geostrophic eddy shear. Haynes and Anglade (1997) speculated that the coupling between the two cascades would result in tracer filaments with a slope proportional to the ratio of the strain and shear in the first few vertical modes, which contain the bulk of the energy in QG simulations (and the real ocean). However, we found that tracer filaments are generated by small-scale strain and shear. At small scales, the shear–to–strain ratio is proportional to *N*/*f* (Charney 1971). Hence, geostrophic stirring generates vertical gradients that are *N*/*f* times larger than the horizontal gradients. It seems to have gone unnoticed that eddy stirring is so effective at generating vertical gradients. In addition, this implies that the slope of tracer filaments is much steeper than typical isopycnal slopes. Nandi et al. (2004) report compensated *T*–*S* filaments with a slope steeper than the isopycnal slope in seismic profiles collected in the Norwegian Sea, consistent with our results.

Next we addressed the question of what mixing process arrests the formation of compensated *T*–*S* gradients. The traditional explanation is that there is some small-scale, three-dimensional instability, involving breaking gravity waves, that halts the formation of gradients before the molecular scale. We showed that vertical mixing, in the form of a constant turbulent diffusivity *κ _{V}*, is sufficient to arrest the cascade of

*T*–

*S*variance. The dissipation of thermal variance is through the mixing of vertical tracer gradients, so that

*χ*

_{meso}=

*κ*; however, it can also be expressed in terms of horizontal gradients because of the

_{V }*N*/

*f*aspect ratio of tracer filaments, so that

*χ*

_{meso}≃

*κ*

_{V}f^{−2}

*N*

^{2}(where

**∇**

*T*includes only horizontal derivatives in our notation). One might be tempted to conclude that eddy stirring of tracers in the ocean can, therefore, be rationalized as a horizontal advection–diffusion problem with an enhanced horizontal diffusivity

*κ*≃

_{H}*κ*

_{V}f^{−2}

*N*

^{2}(Haynes and Anglade 1997). This is misleading, because a horizontal diffusivity is scale selective (it sets the horizontal scale of tracer filaments). Vertical mixing is not scale selective in the horizontal. Tracer filaments have a mean aspect ratio of

*N*/

*f*, but there is substantial scatter around the mean value. Hence, the vertical scale at which mixing arrests the tracer cascade corresponds to a wide range of horizontal scales. Numerical simulations confirm that tracer filaments have a well-defined vertical scale but no characteristic horizontal scale.

Ledwell et al. (1998) estimated that the width of tracer filaments in the NATRE region at 300 m was 2–5 km. Our model predicts that the width of tracer filaments should be in the wider range of 700 m–7 km at the same level, if we were to use a vertical diffusivity of *κ _{V}* = 1 × 10

^{−5}m

^{2}s

^{−1}. It is hard to say whether the two results are consistent or not. The horizontal resolution of the observations was limited to scales larger that

*O*(1) km, whereas the theoretical prediction here is that the distribution of possible filament widths should be quite wide, spanning an order of magnitude.

Although our theory seems consistent with existing data, if new observations showed that the range of filament widths is narrow banded, we would speculate on two alternative scenarios that might then be invoked. Sundermeyer et al. (2005) and Polzin and Ferrari (2004) suggest that mixing by gravity wave breaking is not isotropic. The idea is that wave breaking occurs in isotropic patches, but the patches rapidly undergo a geostrophic adjustment and become flat pancakes, referred to as *vortical modes* [see Riley and de Bruyn Kops (2003) for a recent, detailed investigation of stratified turbulence]. Sundermeyer et al. (2005) speculate that vortical modes have lateral scales of a few kilometers and could act as efficient horizontal mixers. In terms of the variance budget in (6.1), Sundermeyer et al. (2005) claim that *κ _{H}* ≫

*α*

^{2}

*κ*. This would further imply that the cascade of tracer variance in the ocean is arrested at scales much larger than the characteristic scale of wave breaking; that is, tracers should become smooth at scales of a hundred of meters instead of tens of meters. Observational evidence of vortical modes has so far remained elusive, and the scenario must be treated as an intriguing, but unproven, hypothesis.

_{V}Alternatively, one could argue that vertical mixing is not *passive*, in the sense that it is not independent of the *T*–*S* variability. As *T*–*S* fronts develop, they might become unstable to secondary instabilities that trigger mixing and limit the further narrowing of the front. In other words, the diffusivity *κ _{V}* might not be a constant and would tend to increase at fronts. This might be termed an

*active*process, in that it produces diapycnal mixing that would not occur in the absence of lateral stirring. There are three (and possibly more) plausible mechanisms for this conclusion.

Double-diffusive instabilities. Lateral eddy stirring generates compensated

*T*–*S*fonts with density ratios close to one—that is,*R*≡*αT*/_{z}*βS*≈ 1. Laboratory experiments suggest that at fronts with_{z}*R*≈ 1, double-diffusive instabilities grow quickly and drive diapycnal diffusivities up to 10 times larger than background values set by internal wave turbulence (Schmitt 1994). Notice that this scenario is quite different from the traditional view that double-diffusive instabilities are triggered by climatological*T*–*S*gradients. Here the instabilities appear along the*T*–*S*filaments generated by eddy stirring.Straining of the internal wave field. Bühler and McIntyre (2005) show that the mesoscale strain field will squeeze internal wave packets, leading to a decrease in their wavenumbers and eventually to wave breaking. The effect is concentrated in the same straining regions where sharp tracer filaments are generated. Hence, there might be a correlation between wave breaking and the generation of tracer filaments.

Loss of balance. Molemaker et al. (2009) show that sharp temperature gradients can undergo a loss of balance in primitive equation simulations. The difference from QG simulations is that kinetic energy begins to dominate over potential energy at small scales, Richardson numbers become small, and Kelvin–Helmholtz instabilities lead to irreversible mixing, thus the prediction that there is a collocation between high mixing and

*T*and*S*fronts.

Whether our hypothesis will prove relevant depends on the efficiency of the described processes. Double-diffusive instabilities can be substantially damped by background internal wave or geostrophic shears (Kunze 1994; Beal 2007). Similarly, it is unclear whether mesoscale strain modulation of the internal wave filed can trigger substantial wave breaking. Finally, Ford et al. (2000) argue that the loss of balance is strongly inhibited for typical oceanographic stratification and velocities. Clearly there is need for a more thorough investigation of these dynamics.

In this work we focused on the cascade of *T*–*S* variance to small scales. However, our arguments apply to potential vorticity as well. Eddy stirring cascades potential enstrophy (potential vorticity variance) to large horizontal and vertical wavenumbers. In the active mixing scenario, the small-scale potential vorticity filaments would be associated with enhanced mixing, suggesting the intriguing possibility that vertical mixing in the ocean is coupled to the potential enstrophy cascade. Such a feedback would have important implications for our understanding of the ocean system.

## Acknowledgments

We acknowledge useful discussions with our colleagues Patrice Klein, Alan Plumb, Kurt Polzin, and two anonymous referees. We thank Jim Ledwell for permission to reproduce Fig. 2. KSS was supported by NSF Grant OCE-0620874 and ONR Grant N0001409010633, and RF was supported by NSF Grant OCE-0825376 and ONR Grant N000140910458.

## REFERENCES

**,**(

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Modeling Temperature and Salinity Advection in Quasigeostrophic Flow

We decompose temperature and salinity as

where the asterisk superscripted variables are the background profiles, the primed variables are the dynamical deviations, and **Γ*** _{T}* and

**Γ**

*are horizontally constant approximations of the mean horizontal gradients*

_{S}**∇**

*T*and

**∇**

*S*, respectively. The conservation equations for each are

and

We will formally assume an eddy diffusion and therefore set *κ* ≡ *κ _{T}* =

*κ*. Adding and subtracting the advection–diffusion equations for temperature and salinity (multiplied by the thermal and haline contraction coefficients

_{S}*β*and

_{T}*β*, respectively) and using (3.1), we find

_{S}where

The linear gradients are

where thermal wind balance has been used to relate the horizontal density gradient to the vertical shear of the mean horizontal velocity. Elimination of *w*′ between the density and spice equations leads to the advection–diffusion equation for *C*:

[stated in Eq. (3.4) without the diffusion term], where *C*′ = *γ*′ − (*N _{γ}*

^{2}/

*N*

^{2})

*ρ*′ and

Recall also that *N* ^{2}, *N _{γ}*

^{2},

**Γ**

_{γ}, and

**u**= (

*u*, ) are functions only of

*z*.

The tracer can be represented as *C*′ = Γ* _{C}^{x}A*′ + Γ

*′, where*

_{C}^{y}B*A*′ and

*B*′ are independent tracers stirred separately by the zonal and meridional flow, solving

respectively, and **Γ*** _{C}* = (Γ

*, Γ*

_{C}^{x}*). In practice, we simulate*

_{C}^{y}*A*′ and

*B*′, and from these we reconstruct

*C*′, given the mean gradient

**Γ**

*.*

_{C}From the definition of *C* in (3.3), the quasigeostrophic relation of density to streamfunction *ρ*′ = −*f*_{0}*g*^{−1}∂* _{z}ψ*′ and the relation of temperature and salinity to spice and buoyancy, we find

and

### APPENDIX B

#### Using Observations to Generate the Mean Fields

The HRP data gives mean background values and linear gradients for the potential temperature *T* and salinity *S*. From these the UNESCO equation of state is used to generate the background profile and mean linear gradient of density, *ρ**(*z*) and **Γ**_{ρ}, respectively. A linear equation of state is also defined from the data, with *T*_{0} = 6 C, *S*_{0} = 35.2 psu, *β _{T}* = 1.25 × 10

^{−4}K

^{−1},

*β*= 7.69 × 10

_{S}^{−4}psu

^{−1}, and

*ρ*

_{0}= 1028 kg m

^{−3}—chosen to best match the potential density measured at the depth of the Mediterranean salinity tongue (about 1100 m). The linear equation of state coefficients are used to define the background and mean linear gradient of spice,

*γ**(

*z*) and

**Γ**

_{γ}, respectively. The mean velocity field

**u**is computed from

**Γ**

_{ρ}via thermal wind balance in (A.1) and with a reference level at 1500 m given by current meter data, and background stratification profiles

*N*

^{2}and

*N*

_{γ}^{2}computed from vertical derivatives of

*ρ**(

*z*) and

*γ**(

*z*), respectively. Modeled temperature and salinity fields are reconstructed from the model tracers using the linear equation of state.

### APPENDIX C

#### Slope of Tracer Filaments in QG Turbulence

Consider the conservation equation for the passive tracer *C* defined in (3.3). This tracer is advected by the two-dimensional geostrophic velocity and not by the vertical velocity: in QG vertical advection acts only on the background gradients and *C* is defined to have no such background gradient. In this appendix, we discuss what sets the slope of *C* filaments generated by the geostrophic stirring velocity.

Tracer gradients are generated by geostrophic stirring at a rate determined by the equations

where *σ _{ij}* ≡ ∂

*u*/∂

_{i}*x*is the horizontal strain tensor and Λ

_{j}*≡ ∂*

_{i}*is the vertical shear vector. For simplicity we neglect tracer dissipation, which is assumed to act on time scales longer than the fast time over which filaments are generated by stirring. The first equation shows that lateral tracer gradients are generated by strain in the horizontal velocity field. The second equation in Eq. (C.1) shows that even though advection is horizontal, vertical gradients are generated by the vertical shear tilting over the horizontal tracer filaments.*

_{z}u_{i}Under the assumption that the velocity field is slowly varying along a Lagrangian trajectory (*D σ*/

*Dt*∼ 0 and

*D*

**Λ**/

*Dt*∼ 0), the tracer gradient equations can be integrated yielding

where *σ* = ±−det** σ**. Whether tracer gradients increase along Lagrangian trajectories therefore depends on the sign of the determinant of

**. In regions where det**

*σ***> 0, then**

*σ**σ*is imaginary, and tracer gradients are not actively generated in these regions. However, in regions where det

**< 0, then**

*σ**σ*is real, and the flow is locally hyperbolic. Tracer contours are compressed by the strain leading to an exponential increase of the gradients. These are the regions where tracer filaments develop. According to (C.2) the aspect ratio of the tracer filaments, defined by

*C*=

_{z}

*α*∇*C*, is given by

The prediction of this simple analysis is that the aspect ratio of a tracer filament (inverse slope) is given by the ratio of the shear acting across the tracer gradient and the strain rate (Haynes and Anglade 1997).

One final comment pertains to the assumption that the strain and shear are slowly varying along a Lagrangian trajectory. Klein et al. (1998) show that if one relaxes the approximation, new terms related to accelerations/decelerations of the velocity field appear in the equations for the tracer gradients in (C.1). These terms account for the fact that tracer filaments cannot lie indefinitely in regions of maximum strain and lead to a reduction in the growth rate of tracer filaments, *σ*. We find that these effects lead to quantitatively appreciable variations in tracer slope, but the scaling for the tracer aspect ratio—that is, *α* ≈ *N*/*f*—is qualitatively very robust.

## Footnotes

*Corresponding author address:* K. Shafer Smith, Courant Institute of Mathematical Sciences, New York University, 251 Mercer St., New York, NY 10012. Email: shafer@cims.nyu.edu

^{1}

The dissipation operator *D _{C}* becomes quite convoluted if the operators

*D*and

_{ρ}*D*involve vertical derivatives and the background stratification

_{γ}*N*

^{2}and

*N*

_{γ}^{2}are depth dependent. Therefore, we will not use the tracer

*C*in section 6 where we consider dissipation operators

*D*and

_{ρ}*D*with vertical derivatives.

_{γ}^{2}

The “strain” rate reported by Ledwell et al. (1998) was computed from a direct estimate of the rate of growth of an individual filament in the along-streak direction. This filamentary growth rate may in fact be significantly less than the rms strain rate computed from the velocity field tensor (the quantity plotted in the right panel of Fig. 5), as found by Sundermeyer and Price (1998).