# An introduction to mixing in a stably stratified fluid

Chantal Staquet

► To cite this version:

Chantal Staquet. An introduction to mixing in a stably stratified fluid. Statistical theories and computational approaches to turbulence, Springer Verlag, pp.60-79, 2003, 9784431670049. 10.1007/978-4-431-67002-5_4. hal-00264750

HAL Id: hal-00264750

https://hal.science/hal-00264750v1

Submitted on 8 Jan 2020

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An Introduction to Mixing in a Stably Stratified Fluid

Chantal Staquet

Laboratoire des Ecoulements Géophysiques et Industriels, BP 53, 38041 Grenoble cdx 9, France

**Abstract.** We provide an account of mixing in stably stratified fluids: we present the basic concepts of mixing, recall how mixing is quantified and analyze mixing in the stably stratified shear layer along these lines. We also discuss the importance of mixing in the ocean and describe how it is modeled in this context, focusing upon the commonest approaches.

## 1 Introduction

General circulation models of the ocean, which provide a numerical description of the mean currents at oceanic basin scales, use resolution of a few km along the horizontal and of tens-to-hundreds of meters along the vertical. It is therefore necessary to model the motions below the grid scale, the so-called subgrid-scale motions, at which kinetic energy is dissipated and mixing occurs. This modeling is crucial because the dissipating processes contribute fundamentally to the large scale dynamics of the ocean. Let us illustrate this important point through two arguments, related to the thermohaline circulation.

The thermohaline circulation is a large scale convective motion formed by the sinking at high latitude of cold water masses; these cold water masses spread throughout the bottom of ocean basins and slowly upwell. This upwelling is due to mixing processes: without deep mixing, the ocean would have turned, within a few thousand years, into a stagnant pool of cold salty water extending to within one meter from the surface. The thermocline, which separate the bottom waters from the surface waters mixed by the wind, would have reduced to a very thin layer. This yields the fundamental idea that the thermohaline circulation is not driven by buoyancy effects but by the mechanical energy that mix the water masses, namely wind and tidal energy (see [15] and [24] and references therein).

A global estimate of the turbulent diffusivity associated with mixing in the ocean was proposed by Munk [23], by assuming a balance between the upwelling water masses and the downgradient mixing of heat. Using the density profile measured in the ocean below the thermocline and an estimate of the vertical velocity of upwelling from the distribution of chemical tracers, Munk [23] obtained a value for the turbulent diffusivity equal to $\kappa_t \simeq 10^{-4} m^2/s$ (that is, one thousand times the molecular diffusivity). The method used by Munk was subjected to many discussions (and some controversy) but $\kappa_t$ is now considered as the canonical value for the global diffusivity in the ocean.
As shown by Bryan [3] already in 1987, a constant diffusivity is certainly
too crude an approximation for the parameterization of mixing in general cir-
culation models. A refined parameterization is needed, namely a simple relation
between a turbulent diffusivity and a dynamical parameter characterizing the
large resolved scales of the model (the same need arises for the turbulent viscos-
ity). This is still an open problem: for instance, Beckmann [2] provides results of
three numerical simulations of the thermohaline circulation using general circu-
lation models differing by the representation of the bottom topography and the
modeling of the subgrid scales. Currents for depths greater than 500 m are plot-
ted in Fig. 1: strong discrepancies between these large scale currents are visible
from one model to the other, which illustrates their strong sensitivity to small
scale processes.

It is necessary at this stage to provide a definition of mixing. Mixing is the
destruction by molecular diffusion of density gradients generated by turbulence.
Some authors (such as Toole [37]) also include in this definition of mixing the
destruction of velocity gradients by viscosity. We shall rather keep the usual
terminology of kinetic energy dissipation for the latter process. In the ocean,
mixing takes place at scales of 1 mm to 1 m.

Mixing may also be referred to as vertical transport: because of the existence
of the restoring force of buoyancy along the vertical direction, no net vertical
transport can occur in a stably stratified fluid without mixing. By contrast, no
restoring force exists along the horizontal direction (ignoring the Coriolis force),
so that horizontal motions are not constrained. Horizontal transport thus leads
to characteristic length scales that are much larger than those associated with
vertical transport, because the former process is controlled by advective effects
while the latter results from dissipating processes. Finally, we point out that
mixing should be made distinct from stirring. Stirring results from the straining
by large scale motions, which increases density gradients without mixing. An
example of stirring is provided in Fig. 2.

This introduction mainly focused upon oceanic processes so that the active
scalar field, which is mixed, is the density. The concept of mixing also applies in
the atmosphere, by considering the potential temperature field in place of the
density field. However, all concepts presented in the present paper rely upon
the Boussinesq approximation so that they must be modified when compressible
effects become important (see [1]).

