## Abstract

Thermohaline interleaving is an important mechanism for laterally fluxing salt, heat, and nutrients between water masses. Interleaving is driven by a release of potential energy resulting from the differing diffusivities of heat and salt in seawater. The flows are composed of stacked intrusions that flux more and less buoyant water in opposite directions. In this paper, the role of shear instability caused by this juxtaposed motion is investigated. The model described in Walsh and Ruddick is modified to include both the effects of shear-induced turbulence and an improved convective mixing parameterization. Shear and convective mixing play a similar and significant role in interleaving dynamics. In the absence of either instability, cross-front fluxes are increased by approximately 30%. While in situ observations of horizontal diffusivity resulting from interleaving are not yet precise enough to calibrate the parameterizations independently, parameter values based on independent laboratory and numerical studies lead to diffusivity predictions that are within the error of the observations.

## 1. Introduction

Observational and theoretical studies confirm the importance of thermohaline interleaving in mixing water mass properties across fronts. For example, studies of Mediterranean Sea eddies (meddies; McDowell and Rossby 1978) identify thermohaline interleaving as the predominant mechanism of meddy decay (Armi et al. 1988; Hebert et al. 1990). In the Arctic, thermohaline intrusions are coherent over essentially the entire basin (Walsh and Carmack 2003). These intrusions may be an important mechanism for mixing (relatively) nutrient-rich Pacific Ocean water into the Canada Basin (McLaughlin et al. 2003), in addition to mixing Atlantic Ocean water into the Arctic. The goal of this paper is to present an improved theoretical description of interleaving, which will lead to an improved ability to predict mixing in frontal regions.

The term thermohaline indicates that the difference in the molecular diffusivities of heat and salt is responsible for the interleaving instability. The nonuniform distribution of heat and salt within the ocean, combined with the vastly different diffusivities of heat and salt, allows these structures to form. To understand the mechanism, imagine a front separating two water masses of different temperatures and salinities where natural variations in the ocean cause fluctuations in the front. Temperature in the fluctuation will adjust to the temperature of the surrounding water more quickly than salinity will, altering the density of the perturbations and causing buoyancy-driven flow. In the ocean, thermohaline interleaving is not driven directly by molecular diffusion but rather by buoyancy fluxes resulting from microscale mixing processes, such as salt fingering and diffusive convection. Horizontal motions across the front affect the vertical temperature and salinity gradients that govern microscale fluxes such that the motion is reinforced (Kunze 2003; Ruddick 2003).

After the initial formation of the interleaving structure, shear instability and convective overturning are introduced by flux divergences that result from differences in the way in which salt fingering and diffusive convection transport heat and salt. Shear instability occurs when the stabilizing effect of the buoyancy gradient is overcome by the destabilizing effect of shear. While shear instability is affected both by vorticity and the buoyancy gradient, convective instability depends solely on the degree to which the negative buoyancy gradient overpowers molecular diffusivities. Parameterization of these mixing processes in the context of interleaving is the primary aim of this paper.

Our work builds upon previous studies, beginning with Stern (1967), who hypothesized that the temperature and salinity inversions observed along fronts of strong lateral temperature and salinity gradients could be caused by a double-diffusive instability. Stern’s linear theory posited a uniform eddy diffusivity for salt and related the double-diffusive thermal and saline buoyancy fluxes, *F _{T}* and

*F*, by the uniform buoyancy flux ratio

_{S}*γ*, such that

_{f}*F*=

_{T}*γ*

_{f }*F*. His model validated the hypothesis, although it predicted infinite growth at small wavenumbers, a problem that was resolved later by the inclusion of viscosity (Toole and Georgi 1981). Walsh and Ruddick (1995) introduced a dependence on the density ratio

_{S}*R*, allowing diapycnal diffusivities to vary with stratification.

_{ρ}The nonlinear initial value problem was first addressed by Walsh and Ruddick (1998, hereinafter WR98), who assumed that interleaving motions are unidirectional, thus casting the problem into the numerically efficient one-dimensional (1D) form. In the WR98 model, density gradients between layers of salt fingering and diffusive convection are smoothed by the action of background turbulence (representing random internal wave breaking) and convective overturning. Background turbulence is applied uniformly, while convective turbulence is activated only in the regions where density is unstably stratified. In the convectively unstable regions, a predetermined convective diffusivity turns on and off without accounting for the degree of static instability. For a more comprehensive overview of the theory and laboratory studies, and observations of thermohaline interleaving, the reader is referred to Ruddick and Kerr (2003), Ruddick (2003), and Ruddick and Richards (2003).

In this paper, we explore interleaving in a 1D, nonlinear model based on WR98. Our goal is to advance the theoretical description of thermohaline interleaving so as to improve the ability to predict mixing in frontal regions. To this end, we employ a more advanced convective turbulence parameterization, and we introduce a parameterization of turbulent fluxes resulting from the shearing motion of the intrusions.

Interleaving scales range from *O*(10^{−3} m) for microscale mixing to *O*(10^{6} m) for the macrostructure of the front. The interleaving process is intrinsically 3D; however, the computational challenges for modeling all scales of motion in 3D are currently insurmountable. The main advantage of 1D models is the reduced computational cost, which renders them easier to apply to multiple scenarios. They are limited, however, in their ability to capture all the complexity of interleaving dynamics, for example, the dynamics in the nose of intrusions and feedback mechanisms between intrusion mixing and the background environment. Baroclinic shear can be included in a 1D model, but is not considered here. In this paper, we simplify the interleaving processes in order to isolate the effects of shear and convective mixing. To this end, a 1D model is sufficient.

Methodology, including parameterizations of shear instability and convective overturning, is covered in section 2. In section 3, we describe results for the base case, in which environmental parameters are set to model a meddy and flux parameters are given their “best” values, based on laboratory experiments and direct numerical simulations (DNS). In section 4, we vary the flux parameter values in order to better understand the effects of the various processes and to assess uncertainty in the model predictions. Conclusions are summarized in section 5.

## 2. Method