The outline of the paper is the following. In the next section we present the
basic concepts of mixing along with their modeling in the ocean. We provide
an example of mixing analysis in section 3, by considering three-dimensional
numerical simulations of the stably stratified shear layer (most results presented
in this section stem from [32] and [33]. We also briefly address the question of
parameterization of mixing resulting from shear instability, both in the numerical
simulations and in large scale oceanic models. Conclusions are drawn in the final
section.
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Time-mean bottom layer circulation in the subpolar North Atlantic from the 1/3 degree DYNAMO-SPEM general circulation model. (b) Same as (a), for the DYNAMO-MICOM general circulation model. (c) Same as (a), for the DYNAMO-MOM general circulation model (from [2]).</figcaption></figure>

## 2 Basic concepts of mixing

### 2.1 Background and available potential energy

Energy is required to mix a stably stratified fluid. Even if we consider a fluid at rest, made of two layers of different density, the denser one being at the bottom, the interface mixes through molecular effects : the density gradient smooths, the source of energy being the internal energy of the fluid. The diffusion of
<figure><img src="image_3.png" /></figure>

**Fig. 2.** The large scale vortex that develops from the Kelvin-Helmholtz instability induces a strain field which results in a thin layer of strong density gradient outside the vortex. The wrapping of dense and light fluid as the vortex forms creates intertwined thin layers of different density, thus yielding strong density gradients within the vortex. The large scale motion associated with the vortex mainly stirs the fluid and hardly mixes it (from [33]).

density at the interface occurs through a diapycnal flux of mass ( i.e. across the constant density surfaces, or isopycnals). This is expressed by the standard diffusion equation:

$$ \frac{\partial \rho_b}{\partial t} + \frac{\partial}{\partial z} \phi_{lam} = 0. \qquad (1) $$

$\rho_b(z)$ is the density profile of the fluid and

$$ \phi_{lam}(z) = -\kappa \frac{d\rho_b}{dz} \qquad (2) $$

is the local diapycnal mass flux. This flux results in an increase of the center of mass of the fluid so that the potential energy of the fluid $E_b$ (usually defined per unit volume) increases as well. The expression for $E_b$ is

$$ E_b = \frac{1}{V} \int_V \rho_b(z) gz dV, \qquad (3) $$

where $V$ is the volume of the fluid, assumed in the following to be contained in a fixed parallelepipedic rigid box, bounded by a surface $S$. The evolution equation for $E_b$ is easily inferred from Eq. (1):

$$ \frac{d}{dt} E_b = \Phi_{lam}, \qquad (4) $$

with

$$ \Phi_{lam} = \frac{\kappa g}{H} \Delta \rho = g \frac{1}{V} \int_V \phi_{lam} dV > 0. \qquad (5) $$
$\kappa$ is the molecular diffusivity for density changes, $g$ is the acceleration of gravity, $H$ is the height of the fluid volume and $\Delta\rho = \rho_b(0) - \rho_b(H)$ is the (positive) density difference between the bottom and the top of the container. $\Phi_{lam}$ is the rate of production of the potential energy of the fluid at rest. Equivalently, $\Phi_{lam}$ is the conversion of internal energy into potential energy. We shall refer to this rate as the laminar mixing rate for simplicity.

In a moving fluid, mixing still occurs through a local diapycnal *diffusive* flux of mass, which we refer to as $\phi_d$. The *advective* flux of mass $g\rho u$ only distorts the isopycnals without changing their value. In order to access the diffusive mass flux $\phi_d$, let us do a thought experiment: we assume that the fluid relaxes instantaneously toward a rest state through an adiabatic redistribution of density. Each fluid element is thus brought back to its equilibrium position without changing its density. This transformation amounts to stretching out each isopycnal horizontally; it yields a stable density profile (denoted $\rho_b$ again), whose potential energy ($E_b$) is the minimum potential energy of the fluid at that time. These concepts have been introduced by Lorenz [22]. Hence, the above situation of a fluid at rest is recovered instantaneously through the adiabatic transformation so that the same equations should hold formally: the changes in $\rho_b$ are only due to the diapycnal mass flux $\phi_d$

$$ \frac{\partial \rho_b}{\partial t} + \frac{\partial}{\partial z} \phi_d = 0. \qquad (6) $$

This result has been demonstrated rigorously by Winters & D'Asaro [39]. The expression for $\phi_d$ is:

$$ \phi_d(z) = \kappa \frac{\langle |\nabla \rho|^2 \rangle_I}{-d\rho_b/dz} > 0, \qquad (7) $$

where $\langle\rangle_I$ denotes an average along an isopycnal and $\rho$ refers to the total density field. The equation for the background potential energy is analogous to Eq. (4)

$$ \frac{dE_b}{dt} = \Phi_d + Fb_{surf}^{diff} + Fb_{surf}^{adv}, \qquad (8) $$

up to diffusive ($Fb_{surf}^{diff}$) and advective ($Fb_{surf}^{adv}$) mass fluxes across the bounding surface $S$ (exact expressions of the latter fluxes are provided in [38]). As before, $\Phi_d$ is the volume averaged value of $\phi_d$ times $g$:

$$ \Phi_d = g \frac{1}{V} \int_V \phi_d \, dV \qquad (9) $$

and is the rate of production of the background potential energy. Hereafter, we simply refer to $\Phi_d$ as the mixing rate. When the fluid is at rest, $|\nabla\rho| = -d\rho_b/dz$ and $\Phi_d$ reduces to $\Phi_{lam}$. In general $\Phi_d \ge \Phi_{lam}$ ($\ge 0$), expressing the well-known fact that turbulent motions (more precisely: density fluctuations) increase mixing. Note that density fluctuations refer here to deviations about the background density field $\rho_b$.

Lorenz' concepts were first applied by Thorpe [35] for the case of one-dimensional density profiles measured in a lake: the idea was to sort out the elementary
fluid volumes contributing to the fluid column by order of increasing density, so as to get an estimate of the vertical length of turbulence. Indeed, the sorting method only operates when the fluid has overturned, fluid volumes being left at their position when the fluid is locally stable. Thorpe recorded the distance each particle has to be displaced to reach its position in the stable profile and defined a length scale from the root mean square value of the distances. This length scale, now referred to as the Thorpe scale, is thus a measure of the vertical extent of overturning events in the original density profile. It also provides an estimate of the size of the turbulent eddies.

In a three-dimensional fluid volume, the use of an adiabatic transformation to analyse mixing has been performed only recently, starting with the work [38]; references are provided in the next section. The practical method to compute the stable density profile $\rho_b$ is described in this section as well; we also show that $\rho_b(z)$ is distinct in general from a mean (horizontally averaged) density profile.

The adiabatic transformation that leads to $\rho_b$ releases the energy associated with the density fluctuations. This energy is the available potential energy, introduced by Lorenz [22] also:

$$E_a = E_p - E_b, \quad (10)$$

where $E_p$ is the potential energy per unit volume of the fluid before the transformation is applied. (Closed expressions for $E_a$ are provided by Holliday & McIntyre [14]). The minimum potential energy state can also be characterized by $E_a = 0$.

The evolution equation for $E_a$ is simply derived from that for $E_b$ and $E_p$. We recall the evolution equation for the total potential energy per unit volume (see f.i. [38])

$$\frac{d}{dt}E_p = g \frac{1}{V} \int_V \rho w dV + F p_{surf}^{adv} + \Phi_{lam} + F p_{surf}^{diff}, \quad (11)$$

the first term on the right-hand-side is the volume averaged buoyancy flux, the second term is the advective mass flux through the boundaries, $\Phi_{lam}$ is the laminar mixing rate and the last term is the mass flux across the boundaries due to molecular effects. Using (8), (10) and (11), one gets the evolution equation for $E_a$:

$$\frac{d}{dt}E_a = \frac{1}{V} \int_V g \rho w dV - (\Phi_d - \Phi_{lam}) + (Fp_{surf}^{diff} - Fb_{surf}^{diff}) + (Fp_{surf}^{adv} - Fb_{surf}^{adv}). \quad (12)$$

This equation shows that the dissipation rate of $E_a$ within the fluid volume, $\Phi_d - \Phi_{lam}$, is the source for turbulent fluid mixing. The study of mixing thus amounts to decomposing the total potential energy $E_p$ into a reversible part $E_a$, available for mixing, and an irreversible part $E_b$, which increases through mixing.
## 2.2 Quantification of mixing: diapycnal diffusivity and mixing efficiency

The rate of change of the background potential energy $E_b$ per se does not provide any usable information on mixing. $dE_b/dt$ needs to be related to another quantity to provide some quantitative information. Equation (12) for $E_a$ shows that the mixing rate due to density fluctuations, $\Phi_d - \Phi_{lam}$, should rather be used in place of $\Phi_d$ when mixing resulting from the flow dynamics is to be studied. $\Phi_d - \Phi_{lam}$ is referred to as the turbulent mixing rate hereafter.

One common measure of mixing is provided by a flux Richardson number, $Rf_b$, which relates the rate of production of potential energy due to turbulent mixing within the fluid volume $\Phi_d - \Phi_{lam}$ to the rate of energy input to the fluid $\dot{M}$:

$$Rf_b|_{forced} = \frac{\Phi_d - \Phi_{lam}}{\dot{M}}. \qquad (13)$$

When the fluid is not forced, the total sink of energy within the fluid volume is used in place of $\dot{M}$. The energy is lost through mixing and through kinetic energy dissipation $\epsilon$ so that, in this case, the flux Richardson number within the fluid volume is expressed as:

$$Rf_b = \frac{\Phi_d - \Phi_{lam}}{(\Phi_d - \Phi_{lam}) + \epsilon}. \qquad (14)$$

The changes of energy from the beginning of the flow evolution may also be used in place of the instantaneous rates of change, leading to the definition of a global flux Richardson number:

$$Rf_b|_{global} = \frac{\Delta E_b}{\Delta E_b + \Delta E_k}. \qquad (15)$$

The interest in the latter definition is that it leads to a smoother evolution of the flux Richardson number and accounts for the boundary mass fluxes. $E_k$ is the volume averaged kinetic energy. Note that either definition of the flux Richardson yields a value comprised between 0 and 1. A very small value ($\Delta E_b \ll \Delta E_k$) implies that very weak mixing occurs while a value close to 1 ($\Delta E_k \ll \Delta E_b$) implies that mixing is very efficient. In the latter case indeed, much more energy is lost into mixing than into heat. In the literature, the important concept of mixing efficiency is rather defined by the quantity $\gamma_b = 1/(Rf_b^{-1} - 1)$, instead of $Rf_b$ itself.

Mixing is also most commonly characterized by a turbulent diffusivity $K_\rho$. When the scale of the background density gradient is much larger than the overturning scale (as measured by the Thorpe scale for instance), the turbulent flux is assumed to depend linearly on the gradient

$$\phi_d(z) = -K_b(z) \frac{\partial \rho_b}{\partial z}, \qquad (16)$$

the coefficient defining the turbulent diffusivity. The definition ensures that $K_b(z)$ is positive at any time and for any value of $z$. When the background
density gradient is linear, a volume averaged turbulent diffusivity can be defined

$$ \frac{1}{V} \int_V K_b(z) dV = \kappa \frac{\Phi_d}{\Phi_{lam}}. \qquad (17) $$

One important issue is the dependency of the turbulent diffusivity $K_b$ upon the molecular diffusivity $\kappa$. When fluid motions are strictly horizontal, $K_b = \kappa$. When fluid motions displace the isopycnals but remain laminar, as for internal gravity waves of low amplitude, $K_b$ is proportional to $\kappa$. When fluid motions are turbulent, a prediction may be conjectured when $\kappa \to 0$, the Prandtl number being set to a constant value, from analogy with ordinary turbulence. In the latter case, when $\nu \to 0$, the dissipation rate of kinetic energy remains finite and $\nu$-independant. The analogy suggests that, as $\kappa \to 0$, the mixing rate $\Phi_d$ remains finite and $\kappa$-independent, and so is the volume averaged value of $K_b$, from Eq. (17). When the flow is not laminar and $\kappa$ is not small, no theory exists, to our knowledge, and one should rely upon (numerical or laboratory) experiments to compute $K_b$.

## 2.3 Characterization of mixing in geophysical fluids: traditional approach

Mixing has been investigated for several decades in the ocean. As we discuss it in the present subsection, the major difficulty in such studies is the computation of the turbulent mixing rate ($\Phi_d - \Phi_{lam}$) so that models have always been used to estimate it. Hence, the exact computation of the mixing rate from the background density profile described above is a novel approach of mixing.

In oceanic models, the turbulent mixing rate is usually estimated from the average value of the advective flux of density:

$$ \Phi_d - \Phi_{lam} \approx g \overline{\rho' w'}, \qquad (18) $$

the average being a temporal, a spatial or an ensemble one; it is denoted by an overbar in the following (thus $\rho' = \rho - \bar{\rho}$, $w' = w - \bar{w}$, etc). The idea of taking the averaged advective density flux to estimate the turbulent mixing rate is that the average filters out the oscillations of the advective flux and yields the residual contribution of diffusive processes. Unlike in the atmosphere, the direct measurement of $\overline{\rho'w'}$ has been attempted only rarely in oceanography. The main reason is that, in addition to the smallness of this residual contribution, most instruments profile vertically in the ocean and therefore cannot measure vertical velocities (see, e.g., [11] and [37] for a review). Therefore the oceanic turbulent mixing rate is usually inferred from microstructure measurements (at scales of about 1 m) through statistical models. In this approach, the turbulent diffusivity is defined by

$$ K = \frac{\overline{\rho'w'}}{-d\overline{\rho}/dz} = \frac{\overline{\rho'w'}}{N^2}, \qquad (19) $$

where $N^2 = -(g/\rho_0)d\overline{\rho}/dz$ is the Brunt-Väisälä frequency.
A now classical statistical model was proposed by Osborn & Cox [26]. It relates the turbulent mixing rate to the dissipation rate of the variance of the density fluctuations. Assuming that the mean density gradient is about constant over the integral scale of the turbulence and that a statistically stationary regime has been reached, the equation for the variance of the density fluctuations reduces to:

$$ \overline{\rho'w'} \frac{d\overline{\rho}}{dz} = -\kappa |\overline{\nabla\rho'}|^2. \qquad (20) $$

This equation states that the production of fluctuating density variance by the action of the turbulent vertical flux of density $\overline{\rho'w'}$ on the mean gradient $\overline{d\rho/dz}$ is balanced by the dissipation rate of this variance $\kappa |\overline{\nabla\rho'}|^2$. Using definition (19) for the turbulent diffusivity, Eq. (20) yields the expression:

$$ K' = \kappa \frac{|\overline{\nabla \rho'}|^2}{(\overline{d\rho/dz})^2}. \qquad (21) $$

This relation was originally established by Osborn & Cox [26] for the variance of the temperature fluctuations. It can be expressed in terms of the density field, as done here, when the Boussinesq approximation is valid.

The same approach was used by Osborn [25] when the turbulence is sustained by a mean shear flow $\overline{U}(z)$. Under the same assumption of a statistically steady state, the equation for the velocity variance is (written per unit mass):

$$ \overline{u'w'} \frac{d\overline{U}}{dz} = -\epsilon - \frac{g}{\rho_0} \overline{\rho'w'}. \qquad (22) $$

As above, this equation states that the production of turbulent kinetic energy per unit mass by the action of the Reynolds stress $\overline{u'w'}$ onto the mean velocity gradient $\overline{d\overline{U}/dz}$ is balanced in the mean by the sinks of turbulent kinetic energy; these are the dissipation rate of kinetic energy $\epsilon$ and the turbulent mixing rate per unit mass $g/\rho_0 \overline{\rho'w'}$ (which increases the potential energy resulting from mixing). Introducing a flux Richardson number:

$$ Rf = \frac{\frac{g}{\rho_0} \overline{\rho'w'}}{-\overline{u'w'} \frac{d\overline{U}}{dz}}, \qquad (23) $$

Eq. (22) can be rewritten in a simplified way:

$$ \frac{g}{\rho_0} \overline{\rho'w'} = \frac{Rf}{1-Rf}\epsilon. \qquad (24) $$

Using the definition of the mixing efficiency $\gamma = \frac{Rf}{1-Rf}$, the turbulent mixing rate (per unit mass) becomes related to the dissipation rate of kinetic energy (instead of the dissipation rate of fluctuating density variance):

$$ \frac{g}{\rho_0} \overline{\rho'w'} = \gamma\epsilon. \qquad (25) $$
The turbulent diffusivity within this model is thus defined by, using (19) and (25)

$$K_{\rho}'' = \gamma \frac{\epsilon}{N^2} = \frac{R_f}{1 - R_f N^2} \quad (26)$$

From a practical point of view, $\epsilon$ is measured from a single component, assuming isotropy (the validity of this assumption has been addressed in [40] and [36]). A constant value of 0.2 for $\gamma$ is generally adopted. This value has been inferred from a conjecture by Osborn [25] (stemming from observations) that $R_f \le 0.15$. Because the maximum value for $R_f$ is not firmly established, an empirical value for $\gamma$ equal to 0.2 is used. As pointed out by Caldwell & Moum [4], $\gamma$ is not constant: values comprised between 0.05 and 0.48 have indeed been computed from microstructure measurements in the ocean (assuming $K_{\rho}'' = K_{\rho}'$), which seem to depend upon the instability that creates the turbulence. In the atmosphere, $\gamma$ displays the same variability about the 0.2 value.

# 3 Mixing in a stably stratified shear layer

In section 2, we introduced basic concepts of mixing along with two estimators of mixing (the flux Richardson number and the turbulent diapycnal diffusivity). We now illustrate these theoretical points and consider for this purpose the stably stratified shear layer. Most of the results presented in this section stem from [32] and [33].

## 3.1 The stably stratified shear layer: a brief introduction

The shear layer is the simplest case of a inhomogeneous shear flow. It is formed by two layers of fluids having different densities and horizontal velocities. The interface between the layers, which thickens by molecular diffusion, enlarges most efficiently if the interface is unstable. A necessary condition for instability is the Richardson number

$$J = \min_z \frac{g \frac{d\bar{\rho}}{dz}}{\rho_0 (d\bar{U}/dz)^2} \quad (27)$$

to have a value smaller than 0.25 somewhere in the fluid (see e.g. [10]). The overbar in Eq. (27) refers to a horizontal average. If this condition is met, the Kelvin-Helmholtz instability develops in the vertical plane of the shear flow. It is therefore a two-dimensional instability. As the amplitude of the perturbation grows, nonlinear effects gain importance and makes the instability saturate while the flow vorticity organizes into quasi-horizontal vortices (with axis parallel to the mean shear). The vortices are most unstable to a (two-dimensional again) subharmonic instability which makes them pair (see e.g. [8] and, for a review, [13]). Secondary three-dimensional instabilities may also occur, at smaller scales than the two-dimensional instabilities. Further detail about the three-dimensional instabilities in a stably stratified shear layer may be found in, e.g., [28], [16] or [5].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Initial density profile (a) and velocity profile (b) of the runs. $\delta_i$ is half the (so-called) vorticity thickness, $U$ and $\Delta\rho/\sqrt{Pr}$ respectively are half the velocity and density difference across the horizontal boundaries; $Pr$ is the Prandtl number, equal here to 0.7 (from [32]).</figcaption></figure>

## 3.2 Mathematical model and numerical method

To investigate the shear layer dynamics and resulting mixing properties, we solve the three-dimensional Boussinesq equations:

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u = -\nabla P + J \rho + \frac{1}{Re} \nabla^2 u \quad (28) $$

$$ \frac{\partial \rho}{\partial t} + (u \cdot \nabla)\rho = \frac{1}{PrRe} \nabla^2 \rho \quad (29) $$

$$ \nabla \cdot u = 0. \quad (30) $$

$u = (u, v, w)$ denotes the velocity field with $u$, $v$ and $w$ being respectively the streamwise, spanwise and vertical components; $P$ and $\rho$ are the dynamical pressure deviation and the density deviation about a hydrostatic balanced state of reference density $\rho_0$; $v$ is the kinematic viscosity.

The initial velocity and density fields, $U_i(z)$ and $\rho_i(z)$ respectively, are of the error function type, the thickness of the density profile being $\sqrt{Pr}$ larger than that of the velocity profile (Fig. 3). A two-dimensional perturbation is superimposed upon $U_i$ to destabilize the flow (the Kelvin-Helmholtz instability is thus able to grow). Three-dimensional perturbations are also added to create three-dimensional motions. The boundary conditions are periodic along the streamwise ($x$) and spanwise ($y$) directions and of free slip type along the vertical ($z$).

The Boussinesq equations (28–30) are written in non-dimensional form. We take half the thickness of the initial velocity profile as a length scale, which we denote as $\delta_i$, half the velocity difference across the flow as a velocity scale ($U$) and half the density difference across the flow times $\sqrt{Pr}$, as a density scale
($\Delta\rho\sqrt{Pr}/2$). The time unit is $\delta_i/U$. This scaling sets three parameters $J$, $Re$ and $Pr$, which control the flow dynamics. They are defined by:

$$J = \frac{g \Delta \rho \sqrt{Pr} \delta_i}{\rho_0 2 U^2} \qquad Re = \frac{U \delta_i}{\nu} \qquad Pr = \nu/\kappa. \qquad (31)$$

A pseudo-spectral method is used to solve the Boussinesq equations: the equations are solved in Fourier space except for the nonlinear terms, which are computed in physical space; aliasing errors are removed by a truncation scheme; symmetry conditions are imposed along the z-direction so that the fields can be expanded in cosine or sine series along this direction. A third-order Adams-Bashforth scheme is used to advance the equations in time.

We have performed several runs in which $J$ ranges from 0.04 to 0.16 and $Pr$ has a value of 0.7 or 1.4. The resolution is $(N_x, N_y, N_z) = (128, 128, 257)$ and the domain size is $(\lambda, \lambda, 2\lambda)$, where $\lambda$ is the most amplified wavelength predicted by linear theory ($\lambda$ is nearly independent upon $J$, as long as $J < 1/4$ [12]). Hence, only one Kelvin-Helmholtz vortex can develop in the numerical domain. The Reynolds number $Re \approx 300$. The approximate value for $Re$ stems from the fact that a uniform but time varying viscosity is used. Indeed, in order not to have the flow dominated by both stratification and viscous effects at the end of the runs, the viscosity is allowed to vary in time according to the constraint that the Kolmogoroff scale $\eta = (\nu^3/\epsilon)^{1/4}$ remains equal to the grid size at any time. Computing $\epsilon$ (the dissipation rate of kinetic energy) at each time step provides the value of $\nu$ at that time. In practice, $\nu$ smoothly decreases by a factor 3 from its initial value (equal to 1/300) up to the time $\epsilon$ has reached a maximum value (the flow has become strongly turbulent) and eventually reaches a value lower than $\nu(0)$ at the end of the run (equal to 1/500).

An example of the flow development is displayed in Fig. 4. At $t = 34$ (Fig. 4a), the Kelvin-Helmholtz instability has grown and saturated: it is manifested as one Kelvin-Helmholtz vortex, already distorted by the growth of the three-dimensional perturbation. The vortex is located in between sheets, which are the usual braids reinforced by the baroclinic torque (see [7], [30]). The saturation of the Kelvin-Helmholtz instability is followed by the growth of a three-dimensional secondary instability. Fig. 4b displays the total vorticity field at the time the dissipation rate of kinetic energy $\epsilon$ is maximum: the three-dimensional instability has led to small scale structures, whose dominant vorticity component lies along the streamwise direction. This instability is thus of the Rayleigh-Taylor type, having occurred when heavy fluid lies over light fluid ([16]). The end of the run is visualized in Fig. 4c: the central part of the shear layer is relaminarizing and the turbulent activity, while decaying, is expulsed further outward of the central region.

In the following, the Kelvin-Helmholtz vortex will be referred to as a *large* scale structure, as opposed to the *small* scale three-dimensional structures. These adjectives may be quantified by the Reynolds number associated with each structure: $Re \approx 800$ in the former case and $Re \approx 25$ in the latter.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Temporal evolution of a three-dimensional shear layer simulation with $J = 0.125$, $Re \approx 1/300$ and $Pr = 1.4$. A surface equal to $0.5 \max_D(|\omega_x| + |\omega_y| + |\omega_z|)$, where $|\omega_i|$ is the $i$-component of the vorticity and $D$ the fluid domain, is shown at successive times: (a) $t = 34$. (b) $t = 108$. (c) $t = 227$, at the end of simulation; only 3/4 of the numerical domain along the vertical direction is shown for clarity. The time unit is $\delta_i/U$ ($\delta_i$ and $U$ are defined in section 3.2) (from [32]).</figcaption></figure>

## 3.3 Analysis of mixing

The mixing properties of a stably stratified shear layer have been studied by a few authors. To our knowledge, the only laboratory experiments are those of Thorpe [34] and Koop & Browand [17]. All the other works are numerical ([29], [9], [31], [32], [6]). The method derived by Winters et al. [38] is used in the two
latter papers only (we still note that background density profiles are computed
in [29]).

* *Computational method of the background density profile*

The method to compute the background density profile $\rho_b$ from the three-
dimensional density field is as follows. The fluid domain is partitioned into el-
ementary volumes of side equal to the grid size. Each elementary volume rep-
resents one fluid particle. These fluid particles are next sorted according to the
value of their density. Two equivalent methods can be used to construct the
stable density profile. In the first method ([38], [32]), the profile is constructed
by filling the successive horizontal planes with the ordered particles, the bottom
plane being filled with the heaviest particles. The value of the density slightly
varies from particle to particle in a given plane (by a few percents at most in
our computations) so that the value of the background density at each plane is
computed by a horizontal average. In the second method ([6]), each particle is
stretched horizontally so as to fill successively each horizontal plane; no horizon-
tal average is therefore required. In numerical computations however, the density
field needs to be interpolated if combined with other quantities, to match the
vertical resolution.

* *Horizontally averaged and background density profiles*

The horizontally averaged and background density profiles, $\bar{\rho}(z)$ and $\rho_b(z)$ respectively, are plotted in Fig. 5 at the same key times as in Fig. 4, but for a different run (see caption). Figure 5a shows that $\rho_b(z)$ is stable, consistently with its definition, and remains close to the initial density profile. By contrast, $\bar{\rho}(z)$ is not monotonic and strongly differs from $\rho_b(z)$. This indicates that (i) the large scale Kelvin-Helmholtz vortex hardly mixes the fluid (because $\rho_b \approx \rho_i$); (ii) a horizontal average does not remove reversible effects (because $\bar{\rho} \neq \rho_b$). Hence, the use of the horizontally averaged density profile is not a good approximation to the background density profile. In Fig. 5b, the two profiles are displayed when ubiquitous small scale three-dimensional motions have developed. The striking feature of the figure is the broadening of the sorted density profile: this reveals the dominant contribution to fluid mixing by the three-dimensional small scale structures. The horizontally averaged density profile displays the same trend but it fluctuates about the background profile and its gradient changes sign. Only at the end of the run (Fig. 5c), when the flow starts to relaminarize in its central part and turbulence decays in its outer part, do both profiles coincide. They have nearly become linear and strongly depart from the initial profile.

* *Turbulent mixing rate*

The temporal evolution of the turbulent mixing rate $\Phi_d - \Phi_{lam}$ offers a syn-
thetic view of the respective contribution to mixing by the Kelvin-Helmholtz
vortex and the three-dimensional structures. Fig. 6 thus shows that $\Phi_d - \Phi_{lam}$
displays a relative maximum at $t = 35$, when the Kelvin-Helmholtz instability
saturates; $\Phi_d - \Phi_{lam}$ further increases when three-dimensional instabilities de-
velop and, at the time turbulence is the strongest ($t \approx 85$), reaches an absolute
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Horizontally averaged $\bar{\rho}(z)$ (dashed line) and background $\rho_b(z)$ (full line) density profiles. In (a) and (c), the initial density profile is plotted for comparison (dotted line). (a) $t=28$ (when the Kelvin-Helmholtz vortex has developed); (b) $t=108$ (when three-dimensional small scale structures have formed); (c) $t=227$ (end of simulation). The parameters of this run are $J = 0.04$, $Pr = 0.7$, $Re \approx 300$ (from Staquet [32]).</figcaption></figure>

maximum which is nearly three times higher than the relative maximum (note that $\Phi_d/\Phi_{lam} \approx 20$ when these maxima are reached).

### *Mixing efficiency*

The previous results demonstrate that mixing is promoted by small scale three-dimensional structures, and not by a large scale Kelvin-Helmholtz vortex. But what about the efficiency of mixing? To address this important point, and in order to interpret it properly, we first present in Fig. 7 the dissipation rate of
<figure><img src="image_7.png" /><figcaption>Fig. 6. Turbulent diapycnal flux $\Phi_d - \Phi_{lam}$ versus time for the same run as in Fig. 5.</figcaption></figure>

kinetic energy $\epsilon$ for the same run as before. The behavior of $\epsilon$ is not surprising: it remains close to its initial value when the Kelvin-Helmholtz instability grows and saturates, because this instability, and the associated Kelvin-Helmholtz vortex, are two-dimensional. As is well-known (e.g. [20]), the material conservation of vorticity in a two-dimensional flow prevents energy from being transferred toward small scales, where dissipation is acting. By contrast, Fig. 7 shows that the growth of a three-dimensional instability, and its nonlinear manifestation as small scale structures, makes $\epsilon$ strongly increase, up to a value equal to nearly 20 times its initial value.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Dissipation rate of the turbulent kinetic energy, $\epsilon$, versus time for the same run as in Fig. 5.</figcaption></figure>
The behavior of the mixing efficiency or, equivalently, of the flux Richardson number $Rf_b$ can be immediately inferred from Figs. 6 and 7. Figure 8, where $Rf_b$ is displayed versus time, thus shows that the efficiency of mixing is high (close to 0.7) during the initial two-dimensional stage of the flow and moderate (close to 0.3) when three-dimensional secondary structures have developed. In other words, the large scale Kelvin-Helmholtz vortex contributes weakly to mixing, but in an efficient manner (because nearly no energy is dissipated) while the three-dimensional structures are responsible for the mixing of the fluid, but much less efficiently (because they occur at scales close to dissipative scales). Hence, the Kelvin-Helmholtz vortex mainly stirs the fluid while the three-dimensional structures effectively mix it.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Flux Richardson number $Rf_b$ averaged over seven runs differing by the value of $J$ (from 0.04 to 0.16) and $Pr$ (equal to 0.7 or 1.4) (from [33]).</figcaption></figure>

• *A parameterization of mixing: turbulent diffusivity versus Richardson number*

As said in the Introduction, the parameterization of mixing consists of a func-
tional dependency of the diapycnal diffusivity $K_b$ (defined by Eq. (16)) upon a
dynamical parameter characterizing the largest scales; in a general circulation
model, the latter scales belong to the resolved scales. Such a functional depen-
dency may be seeked for in the numerical simulations of a stably stratified shear
layer we presented because the largest scale is a shear flow (that is, with no
horizontal scale) and is well separated from the scale of the three-dimensional
motions that mix the fluid. One may thus assume the former motion to be as-
sociated with a current that would lie among the resolved scales in an oceanic
circulation model while the latter motions are at subgrid scale.

The presence of a mean shear implies that a Richardson number should be
taken for the dynamical parameter. We make use of the monotonicity of the
background density profile $\rho_b$ to introduce a Richardson number that remains
one-signed at any time:

$$Ri_b(z) = -J \frac{d\rho_b/dz}{(d\bar{U}/dz)^2} \qquad (32)$$

The normalized diapycnal diffusivity $K_b(z)/\kappa$ is plotted as a function of $Ri_b(z)$ in Fig. 9a, when the turbulence has started to decay (for $t \ge 196$; no clear dependency could be found before this time). Two parts may be distinguished in Fig. 9a, where a log-log scale is used: for $Ri_b < 0.2$, $K_b(z)/\kappa$ decays quasi-linearly with $Ri_b(z)$. This decaying behavior is expected: when $Ri_b$ decreases, stratification effects weaken relative to nonlinear effects so that the flow becomes more turbulent and mixing increases. This part of the curve is contributed by the outer region of the shear layer, where turbulence remains. Figure 9a shows that values of the turbulent diffusivity $K_b$ up to 200 times the molecular diffusivity are found in this region. For $Ri_b(z) > 0.2$ by contrast, that is, in the central relaminarizing part of the shear layer, no clear dependency is observed.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Diapycnal diffusivity $K_b(z)$ (defined by Eq. (16)) normalized by the molecular diffusivity $\kappa$ as a function of a gradient Richardson number $Ri_b(z)$ (defined by Eq. (32)), for $t \ge 196$ (each curve corresponds to a different time). (a) Log-log scale; (b) Lin-lin scale. For this run, $J = 0.125$, $Pr = 1.4$.</figcaption></figure>

The quasi-linear dependency of $K_b(z)/\kappa$ with $Ri_b(z)$ displays a -2 slope. This is confirmed in Fig. 9b, where $K_b(z)/\kappa$ is plotted versus $1/Ri_b^2(z)$: the linear dependency is recovered but the slope varies in time. In the ocean, below the upper mixed layer, several parameterizations are used depending upon the process that creates the turbulence: shear instability, internal gravity wave breaking and double diffusion (due to the simultaneous diffusion of heat and salt). For the former process, a parameterization of the form $K/K_{max} = [1 - (Ri/Ri_0)^2]^3$
is commonly used in large scale circulation models for $0 < Ri < Ri_0$, with $K_{max} = 5 \ 10^{-3} m^2/s$ and $Ri_0 = 0.7$. The high value of $Ri_0$ stems from the fact that the measurements of $Ri$ rarely fall as low as the theoretical value of 1/4. According to Large et al. [18], this parameterization is supported by limited observations ([27]). Our numerical results suggest that this is a too sharp dependency.

# 4 Conclusion

As discussed in this article, basic concepts of mixing have been known for a long time (dating from Lorenz [22]) and are now properly formalized. The parameterization of mixing still remains a difficult problem. Parameterizing implies a fundamental understanding between mixing and flow dynamics. The existence of a link between mixing, which occurs at dissipative scales and the flow dynamics, whose energy reservoir is at large scale, is not surprising: mixing is controlled by these large scales, as is kinetic energy dissipation in (three-dimensional) ordinary turbulence. In the latter case, the problem is simplified by the fact that, at high enough Reynolds number, the kinetic energy dissipation rate depends only upon the rate of input of energy at large scale, and not upon the details of forcing (assumed to be isotropic). This is not true when a stable stratification exists: the memory of forcing is all the more important the higher the stratification is. The details of forcing is most often associated with an instability that triggers the (kinetic and potential) energy transfer toward small scale. To derive the parameterization, the source of the instability (an unstable shear flow for instance) should rather be considered. Computing a dynamical parameter from the characteristics of this source may indeed lead to a more general behavior than considering the instability itself, as is commonly done. This approach would also provide a usable parameterization of mixing, because the source of turbulence belongs to the resolved scales in large scale circulation models, while the instability that leads to fluid mixing does not in general.

# Acknowledgments

I thank A. Beckman for having kindly provided Fig. 1 and A.-M. Treguier for useful discussions about parameterization of mixing in the ocean.

# References

1. D.G. Andrews: *J. Fluid Mech.* **107**, 227 (1981)

2. A. Beckmann: 'Turbulent mixing in the ocean: intensity, causes and consequences'. In: *Proceedings of Les Houches Winter School, France, January 1998*, ed. by E. Chassignet, J. Verron (Kluwer Academic Publishers, 1998) pp. 135-154

3. F. Bryan: *J. Phys. Oceanogr.* **17**, 970 (1987)

4. D.R. Caldwell, J.N. Moum: *Rev. Geophys. Suppl.*, **1385** (1995)
5. C.P. Caulfield, W.R. Peltier: *Phys. Fluids* **6**(12), 3803 (1994)

6. C.P. Caulfield, W.R. Peltier: *J. Fluid Mech.* **413**, 1 (2000)

7. G.M. Corcos, F.S. Sherman: *J. Fluid Mech.* **73**, 241 (1976)

8. G.M. Corcos, F.S. Sherman: *J. Fluid Mech.* **139**, 29 (1984)

9. A.B. Cortesi, B.L. Smith, G. Yadigaroglu, S. Banerjee: *Phys. Fluids* **11-1**, 162 (1999)

10. P.G. Drazin, W.H. Reid: *Hydrodynamic stability* (Cambridge University Press, Cambridge 1981)

11. M.C. Gregg: *J. Geophys. Res.* **92-C5**, 5249 (1987)

12. P. Hazel: *J. Fluid Mech.* **51**, 39 (1972)

13. H. Ho, P. Huerre: *Ann. Rev. Fluid Mech.* **16**, 365 (1984)

14. D. Holliday, M.E. McIntyre: *J. Fluid Mech.* **107**, 221 (1981)

15. R.X. Huang: *J. Phys. Oceanogr.* **29**, 727 (1999)

16. G.P. Klaassen, W.R. Peltier: *J. Fluid Mech.* **227**, 71 (1991)

17. C.G. Koop, F.K. Browand: *J. Fluid Mech.* **93**, 135 (1979)

18. W.G. Large, J.C. McWilliams, S.C. Doney: *Reviews of Geophysics* **32-4**, 363 (1994)

19. J.R. Ledwell, A.J. Watson, C.S. Law: *Nature* **364**, 701 (1993)

20. M. Lesieur: 'Turbulence in fluids'. Kluwer Academic Pub. (1990)

21. P.N. Lombard, D.D. Stretch, J.J. Riley: 'Energetics of a stably stratified mixing layer'. In: *Ninth Symposium on Turbulence and Diffusion, Roskilde, Danemark, April 30-May 3, 1990*, (A.M.S. Publ. 1990), pp. 202-206

22. E.N. Lorenz: *Tellus* **7**, 157 (1955)

23. W. Munk: *Deep-Sea Research* **13**, 707 (1966)

24. W. Munk, C. Wunsch: *Deep-Sea Research I* **45**, 1977 (1998)

25. T.R. Osborn: *J. Phys. Oceanography* **4**, 83 (1980)

26. T.R. Osborn & C.S. Cox: *Geophys. Fluid Dyn.* **3**, 321 (1972)

27. H. Peters, M.C. Gregg, J.M. Toole: *J. Geophys. Res.* **93**, 1199 (1988)

28. D.G. Schowalter, J.C. Lasheras, C.W. van Atta: *J. Fluid Mech.* **281**, 247 (1994)

29. J.F. Scinocca: *J. Atmos. Sciences* **52-14**, 2509 (1995)

30. C. Staquet: *J. Fluid Mech* **314** (1995)

31. C. Staquet, K.B. Winters: 'Mixing in a stably-stratified shear layer'. In: *Proc. of the XIth Symposium on Turbulent Shear Flows, Grenoble, France, September 1997*, ed. by G. Binder, vol. 2, session 20, p. 25-30

32. C. Staquet: *Fluid Dyn. Research* **27**, 367 (2000)

33. C. Staquet, P. Bouruet-Aubertot: *Dyn. Atmos. Oceans* **34**, 81 (2001)

34. S.A. Thorpe: *J. Fluid Mech.* **61**, 731 (1973)

35. S.A. Thorpe: *Phil. Trans. R. Soc. Lond.* **A 286**, 125 (1977)

36. S. Thorroddsen, C.W. van Atta: *J. Geophys. Res.* **97-C3**, 3647 (1992)

37. J. Toole: 'Turbulent mixing in the ocean: intensity, causes and consequences'. In: *Proceedings of Les Houches Winter School, France, January 1998*, ed. by E. Chassignet, J. Verron (Kluwer Academic Publishers, 1998) pp. 171-190

38. K.B. Winters, P.N. Lombard, J.J. Riley, E.A. D'Asaro: *J. Fluid Mech.*, **289**, 115 (1995)

39. K.B. Winters K., E.A. D'Asaro: *J. Fluid Mech.* **317**, 179 (1996)

40. H. Yamazaki, T. Osborn: *J. Geophys. Res.* **95-C6**, 9739 (1990)