In this section, we review the model developed by WR98 and Merryfield (2000) and describe the new parameterizations used in the present model to represent the effects of shear and convective instabilities.

Thermohaline fronts may be classified by their vertical stratification. When warmer, saltier water overlies cooler, fresher water (e.g., the lower flank of warm-core eddies), the water column is unstable to salt fingering. In the reverse case (e.g., the polar regions), diffusive convection ensues. In both cases, temperature and salinity can compensate for one another in the horizontal and thus allow equilibrium without a thermal wind shear. While both environments support thermohaline interleaving, this paper focuses on a case where the background gradients are unstable to salt fingering.

### a. The one-dimensional model

A Cartesian coordinate system is used where *x* represents the across-front direction, *y* represents the alongfront direction, and *z* represents the vertical direction. The front is assumed to be much larger than the intrusion, allowing end effects to be neglected and initial property gradients to be approximated by uniform values. Variation in the *y* direction is neglected, and there is no flow in the alongfront direction, that is, ∂/∂*y* = 0, *υ* = 0. Buoyancy is defined as *b* = −*g*(*ρ* − *ρ*_{0})/*ρ*_{0}, where *g* is the acceleration resulting from gravity and *ρ*_{0} is a uniform reference density. Buoyancy is split into an initial, or background, component and a fluctuating component, whereas velocity has only a fluctuating component. Throughout this paper, the fluctuating velocity and buoyancy components will be denoted by lower-case letters, and background values by upper-case letters. The definitions of velocity and buoyancy are, therefore,

The unit vectors **î**, **ĵ**, and **k̂** represent the *x*, *y*, and *z* directions. The contribution to buoyancy by temperature or salinity is denoted with a superscript *T* or *S*. The front is assumed to be compensated in temperature and salinity with constant horizontal and vertical gradients, that is, *B*^{S}_{x}, *B*^{S}_{z}, *B*^{T}_{x}, *B*^{T}_{z} are constants and *B*^{S}_{x} + *B*^{T}_{x} = 0. In the fingering-favorable case on which we focus here, *B*^{T}_{z} > 0 and *B*^{S}_{z} < 0.

The Boussinesq equations describing the evolution of buoyancy and velocity are

Here *D*/*Dt* symbolizes the material derivative (*D*/*Dt* = ∂/∂*t* + **u** · **∇**), *π* = *p*/*ρ*_{0} is the scaled pressure, and *A* and *K _{i}* represent eddy viscosity and diffusivities. Equation (7) represents two equations for the evolution of the thermal and saline buoyancies, where “

*i*” symbolizes

*T*or

*S*. Note that parameterized microscale fluxes act in the vertical direction only.

The curl of (6) yields

where Ψ is a streamfunction, such that *u* = −∂Ψ/∂*z* and *w* = ∂Ψ/∂*x*, and *ω* = −∇^{2}Ψ is the vorticity. Expanding *b*^{i} in (7) into the background and perturbation components and replacing the velocity components with the streamfunction Ψ results in the following scalar equations:

Following WR98, the coordinate frame is now tilted by the intrusion slope *θ* such that *χ* represents the tilted *x* coordinate (along-intrusion direction) and *ξ* represents the tilted *z* coordinate (across-intrusion direction). The rotation operator is

where *C* = cos(*θ*) and *S* = sin(*θ*). The intrusion slope is typically small (<1°) and is calculated using linear theory (Walsh and Ruddick 1995). Perturbations vary only in the across-intrusion (*ξ*) direction, so that ∂/∂*χ* = 0 for perturbation quantities. The rotated equations are

where subscripts indicate partial derivatives. Equations (11)–(13) form a closed set of evolution equations for *ω*, *b*^{T}_{ξ}, and *b*^{S}_{ξ}.

Boundary conditions are periodic in *ξ*. Each model run was initialized using values of the intrusion slope, vertical wavenumber, and initial eigenfunction amplitudes derived from linear theory (Walsh and Ruddick 1995). For the cases described here, the vertical domain accommodates two wavelengths of the fastest-growing mode, which allows for the possibility of subharmonic instability.

All integrations and derivatives were performed in Fourier space using 256 grid points (128 Fourier modes) per vertical wavelength. Time was discretized using the third-order Adams–Bashforth method.

### b. Mixing parameterizations

The diffusivities *A*, *K _{T}*, and

*K*represent combinations of several parameterized microscale mixing processes. Diffusivities resulting from particular processes will be identified by a superscript, indicating salt fingering (sf), diffusive convection (dc), shear-driven turbulence (shr), convection (cnv), or ambient turbulence (at). Ambient turbulence represents the effect of random internal gravity wave breaking events, and is modeled by a constant diffusivity

_{S}*K*

^{at}= 1 × 10

^{−5}m

^{2}s

^{−1}that is applied equally to velocity and thermal and saline buoyancies. Molecular effects are also represented by constant diffusivities for velocity

*ν*= 1.0 × 10

^{−6}m

^{2}s

^{−1}, thermal buoyancy

*κ*=

_{T}*ν*/7, and saline buoyancy

*κ*=

_{S}*κ*/100. Diffusivities resulting from other microscale processes vary in space and time in accordance with parameterizations to be described in the remainder of this section.

_{T}#### 1) Double-diffusive mixing

Salt fingering and diffusive convection parameterizations are dependent on the instantaneous density ratio, defined by

When *R _{ρ}* is negative, the stratification is either stable in both temperature and salinity or unstable to convection; in either case, the salt fingering and diffusive convection instabilities are inactive. When 0 ≤

*R*< 1 the stratification is unstable in temperature and able to support diffusive convection. Salt fingering is characterized by

_{ρ}*R*≥ 1, in which case the stratification is unstable in salinity. The strength of both the salt fingering and the diffusive convection instabilities increase as

_{ρ}*R*approaches unity.

_{ρ}Following Merryfield and Grinder (2000, the diffusivities of temperature, salinity, and momentum resulting from salt fingering are defined by

where *γ _{f}* is the salt fingering flux ratio,

*τ*=

_{f}*κ*/

_{S}*κ*is the inverse Lewis number, and Sc

_{T}*=*

_{f}*A*/

*K*is the Schmidt number. Default values are

_{S}*γ*

*= 0.6,*

_{f}*τ*= 0.01, and Sc

_{f}*= 1. Note that (15) pertains strictly to two-dimensional salt fingers, but is expected to be appropriate for interleaving because salt fingers become two-dimensional in the presence of shear (Linden 1974).*

_{f}The diffusive convection parameterization follows Merryfield (2002) and Kelley (1990):

where Pr* _{d}* =

*A*/

*K*denotes the turbulent Prandtl number in a diffusive convection environment and has a default value of 1.0. Effects of varying both Sc

_{T}*and Pr*

_{f}*are considered in section 4d.*

_{d}#### 2) Shear-driven turbulence

Turbulence driven by the shear of the interleaving motions has not been included in previous modeling studies, so we describe its parameterization in detail. We employ a mixing length model with stratification represented via the Richardson number. This model is relatively simple and does not involve the solution of any additional equations, but it is accurate enough for the present purpose.

In the absence of stratification, an eddy viscosity resulting from sheared turbulence may be written in the nonlocal mixing length form *K*^{shr} = *a*_{0}*H*^{2}_{s}|*ω*|, where *H _{s}* is the vertical thickness of a sheared layer,

*a*

_{0}is a constant, and |

*ω*| is the magnitude of the shear

*U*. [In the interleaving context,

_{ξ}*U*(

*ξ*,

*t*) is the along-intrusion velocity.] Without forcing, such a layer tends to spread in time because of turbulent entrainment, a process that can be described by the diffusion equation

*U*= (

_{t}*K*

^{shr}

*U*)

_{ξ}*. A similarity solution in the form*

_{ξ}*U*=

*U*[

*ξ*/

*H*(

_{s}*t*)] results in

*dH*/

_{s}*dt*= (64/9)

*a*

_{0}Δ

*U*, where Δ

*U*is the velocity difference across the layer. Details may be found in section 5.4.2 of Pope (2000). This constant spreading agrees with laboratory experiments (Bell and Mehta 1990) and DNS (Rogers and Moser 1994), which yield

*dH*/

_{s}*dt*= ΣΔ

*u*with spreading rate Σ in the range 0.06 < Σ < 0.11. Equating these results, we find

*a*

_{0}in the range 0.0085–0.0155. The higher values in the DNS and laboratory experiments came from spatially developing turbulence downstream of a splitter plate, while the lower values pertain to temporal development between two oppositely directed flows. Because the latter case most closely resembles shear-driven mixing between interleaving layers, we anticipate that a value at the low end of the range will be most appropriate. We therefore choose the compromise value

*a*

_{0}= 0.01.

The diffusion equation used above to describe the spreading of an unforced shear layer is an incomplete model for interleaving flow, where buoyancy forces oppose shear layer spreading (and must ultimately halt it if an equilibrium interleaving state is to be attained). Nevertheless, the simpler case of the unforced shear layer embodies the physics of shear-driven mixing well enough, we believe, to provide a reasonable estimate of the parameter *a*_{0}. We test for sensitivity to the value of *a*_{0} in section 4a.

Shear-driven mixing between interleaving layers is expected to be modified by the action of density stratification, which may be either stable or unstable. The role of stratification is parameterized in terms of the Richardson number Ri = b* _{ξ}*/

*ω*

^{2}. When Ri exceeds some critical value Ri

*, turbulence will be suppressed. In laboratory experiments (Thorpe 1973) and DNS (Smyth and Moum 2000), mixing exhibits a hysteresis effect. As Ri decreases from large values (resulting from accelerating shear, e.g.), instability begins as Ri drops below ¼ (Miles 1961; Howard 1961). Subsequent mixing acts to increase Ri so that, absent additional forcing, Ri eventually becomes large enough that turbulence is extinguished. This occurs at a value of Ri roughly between 0.3 and 0.5, suggesting that the best value of Ri*

_{c}*will be greater than ¼. Some models have used much larger values (e.g., Ri*

_{c}*= 0.70; Large et al. 1994), but this is in compensation for low vertical resolution, which is not an issue here. We choose the value Ri*

_{c}*= 0.32 based on the laboratory experiments of Thorpe (1973) and our own DNS (Smyth and Moum 2000). We will see that our results show very little sensitivity to the value of Ri*

_{c}*.*

_{c}For 0 < Ri < Ri* _{c}*, we expect a smooth transition between strong mixing at Ri = 0 and no mixing at Ri = Ri

*. This transition has been modeled in various ways; for example, Large et al. (1994) chose diffusivity to be a cubic function of Ri. Here, we find that a linear function works as well as anything, and therefore choose it for simplicity (cf. Ellison and Turner 1959; Xu et al. 2006).*

_{c}A final issue, and a difficult one, is mixing driven by shear when Ri < 0. In this case, convection provides an alternative mixing mechanism that interacts with shear-driven turbulence in complex and poorly understood ways. For the present study, we make the conservative assumption that unstable stratification neither impedes nor amplifies shear-driven turbulence when Ri < 0 (although convection is assumed to drive mixing of its own accord, as described in the following subsection). This issue merits significant further study, because it transpires (section 3b) that much of the shear-driven mixing in thermohaline interleaving coincides with unstable stratification.

Taken together, these considerations lead to a mixing length parameterization with stratification effects represented via Ri,

with default parameter values *a*_{0} = 0.01 and Ri* _{c}* = 0.32. Here

*H*is defined as the vertical span over which Ri < Ri

_{s}*. The same diffusivity is applied to velocity and thermal and saline buoyancies, that is, the turbulent Prandtl and Schmidt numbers are taken to be unity.*

_{c}#### 3) Convective turbulence

The static instability is determined by a Rayleigh number,

where *H _{c}* is the convective layer thickness. In this study, the step function parameterization used in WR98 and Merryfield (2000) is replaced with a parameterization of the Nusselt number, Nu =

*K*

^{cnv}/

*κ*,

_{T}where Nu_{0} is a constant. The Rayleigh number dependence is established on dimensional grounds (e.g., Linden 2000). A default value of Nu_{0} = 0.09 is chosen by using a free boundary form of the empirical relationship described in Linden (2000). The large uncertainty in this value of Nu_{0}, as well as in the appropriate form for the Rayleigh number, motivates sensitivity tests to be described in section 4b. As with shear-driven turbulence, convective diffusivities of velocity and thermal and saline buoyancies are chosen to be equal.

### c. Environmental parameters

Because one aim of this study is to test mixing parameterizations, we sought parameter values from an extensively studied interleaving environment. Meddy Sharon offers such an environment (Armi et al. 1988). The data are pooled from four surveys over 2 yr that include velocity, temperature, and salinity profiles. The Meddy Sharon salt lens was located at 1000-m depth and was 2.5°C warmer and 0.6 psu more salty than ambient water (Hebert 1988), providing strong lateral temperature and salinity gradients. The upper flank of the lens was unstable to diffusive convection and the lower flank was unstable to salt fingering. Parameter values chosen to approximate the lower flank of Meddy Sharon include the squared Brunt–Väisälä frequency *N* ^{2} = (3 × 10^{−3})^{2} s^{−2}, the isohaline slope *ϕ* = 0.012, and the background buoyancy ratio *R*^{b}_{ρ} = 1.6 (Ruddick 1992).

The model run referred to as the base case uses the best estimates of model parameter values (*a*_{0} = 0.01, Ri* _{c}* = 0.32, Nu

_{0}= 0.09, S

*c*= 1.0, and Pr

_{f}*= 1.0) together with environmental parameter values reflecting the lower flank of Meddy Sharon, as described above. Subsequent runs are used to test the sensitivity of the results to variations in the model parameters.*

_{d}For these base-case parameter values, linear theory (Walsh and Ruddick 1995) yields an intrusion tilt angle of *θ* = 0.073°, a cross-intrusion wavelength of *λ _{ξ}* = 21 m, and a growth rate of 2.5 × 10

^{−6}s

^{−1}. The tilt angle and wavelength are built into the initial conditions; the growth rate is then used to test the accuracy of the numerics.

### d. The kinetic energy budget

A kinetic energy budget was calculated for analysis of model results. The derivation begins with the vertically integrated form of (11):

where *U* = ∂Ψ/∂*ξ* is the along-intrusion velocity. The terms on the right-hand side represent, first, the buoyancy forces and, second, viscosity resulting from parameterized microscale mixing processes. Because perturbations do not vary in *χ*, there are no advection or pressure gradient terms.

Multiplying (24) by *U* and integrating in *ξ* yields the evolution equation for kinetic energy:

The time rate of change of the total kinetic energy, KE = ∫(1/2)*U*^{2 }*d**ξ*, is balanced by the net release of potential energy resulting from buoyancy fluxes, *S*∫*Ub **d**ξ*, and the total energy dissipation resulting from eddy viscosity, *C*^{2}∫*A**ω*^{2 }*d**ξ*. We refer to these terms as KE* _{t}*, KE

*, and KE*

_{b}*, respectively. Because*

_{d}*A*is the sum of viscosities resulting from six different mixing processes (section 2b), there are six components of KE

*. These components will be separated to investigate the role of the different mixing processes. As with diffusivities, superscripts identify contributions from salt fingering KE*

_{d}^{sf}

_{d}, diffusive convection KE

^{dc}

_{d}, shear-driven turbulence KE

^{shr}

_{d}, convective overturning KE

^{cnv}

_{d}, ambient turbulence KE

^{at}

_{d}, and molecular viscosity KE

^{nu}

_{d}.

## 3. Base-case results

Model results are presented for the base case described in section 2c. To begin with, the kinetic energy growth rate and budget are used to classify and describe different stages of growth. This understanding is used as a basis to describe the effects of shear and convective instabilities.

Model intrusions exhibit the following four stages of growth: exponential growth (the linear regime), saturation, braking, and equilibration (Fig. 1). These stages can be identified by changes in the kinetic energy growth rate;

In the exponential growth stage, the growth rate of the kinetic energy is constant and closely resembles the linear theory prediction of 5.0 × 10^{−6} s^{−1} (Fig. 1b). During the saturation stage, KE* _{t}* decreases because the dissipation term grows more rapidly than the buoyancy term. Around 40 days, dissipation exceeds the buoyant forcing, the growth rate becomes negative, and the intrusions begin to decay. During this braking stage, KE

*and KE*

_{b}*fluctuate until they reach the equilibrium stage, at which point the change in kinetic energy approaches zero.*

_{d}The four stages of intrusion growth are explained, in part, by flux convergence and divergence, shown in Fig. 2. This provides a framework in which to discuss the individual and combined roles of the various microscale mixing processes. This understanding can then help explain how changes in the shear and convective parameterization affect interleaving (section 4).

Molecular and ambient diffusivities are independent of interleaving, and are therefore time invariant and present uniformly throughout the domain. When referring to a particular mixing region, we use the dominant instability type (aside from molecular and ambient); for example, a “salt fingering region” implies molecular and ambient turbulence diffusivities as well. Regions where only molecular and ambient turbulence are active are referred to as ambient turbulence.

### a. Exponential growth (days 0–28)

At this stage in flow development, salt fingering fluxes dominate the entire profile. Alternating layers of weaker and stronger salt fingering accelerate flow through buoyancy convergences and divergences. The strengthening and weakening of the salt fingering instability in different layers is thereby reinforced, resulting in a positive feedback loop and exponential growth.

### b. Saturation (days 28–42)

As the intrusions grow to finite amplitude, the positive feedback that sustains exponential growth at smaller amplitude is disrupted. First, the salt fingering fluxes no longer depend linearly on *R _{ρ}* as they do at smaller amplitudes. In addition, the following three new mixing processes become active: shear instability, convective overturning, and diffusive convection.

In the early saturation stage, a brief increase in the growth rate occurs when the ambient turbulence layer develops and before diffusive convection begins (the mixing regions are shown in Fig. 2, days 31.5 and 34.3; the increase in *σ*_{KE} is shown in Fig. 1b). This is followed by a much longer interval of decreasing growth rate.

Layers where initial *R _{ρ}* is small (Fig. 3, dashed curve) exhibit increasingly strong salt fingering as

*R*decreases to its asymptotic value near 1.3. Layers where initial

_{ρ}*R*is large exhibit considerably more complex evolution (Fig. 3, solid curve). It is here that diffusive convection, shear instability, and convective overturning develop. In this layer,

_{ρ}*b*

^{S}

_{ξ}and

*b*

^{T}

_{ξ}have signs opposite to those of the corresponding background gradients. As well,

*b*

^{S}

_{ξ}becomes larger in magnitude and changes more quickly than

*b*

^{T}

_{ξ}, causing

*R*to grow toward infinity before overcoming the background saline gradient and forcing

_{ρ}*R*to become negative (14). The net result is a layer where only ambient turbulence and molecular fluxes are active (Fig. 2,

_{ρ}*t*= 31.5 days, gray region).

The ambient turbulence layer is much less effective at smoothing gradients than the salt fingering layer, allowing gradients in velocity and buoyancy to sharpen. At the same time, *b*^{T}_{ξ} continues to grow. It eventually overcomes *B*^{T}_{ξ}, causing the net thermal buoyancy to switch sign. As a result, the flow becomes unstable in temperature and a diffusive convection layer develops (Fig. 2, day 34.3, yellow region, and Fig. 3).

The stratification now supports both salt fingering and diffusive convection layers, with thin layers of ambient turbulence separating the two (Fig. 2, between *t* = 31.5 and *t* = 34.3 days). In the ambient turbulence layer, thermal and saline diffusivities are equal, so the net buoyancy convergences and divergences are established by the neighboring salt fingering and diffusive convection layers. Velocity flux convergence and divergences lead to an increase in vorticity (Fig. 4, *t* = 34 days). Initially, this increase in vorticity leads to shear instability (Fig. 2, *t* = 35.8 days, light blue regions, and Fig. 4); however, the velocity flux convergences that lead to this increase in vorticity are accompanied by buoyancy fluxes that lead to a reduced buoyancy gradient (Fig. 2, *t* = 35.8 days). As the total buoyancy gradient becomes negative, the layer of shear instability becomes a layer of combined shear and convective instability (Figs. 4 and 2, *t* = 37.6 days, red regions). As the region of shear and convective instabilities grows, the kinetic energy growth rate decreases (Fig. 1b). This result is in contrast to Simeonov and Stern (2004), who employed a two-dimensional model of intrusions on a finite-width front and found that exponential growth continued after the development of overturns.

Throughout the flow development, thin layers of “pure” shear instability (i.e., in the absence of convective instability) persist along the borders of sheared layers where the buoyancy gradient is statically unstable. Shear instability is not strong enough, however, to prevent the development of the convective layers.

### c. Braking (days 42–60)

The braking stage is marked by a negative kinetic energy growth rate, beginning near day 42. At this stage, the spatial relationship between the different mixing processes is well established, although the layers continue to evolve in extent and strength (Fig. 2, *t* = 41.6 days).

Differential diffusion by salt fingering causes a buoyancy flux convergence (divergence) at the upper (lower) edge of each salt fingering layer. As a result, convective instability is maintained in the adjacent layers. The same buoyancy flux variations drive the along-intrusion velocities that characterize interleaving. Any variation in the strength of salt fingering, diffusive convection, or shear/convective mixing affects the relative strength of the convergent and divergent regions; hence, throughout the braking stage, the flow oscillates between accelerating and decelerating states.

### d. Equilibrium (days 60–end)

The evolution to the equilibrium state depends on the details of the initial condition, which is unlikely to be reproduced exactly in the ocean. Nonetheless, it is reasonable to assume that the equilibrium state represents the end result of a class of initial conditions lying within some basin of attraction in parameter space. The equilibrium state is therefore the most relevant to observable interleaving structures.

Equilibrium is attained when KE* _{b}* is matched by KE

*(Fig. 1c), causing*

_{d}*σ*

_{KE}to approach zero. The velocity profile has evolved from its initially sinusoidal form (Fig. 5a). The divergence of the buoyancy flux, which drives interleaving motions downward and to the left in the layer 22 m <

*z*< 32 m, is strongest at the lower edge of the salt fingering layer (

*z*≈ 29 m in Fig. 2, day 68, lower frames). Similarly, buoyancy flux convergence is strongest at the upper edge of the fingering layer (

*z*≈ 12 m in Fig. 2). These asymmetries account for the distortion of the velocity profile away from its initial sinusoidal form. The equilibrium buoyancy profiles are nonsinusoidal for similar reasons (Fig. 5b).

The velocity and buoyancy profiles are nearly linear in the shear/convective and diffusive convection layers, and more curved in the fingering layers. The maximum thermal and saline buoyancies are 2.8 and 2.4 × 10^{−4} m s^{−2}, respectively. The thermal buoyancy, having the larger amplitude, determines the net buoyancy perturbation, so that lighter layers flow to the right (uphill) and vice versa (Figs. 5a,b).

Shear/convective layers have the largest vertical diffusivities (Fig. 5c), but exhibit relatively weak kinetic energy dissipation (Fig. 5d) resulting from relatively weak shear (Fig. 5a). The strongest dissipation is found in fingering layers, where stable stratification allows strong shear to persist.

The flux characteristics of the equilibrium state are very similar to those found in the late saturation and braking stages (Fig. 2, *t* = 41.6 and 68 days). By equilibrium, however, the shear/convective mixing layer is able to balance the destabilizing effect of the interaction between the salt fingering and diffusive convection layers. This balance is reached by strengthening the buoyancy fluxes in the diffusive convection and shear/convection layers, an increase in the momentum flux in the shear/convection layer, and an increase in the thickness of the shear/convective layer (Fig. 2, *t* = 37.6–68 days).

The relative importance of the various microscale mixing processes in the equilibrium state may be assessed in terms of their respective contributions to the kinetic energy dissipation, which balances buoyancy forcing to maintain equilibrium (Figs. 1c,d). Salt fingering is the largest contributor to KE* _{d}* and accounts for 36% of the total. Ambient turbulence is not strong, but it is active throughout the domain and therefore supplies the second-largest fraction of dissipation, 22%. Dissipation resulting from convection accounts for 15% of the total, while shear instability and diffusive convection each account for 12% and molecular fluxes provide the remaining 3%.

Several properties of the modeled equilibrium state can be compared with results from other studies. These comparisons are summarized in Table 1. The maximum velocity is 4.6 × 10^{−3} m s^{−1}. Combined with the wavelength *λ _{ξ}* = 23 m and

*N*= 3 × 10

^{−3}s

^{−1}, this gives a Froude number of

*U*

_{max}/

*Nλ*= 0.067. For comparison, the Meddy Sharon intrusions had a nose velocity of 1 × 10

_{ξ}^{−3}m s

^{−1}(Ruddick and Hebert 1988). The interior velocity of interleaving layers (which we model here) has been estimated as 3.7 times the nose velocity (Ruddick et al. 1999). Using the same

*N*and

*λ*= 25 m (Ruddick 1992), one therefore obtains

_{ξ}*U*

_{max}/

*Nλ*= 0.049. The modeled and measured Froude numbers agree to within ∼25%. From laboratory results, Ruddick et al. (1999) estimated a Froude number of 0.005 based on the nose velocity, or 0.019 based on interior velocity. Both modeled and measured Froude numbers exceed this value significantly. A still-higher value of 0.14 has been obtained in the two-dimensional modeling experiments of Simeonov and Stern (2007).

_{ξ}The thermal and saline buoyancy amplitudes convert into temperature and salinity amplitudes of 0.19 K and 0.033 psu. The latter is somewhat smaller than the salinity amplitude of 0.05 psu estimated for Meddy Sharon (Ruddick and Hebert 1988).

The horizontal saline diffusivity was calculated from along-intrusion advective fluxes using

with 2*λ _{ξ}* representing the domain height. The cosine factor

*C*delivers the horizontal component of the along-intrusion flux. The thermal diffusivity is essentially equal to

*K*

_{Sx}. For the base case at equilibrium, the horizontal saline diffusivity is 3.7 m

^{2}s

^{−1}. Hebert et al. (1990) give estimates for Meddy Sharon ranging from 0.5 to 5 m

^{2}s

^{−1}. Our model result is within that range. Comparison of horizontal diffusivities is discussed in more detail in section 4d.

The kinetic energy dissipation rate KE* _{d}* has a domain-averaged value of 5.9 × 10

^{−11}m

^{2}s

^{−3}and a maximum value in the fingering layers of 2.4 × 10

^{−10}m

^{2}s

^{−3}. These values are comparable with the mean dissipation rate 2 × 10

^{−10}m

^{2}s

^{−3}inferred by Marmorino (1990) for turbulence in a thermohaline staircase, but are considerably smaller than the estimates made by Oakey (1988) for intrusions in Meddy Sharon. Oakey computed dissipation rates averaged over 50-m depth, which is similar to our vertical domain height 2

*λ*= 46 m. On the lower flank of the meddy (see Oakey’s Fig. 3, 900–1200-m depth), dissipation rates were in the range from 0.2 × 10

_{ξ}^{−9}to 2 × 10

^{−9}m

^{2}s

^{−3}. A potential explanation for this large discrepancy will be discussed in section 5.

The mean dissipation rate, combined with *N* = 3 × 10^{−3} s^{−1}, gives a buoyancy Reynolds number KE* _{d}*/

*νN*

^{2}= 6.6. This value agrees very closely with the value 5.9 found in the two-dimensional simulations in a similar parameter regime Simeonov and Stern (2007). The Cox number

*K*/

_{T}*κ*, using the mean thermal diffusivity (Fig. 5c), is 1070. Cox number measurements for Meddy Sharon are not available; instead, we cite for comparison the smaller values of 180, from the simulations of Simeonov and Stern (2007), and 256, measured in a thermohaline staircase by Gregg and Sanford (1987).

_{T}## 4. Sensitivity to parameterized mixing processes

Model parameters appearing in the microscale flux parameterizations were varied in order to discern how the strength of these processes affect interleaving.

### a. Shear-driven turbulence

Increasing the shear instability parameter *a*_{0} decreases the overall kinetic energy of the system (Fig. 6). Most of the difference between cases occurs during the saturation phase, revealing the influence of *a*_{0} on the development of shear instability and its secondary effect on convection and diffusive convection. The model is insensitive to variations in Ri* _{c}*, as can be seen from the fact that the curves representing Ri

*= 0.32 and Ri*

_{c}*= 1 are indistinguishable. This insensitivity is due to the small spatial and temporal extent of layers with 0 < Ri < Ri*

_{c}*.*

_{c}Based on (22), one might expect the net dissipation resulting from shear-driven turbulence to be proportional to *a*_{0}. Evidently, however, changes in *a*_{0} lead to changes in other parameters that govern the effect of shear instability, such as the magnitude of the shear in the unstable layers. When *a*_{0} is set to zero (i.e., neglecting shear-driven turbulence, as has been done in previous studies), *K*^{shr}_{d} vanishes; this effect is offset by a dramatic increase in the net dissipation resulting from convection (Fig. 7, thin, solid curve) and a more moderate increase in dissipation resulting from ambient turbulence and salt fingering. Because of these compensating effects, intrusion strength is less sensitive to *a*_{0} than it would be otherwise.

Interestingly, increasing *a*_{0} above its default value does not significantly affect shear-driven dissipation; instead, the difference is more striking in convection, salt fingering, and diffusive convection (Fig. 7). When *a*_{0} is doubled to 0.02, ∫KE^{shr}_{d }*dt* increases by a mere 10% while ∫KE^{cnv}_{d }*dt* decreases by 38%. Concurrently, ∫KE^{sf}_{d }*dt* decreases by 9% and ∫KE^{dc}_{d }*dt* increases by 6%. These changes indicate that salt fingering and convective instabilities become weaker while the diffusive convection instability becomes stronger, and vice versa. The overall effect (including the changes in ambient and molecular contributions) is a *decrease* of 8% in the net dissipation rate ∫KE_{d }*dt*. The increase in diffusivities resulting from shear-driven turbulence decreases the buoyancy gradients in the salt fingering and shear/convective instability region while increasing the gradients in the diffusive convection region.

Increasing *a*_{0} by an order of magnitude has a similar effect on ∫KE^{cnv}_{d }*dt*, ∫KE^{sf}_{d }*dt*, and ∫KE^{dc}_{d }*dt* as increasing *a*_{0} by a factor of 2: ∫KE^{cnv}_{d }*dt* decreases, ∫KE^{sf}_{d }*dt* decreases, and ∫KE^{dc}_{d }*dt* increases, resulting in less vertical flux in the salt fingering and shear/convective instability layers, with more vertical flux in the diffusive convection layer. Unlike *a*_{0} = 0.02, however, *a*_{0} = 0.1 exhibits less ∫KE^{shr}_{d }*dt* than the base case. This decrease is caused by a decrease in the vorticity within the shear/convective layer (Fig. 2, *t* = 68 days). Increasing *a*_{0} results in a large velocity flux in the shear/convective region, leading to stronger convergences and divergences at the boundaries with the salt fingering and diffusive convection layers. This combined effect acts to decrease vorticity in the shear/convective layer.

Overall, increasing *a*_{0} has the effect of increasing dissipation and mixing in the shear/convective instability layer. This increase leads to a weakening of the salt fingering instability and a strengthening of the diffusive convection instability. At a certain point, the weakening and strengthening of these instabilities moderates the gradient within the shear/convective instability layer to the point that, in spite of a larger *a*_{0}, dissipation resulting from shear instability is less than it is with *a*_{0} = 0.01 or 0.02. A decrease in the buoyancy gradient in the shear and convective layer is reflected in the decreasing KE^{cnv}_{d} that accompanies an increase in *a*_{0}, as shown in Fig. 7.

### b. Convection

An early hypothesis was that, as soon as interleaving amplified to the point that Ri dropped below Ri* _{c}*, shear instability would act to oppose further sharpening of the buoyancy gradient, and thereby prevent the development of convectively unstable layers. This expectation proved false. With default parameter values, convection remained roughly equal in importance to shear-driven turbulence; and even with

*a*

_{0}set an order of magnitude greater than our best estimate, significant convective layers remained. In this section, we test for sensitivity to the parameterization of convective turbulence.

The effect of varying Nu_{0} is remarkably similar to that of varying *a*_{0} (Figs. 8, 9). The total kinetic energy in the system decreases with an increasing Nu_{0}. In the shear/convective instability layer, the absence of convective dissipation is compensated by dramatically increased shear instability (Fig. 9). Overall, in the absence of the convective instability, kinetic energy dissipation in salt fingering, shear instability, ambient turbulence, and molecular viscosity all increase, with shear instability exhibiting the largest change. The diffusive convection layer is an exception in that the kinetic energy dissipation decreases. Net dissipation *increases* by 23% when convection is eliminated from the model.

When Nu_{0} is doubled from its base value of 0.09 to 0.18, kinetic energy dissipation resulting from convection increases only slightly while that from shear instability is roughly halved. Beyond a value of Nu_{0} = 0.18, the energy dissipation resulting from convection decreases. Dissipation resulting from shear instability also decreases, as do all other terms except KE^{dc}_{d}, which shows a slight increase.

### c. The relationship between shear instability and convective overturning

The close relationship between the role of convective overturning and that of shear-driven turbulence in this model warrants further discussion. In the layers where salt fingering is weak initially, buoyancy fluxes drive the fluid toward a statically unstable state (section 3b). The result is the unstable layers surrounding the diffusive convection layer visible in Fig. 2. Other mixing processes are unable (in the present parameter regime, at least) to prevent these layers from developing.

In the statically unstable layers, double-diffusive mixing is inactive. The remaining mixing processes (shear and/or convection, depending on parameter choices, as well as ambient and molecular mixing) all differ from double diffusion in that they drive purely downgradient fluxes, that is, they always act to smooth gradients. As a result, these processes participate in a negative feedback mechanism that compensates for changes in the strength of any single process. Should any process be weakened or removed, the proximate effect is a sharpening of gradients, which leads to an increase in fluxes driven by the remaining processes. Conversely, an increase in the diffusivity of any process acts to smooth gradients and thus to reduce fluxes. In either case, the effect of a change in any diffusivity, for example, a change in *a*_{0} or Nu_{0}, is automatically mitigated.

Because of this negative feedback, the addition of shear-driven mixing to the model has only a minor effect on the overall strength of interleaving. Instead, the salient effect of shear-driven mixing is to reduce the unstable buoyancy gradients in the convective layers. When shear instability is neglected (*a*_{0} = 0) and all other parameters are set to default values, the minimum vertical buoyancy gradient in the equilibrium state is −3.5 × 10^{−6} s^{−2}. Increasing *a*_{0} to its default value 0.01 changes this minimum buoyancy gradient to −2.7 × 10^{−6} s^{−2}.

### d. Momentum fluxes by diffusive instabilities

There is evidence that the effective viscosity of double-diffusive motions can differ dramatically from the scalar diffusivities. For example, the Prandtl number for diffusive convection has been estimated by Padman (1994) to lie between 1 and 3. Recent linear analyses by Smyth and Kimura (2006) agree with Padman’s values of Pr* _{d}* while suggesting a contrasting result for salt fingering; Sc

*is predicted to be of the order of 10*

_{f}^{−1}or even lower. In the laboratory experiments of Ruddick et al. (1989), diffusive convection yielded an effective viscosity twice that of salt fingering. To assess the significance of these parameter variations for interleaving, we performed four additional sensitivity tests. These were identical to the base case except that both Pr

*and Sc*

_{d}*, in turn, were increased to 3 then reduced to 0.3.*

_{f}Varying the friction because of salt fingering had a significant impact on the growth rate in the linear regime (Fig. 10, dashed curves), whereas varying Pr* _{d}* did not (dotted curves). This is to be expected, in this case, because linear growth is driven by salt fingering. Equilibration occurred at a higher energy level when either Sc

*or Pr*

_{f}*was decreased and at a lower energy when either parameter was increased.*

_{d}### e. Horizontal diffusivity

The horizontal saline diffusivity resulting from interleaving in Meddy Sharon has been estimated under various plausible assumptions, yielding values between 0.5 and 5 m^{2} s^{−1} (Hebert 2001). The uncertainty in this estimate is significant; nonetheless, the estimate is very precise when viewed in the context of interleaving diffusivities in the full range of oceanic parameter regimes, which range from zero to several thousand squared meters per second (Ruddick and Kerr 2003). Hebert et al.’s (1990) quantification of the Meddy Sharon diffusivities therefore provides a valuable opportunity to assess both the accuracy of our one-dimensional model and the impact of uncertainty in the microscale flux parameterizations discussed previously in this section.

Diffusivity evolution for the base case and all of the sensitivity tests are shown in Fig. 11. The observational bounds are represented by horizontal lines. Parameter variations lead to changes on the order of a few tens of a percent, but in every case the diffusivity is within the range of the observational estimates. This agreement offers strong support for Stern’s original insight that interleaving is driven by double-diffusive instability, for WR98’s choice to model interleaving as a one-dimensional process, and for the parameterization [(15)–(17)] of the salt fingering fluxes that are the primary drivers of interleaving in this regime. However, the observational estimates of diffusivity are not precise enough to discriminate among the parameter values used in our shear-driven and convective flux parameterizations.

## 5. Conclusions

We have introduced a shear instability parameterization and an improved convection parameterization to the interleaving model described in WR98. Shear instability and convection occur almost simultaneously. In sensitivity experiments, increases in shear-driven turbulence compensate for the absence of convective turbulence, and vice versa. When either instability is eliminated from the model entirely, the horizontal diffusivity is increased by approximately 30%. In other words, both shear and convective turbulence act to decrease the net amount of kinetic energy in thermohaline intrusions and lead to a significant decrease in horizontal fluxes. A criticism of previous 1D interleaving models is that they require unrealistically strong density inversions to equilibrate (Simeonov and Stern 2007). By improving the shear and convective parameterizations, we have found that a 1D model can equilibrate with reduced density inversions. We emphasize, however, that convective layers persist even when the amplitude constant *a*_{0} that governs shear-driven mixing is an order of magnitude above our best estimate.

Viscosity resulting from double-diffusive processes has proven important in determining both initial growth rates and final energy levels. We have seen this by varying both the salt fingering Schmidt number and the Prandtl number for diffusive convection. In this fingering-favorable environment, the former parameter controls the initial growth rate, while both parameters affect the equilibrium state.

A comparison with observations was accomplished mainly via predictions of the horizontal saline diffusivity resulting from interleaving in Meddy Sharon (Hebert et al. 1990). The value for the base case in the equilibrium state is 3.7 m^{2} s^{−1}. Uncertainty in model parameter values leads to significant uncertainty in the model result, but in every case the predicted diffusivity is within the error of the observations. Good agreement was also found between modeled and measured intrusion velocities, Froude numbers, and salinity amplitudes. This correspondence provides strong support for the essential understanding of interleaving that the model represents.

The kinetic energy dissipation rate is smaller than the observed value by at least an order of magnitude. The buoyancy Reynolds number of the base-case equilibrium state of 6.6 is also small but agrees very well with the two-dimensional modeling result of Simeonov and Stern (2007). This suggests that the low dissipation rate may be due to some mechanism not considered either in the present model or in Simeonov and Stern (2007). One possibility is the effects of a sheared alongfront current.

An alongfront current may be sheared either vertically (to maintain thermal wind balance with a baroclinic density gradient), in the cross-front direction, or both. Such a current is excluded in the barotropic model used for the present paper. In the baroclinic case, as we envision it, intrusions grow from the linear instability described by May and Kelley (1997). Alongfront momentum (or more properly for the meddy case, azimuthal momentum) is advected by the intrusions, creating vertical shear of alongfront velocity at intrusion scales. The resulting shear reduces Ri and contributes to turbulent mixing and energy dissipation. In the meddy case, the cross-front momentum flux would cause the meddy to spin down and lose kinetic energy. Simultaneously, geostrophic adjustment would convert potential to kinetic energy, which would in turn be dissipated by the same mechanism until the energy of the meddy was exhausted. This additional mixing mechanism could account for the order-of-magnitude discrepancy between the kinetic energy dissipation rate in the present model and the values measured in Meddy Sharon. These effects are readily added to the 1D model described here, and will be the subject of a future publication.

While this study offers encouraging results, it also reveals opportunities to improve our understanding. The assumption that all the mixing processes are independent is challenged by current research. For example, Kimura and Smyth (2007, manuscript submitted to *Geophys. Res. Lett.*) show that salt fingering fluxes are affected by ambient shear, and that interaction is not accounted for in the present model. The interaction between shear and convective mixing must be better understood. Model results should be compared with observations other than Meddy Sharon. For this, allowance must be made for baroclinic effects and for background environments that are unstable to diffusive convection rather than salt fingering. This study also reveals a need for more detailed observations of horizontal diffusivity resulting from interleaving in all oceanic regimes.

## Acknowledgments

This paper is based on the first author’s M.S. thesis. Valuable input was provided by advisory committee members Jim Moum, Yvette Spitz, and Jim Liburdy. We thank Julian Simeonov for helpful input and for sharing results prior to publication. This project was supported by the National Science Foundation under Grants OCE-0221057 and OCE-0622922.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* William D. Smyth, 104 COAS Admin Bldg., Oregon State University, Corvallis, OR, 97331. Email: smyth@coas.oregonstate.edu