Mixing in Turbulent Overturns
W.D. Smyth College of Oceanic and Atmospheric Sciences Oregon State University 

Introduction *
1. The evolution of turbulence in localized shear *
2. The state of maximum turbulence *
3. The decay phase *
4. Do these results represent ocean turbulence? *
5. Mixing in laminar and turbulent flow *
6. Implications for the broader study of turbulent mixing *
Bibliography *
The past decade has brought tremendous insights into the physics of turbulence, due largely to direct numerical simulations (DNS). This new understanding applies almost entirely to the simplest idealization, i.e. stationary, homogeneous, isotropic turbulence. In nature, turbulence never conforms to this simple picture. In particular, geophysical turbulence is almost always affected by ambient shear and density stratification, which complicate the physics greatly. The turbulence theory and modeling program at COAS aims to extend stateoftheart theories of turbulence to smallscale geophysical flows by accounting for these effects.
Our primary tool in recent years has been direct numerical simulation of turbulent overturns driven by KelvinHelmholtz (KH) instability of an inflectional shear profile. It has often been suggested that this mechanism is what generates the turbulent "patches" that are observed to account for most of the mixing in the ocean thermocline. In this review, I summarize the nature and evolution of turbulence generated via the growth and collapse of KH billows.
Section 1 contains an overview of the evolution of an unstable parallel shear layer through stages of billow growth, turbulence, and ultimate decay to a stable parallel state. Section 2 compares the state of maximal turbulence with current understanding of high Reynolds number flow. Section 3 describes the decay of turbulence. In section 4, I demonstrate that this process represents a valid model for turbulent overturns in the thermocline. Section 5 contains further detail on mixing physics. There, I reconsider the traditional link between mixing and turbulence, and show that understanding of mixing must include consideration of coherent structures. Implications for other types of ocean turbulence are discussed in section 6.
1. The evolution of turbulence in localized shear
The essential parameter governing sheared, stratified turbulence is the gradient Richardson number, defined as the squared ratio of the ambient buoyancy frequency to the ambient shear. In general, the latter two quantities vary in the vertical. It is then convenient to define a bulk Richardson number, , based on values that characterize the region of interest (typically a turbulent layer). In the simplest case in which shear and stratification are uniform and turbulence is homogeneous, the bulk Richardson number, is constant in time. The energy of the turbulence grows in time if is less than a critical value, , which is near 1/4. If >, turbulence decays (e.g. Itsweire et al., 1993).
When turbulence grows at a local shear maximum (i.e. an inflection point of the background velocity profile), increases in time as turbulent entrainment causes the sheared zone to spread (figure 1a). Again, there is a critical value that is near 1/4. In this case, however, eventually grows to exceed and turbulence decays. Turbulence in stratified flow is quantified using the buoyancy Reynolds number, (figure 1b). grows to a maximum as approaches , then decreases. Ultimately, approaches a value near 1/3, beyond which entrainment ceases and further spreading of the shear layer is driven only by the action of viscosity on the mean flow (figure 1).
The growth of turbulence at a shear maximum begins with the coalescence of the mean flow vorticity into a periodic train of KH billows (figure 2a). The billows break and become turbulent via a combination of convective and sheardriven secondary instabilities (figure 2b; Klaassen & Peltier, 1991). These will include one or more vortex pairings if the initial value of is sufficiently small. Beyond this point, turbulence decays slowly, and the flow eventually relaxes to a stable, parallel state perturbed by gravity waves (figure 2d).
Figure 1: Evolution of (a) the bulk Richardson number and (b) the buoyancy Reynolds number in a typical simulation. Time is measured in buoyancy periods, i.e. . Horizontal lines represent (a) the approximate critical Richardson number, 1/4, and (b) the critical value =20 for mixing in stratified turbulence. 

Figure 2: Isocontours
of the enstrophy (squared vorticity) at four instants in the life
of a turbulent overturn generated by KH instability.
2. The state of maximum turbulence
When is a maximum, smallscale flow structures exhibit many properties characteristic of turbulence, as expected both from the theory of high Reynolds number flow and from observations of strongly turbulent flow in nature. The probability density function for the turbulent kinetic energy dissipation rate, epsilon, is nearly lognormal, with skewness slightly negative (figure 3).


Figure 3: Probability density of the turbulent kinetic energy dissipation rate at selected instants in the evolution of a turbulent overturn. The state of maximum occurs at =3.4. At =4.2, the darker area represents the statistics of the turbulent region of the flow while the lighter area corresponds to the surrounding laminar regions. 
Figure 4: Second invariant of (a) the dissipation anisotropy tensor and (b) the scalar gradient anisotropy tensor versus buoyancy Reynolds number. Results shown represent several experiments started from different initial conditions. Colors indicate time as in figure 3. 
Although the velocity field is highly anisotropic in consequence of the wavelike nature of the large scales, velocity gradients are nearly isotropic (figure 4). In contrast, temperature gradients exhibit significant anisotropy even at high . This echoes the recent discovery that scalar fluctuations in the presence of a mean gradient are anisotropic, even if the velocity field is isotropic (Sreenivasan 1991).
Velocity spectra show an inertial subrange extending over about one half decade in these simulations. At small scales, the spectra are isotropic and correspond closely to the form derived from measurements in a tidal channel at very high Reynolds number by Nasmyth (Oakey, 1982).
Figure 5: (a) Transverse velocity spectra taken in the streamwise direction in a state of strong turbulence. Blue and green curves correspond to the vertical and spanwise velocity components. The red curve is the spectrum of streamwise velocity transformed into transverse form using the standard isotropic relation. The horizontal line indicates the height of an inertial subrange with Kolmogorov constant equal to 1.6. The black curve is the Nasmyth spectrum, the standard against which oceanographic spectra are compared. (b) Streamwise spectrum of the streamwise scalar gradient (red). The black curves indicate the theoretical forms proposed by Batchelor (1959) and Kraichnan (1968). 

Scalar gradients conform to the theoretical form predicted by Kraichnan (1968) on the basis of the Lagrangian History Direct Interaction Approximation (LHDIA), with Batchelor constant q=7 (figure 5b). This spectral shape is also characteristic of temperature spectra in the ocean (Hill 1978; Smyth, Nash & Moum, unpublished results). The value of q is considerably large than the value 2.0 predicted by Batchelor. The discrepancy is a consequence of the relative ineffectiveness of turbulent strain fields at sharpening scalar gradients, an issue to which we will return in section 5.
All of the properties of the state of maximal turbulence listed in the previous subsection change as the turbulence decays. Some of these changes can be understood as results of the decreasing spectral separation between the production and dissipation subranges. For example, the turbulent kinetic energy dissipation rate loses its lognormal character as the Kolmogorov scale increases to the point where the mean shear exerts a dominant influence (figure 3). The degree of anisotropy in the velocity and temperature gradients increases (figure 4) as the fluctuations decay to a magnitude comparable to the highly anisotropic gradients of the mean flow. The separation between large and small scales can be quantified using the buoyancy Reynolds number, which is the ninefourths power of the ratio of Ozmidov to Kolmogorov length scales. The dependence of velocity gradient anisotropy on the range of scales is demonstrated in figure 4(a) by plotting anisotropy as a function of buoyancy Reynolds number for several runs with different initial conditions and observing how the results collapse.
Not all properties of decaying turbulence are governed by the range of scales. Both , the ratio of Ozmidov to Thorpe scales, and the flux coefficient are primarily a function of event age, i.e. they follow a common evolutionary pattern for overturns with different Reynolds numbers. These properties are discussed in detail below.
Properties that depend on are primarily determined by the smallest scales of the flow. As anticipated by Kolmogorov (1941) and others, the small scales exhibit a structure that is nearly universal, provided that they are insulated from the large scales by an inertial subrange. In that case, the small scales receive energy from the large scales, but do not receive any information of their structure. The large scales, in contrast, have both spatial structure and a pattern of temporal evolution that is highly asymmetrical and is characteristic of the particular initial conditions from which the flow evolved. Thus flow properties that depend specifically on largescale quantities will reflect their particular evolutionary pattern.
4. Do these results represent ocean turbulence?
If these DNS results are to be useful in understanding turbulence as measured in the ocean, it must be established that the model is modeling the correct process. The question of model validity comes in three parts: (1) "Does the numerical model solve the NavierStokes equations plus initial and boundary conditions accurately?", (2) "Do the simulated flows achieve Reynolds numbers high enough to be characteristic of turbulence?", and (3) "Is the simulated flow representative of turbulence as it occurs in the ocean?".
The first question has been addressed by testing conservation properties, by checking spectra for adequate spatial resolution, and by comparing results with laboratory experiments. The second question is answered in part by the magnitude of the buoyancy Reynolds number: the criterion >20 for actively mixing turbulence is well satisfied. More direct confirmation is found in the existence of a short but distinct inertial subrange (figure 5a) and in the near isotropy of the velocity gradients.
The third question is the most difficult to answer owing to the extreme sparseness of ocean turbulence observations. Figure 6 shows the evolution of several turbulence statistics in three different DNS simulations. To the right of each frame is a histogram showing the statistics of the same quantity drawn from 3425 profiles through turbulent overturns in main thermocline off Northern California.

Figure 6: Evolution of turbulence statistics in selected DNS simulations, compared with statistics drawn from observations of turbulent overturns in the ocean. Quantities compared are (a) the buoyancy Reynolds number, (b) the ratio of maximum to RMS Thorpe displacement, (c) the ratio buoyancy scale to Thorpe scale and (d) the ratio of Ozmidov to Thorpe scale. The histogram to the right of each frame gives the corresponding statistics drawn from observations in the main thermocline. 
The buoyancy Reynolds number tends to lie at the low end of the range of the observations. This result illustrates the main limitation of DNS: inability to resolve a wide range of length scales. It is significant, however, that computer size has grown to the point where realistic Reynolds number can be achieved. Also shown are experiments in which the Prandtl number was lowered in order to achieve higher . The ratio of maximum Thorpe displacement to Thorpe scale (figure 8b) characterizes the largescale geometry of the billows. This ratio remains within a very small range of values that coincides closely with the observations. A similarly close correspondence of the ratio of buoyancy scale to Thorpe scale (figure 6c) indicates that the partitioning of kinetic and potential energies in the large scales is consistent with observations. Finally, the ratio of Ozmidov to Thorpe scale (figure 6d) indicates the relationship between the energycontaining scales, the dissipation scales and the background stratification. Again we see a close correspondence with the observational statistics, although the DNS values often lie at the limits of the observed range. The latter result indicates that observational sampling methods miss both the youngest and the oldest events.
These comparisons, and many others in the same vein, have shown that the simulated overturns cannot be distinguished from the observed overturns on the basis of presently available observational data, except that the Reynolds numbers lie at the low end of the observed range. In other words, the data is insufficient to disprove the hypothesis that the observed overturns are KelvinHelmholtz billows.
5. Mixing in laminar and turbulent flow
Figure 7(a) shows the evolution of the total and background potential energies of the flow (denoted and ). The latter quantity is the minimum potential energy achievable by spatial reordering of fluid parcels. Changes to can only be positive, and occur only as a result of mixing. The total potential energy includes both and available potential energy residing in waves and overturns. The early maximum in corresponds to the initial rollup and pairing of the KelvinHelmholtz billows. The subsequent release of available potential energy indicates the turbulent breakdown of the billow (also shown by the gray bar at the top of the figure). Eventually, waves and turbulence decay and the available potential energy drops to zero.

Figure 7: (a) Time evolution of the total and background potential energies in a simulated KH billow. (b) Time evolution of the flux coefficient: , where represents diffusion due to the background state. 
Figure 7(b) shows the evolution of the flux coefficient, , which is defined as the ratio of (adjusted slightly to remove effects of mean flow diffusion) to the TKE dissipation rate, epsilon. is thus closely related to mixing efficiency. In the initial rollup phase, mixing is highly efficient, i.e. significant mixing is occurring while very little energy is being dissipated by friction. then decreases during the transition phase, and finally asymptotes to a value very close to 0.2, the value that is traditionally assumed in the analysis of oceanographic observations.
The high efficiency of mixing during the preturbulent phase challenges our traditional association of mixing with turbulence. Our expectation was that mixing would commence suddenly with the onset of turbulence. This expectation is not entirely false: the high mixing efficiency at early times is due as much to the absence of viscous dissipation as to the presence of mixing. However, inspection of figure 7(a) shows that, by the time the transition to turbulence is complete, about onehalf of the total mixing has already occurred! About one quarter of the net mixing occurs before transition while the flow is still twodimensional.
It is possible that the foregoing result is a function of Reynolds number. At higher Re, the decay phase could last longer and therefore accomplish a larger fraction of the net mixing. This is not obvious, though, since the decay of turbulence is controlled largely by stratification. In the oceanographic context, shear instability is generally driven by transient features, such as internal wave superpositions, so that the duration of the event is once again be limited by factors other than the Reynolds number. In summary, there is no reason to conclude that the relative importance of the preturbulent phase to the net mixing accomplished by overturns is different in general from that suggested in figure 7.
This result has significant implications for observations of ocean mixing. In observational studies of turbulent overturns (the primary mixing process over most of the ocean interior), regions that show no evidence of turbulence are normally ignored. For example, sample selection criteria employed by Moum (1996) effectively filtered out preturbulent billows. Moreover, close inspection of the DNS results shows that most of the mixing in preturbulent overturns occurs in the braids connecting the billows. Observational profiles through these regions would not appear to be associated with overturning at all; in fact, the stratification would be strongly stable!
To understand the extraordinary efficiency of mixing in preturbulent overturns, we must think more closely about the physics of mixing. Mixing occurs entirely through the action of molecular diffusion as it works to reduce gradients of fluid properties. When we assert that turbulence "causes" mixing, we mean that turbulence amplifies gradients, thus accelerating the molecular mixing process. Gradient amplification occurs as a result of compressive strain in the velocity field. The association of strong mixing with turbulence is based on the correct observation that velocity strains are generally much larger in turbulent flow than in laminar flow. However, gradient amplification depends on more than just the magnitude of the strain; it also depends on the orientation between the strain and the scalar gradient. Gradient amplification is most effective when the gradient is aligned parallel to the compressive principal strain. In the opposite extreme, where gradients are aligned randomly with respect to the strain, there is no net amplification at all, however intense the strain may be.
So, turbulent mixing would not occur at all if the gradientstrain alignment was random, and this is clearly not the case. In fact, misaligned gradients tend to rotate toward the direction of perfect alignment over time. If the strain field is steady, the scalar gradients will soon attain nearly perfect alignment. This is exactly what was assumed in early theories of scalar mixing (e.g. Batchelor 1959?): that the turbulent strain field evolves slowly enough that scalar gradients could be assumed to be in this ideal equilibrium state. We now know, however, that this assumption is incorrect: turbulent strain fields evolve far too rapidly for the scalar gradient to maintain anything close to ideal alignment. As a result, turbulent mixing is much less efficient than one would expect based on the magnitude of the strain alone. This is why Batchelor's constant, q, is much larger than the value 2.0 predicted by Batchelor (cf. section 2).
The situation in preturbulent KH billows is very different. While the strain field is relatively weak, it also evolves slowly, so that the scalar gradient attains nearly perfect alignment with the direction of compressive strain. This highly effective compression of scalar gradients occurs mainly in the braids connecting the billows (figure 8). Typically, scalar gradients in those regions are as strong or stronger than those occurring after the transition to turbulence.
Figure 8: Crosssection of the temperature field in KelvinHelmholtz billows before (a) and after (b) the transition to turbulence. The coherent strain field in the preturbulent case generates sharp gradients in the braids between the billows. In the turbulent case, strains are stronger but less coherent, so that mixing efficiency is reduced.
Is the high mixing efficiency of young overturns evident in observational data? That question is difficult to answer because there is no direct way to determine the "age" of an observed overturn. We therefore turn once again to the DNS results for guidance. It has been suggested many times that nondimensional ratios of characteristic length scales of turbulence evolve monotonically in time, essentially because they quantify the geometrical complexity of the flow. As a result, such a ratio may serve as a "clock" indicating the age of an overturn. However, this assumption cannot be checked without some independent estimate of event age (e.g. Wijesekera & Dillon 1997). Here, DNS results provide the needed information. In figure 9, we show that a particular ratio R, involving the Thorpe, Ellison and Ozmidov scales, functions effectively as an observational clock.


Figure 9: Time dependence of the dimensionless length scale ratio in selected simulations. The solid curve indicates the powerlaw approximation . 
Figure 10: Dependence of (approximated using the OsbornCox formula) upon R in observational data from the main thermocline. Each point represents a bin average of about 400 profiles through turbulent overturns. 
In figure 10, we plot estimates of the flux coefficient for observed turbulent overturns as a function of R. As in the simulation, the observed values of range between zero and unity and generally decrease over time. (The fact that the observed values do not asymptote to 0.2 appears to be a result of noise in the measurements of epsilon, which leads to underestimates of in weak turbulence.)
When this dataset was reanalyzed using the revised value for in place of 0.2 the result was a 50% increase in the overall estimate of the turbulent scalar flux in the main thermocline. When the analysis is revised further to include nonoverturning regions (e.g. the braids connecting preturbulent billows), the estimate is expected to increase further.
6. Implications for the broader study of turbulent mixing
The present results are part of a growing realization by the turbulence community of the importance of coherent structures to turbulent mixing. Mixing does not occur randomly throughout a turbulent flow, as was once assumed. In fact, it is concentrated in the strained regions that connect coherent volumes of strong vorticity. It is therefore essential that the dynamics of coherent structures be understood if mixing is to be accurately predicted and parameterized.
While the KH vortices discussed here represent the energycontaining scales of a specific flow geometry, their physics is similar to that of sheardriven vortices in the dissipation range, whose characteristics are to a large degree universal. It follows that the essence of the present results, i.e. that efficient mixing occurs in association not with highly chaotic flow but with coherent structures, is likely to be applicable in all types of turbulence. These results point the way to a fuller understanding not only of turbulent overturns in the main thermocline, but of turbulence in the many forms that it takes in the Earth's oceans.
This report summarizes results from the following publications:
Smyth, W.D., 1999,"Dissipation range geometry and scalar spectra in sheared, stratified turbulence", J. Fluid Mech. 401, 209242.
Smyth, W.D. & J.N. Moum 2000a, "Length scales of turbulence in stably stratified mixing layers", Phys. Fluids 12, 13271342.
Smyth, W.D. & J.N. Moum 2000a, "Anisotropy of turbulence in stably stratified mixing layers", Phys. Fluids 12, 13431362.
Smyth, W.D., J.N. Moum & D.R. Caldwell 2000, "The efficiency of mixing in turbulent patches: inferences from direct simulations and microstructure observations", J. Phys. Oceanogr. (submitted).
Other references cited:
Batchelor, G.K., 1959, "Small scale variation of convected quantities like temperature in turbulent fluid", J. Fluid Mech. 5, 113133.
Itsweire, E., J. Koseff, D. Briggs & J. Ferziger, 1993, "Turbulence in stratified shear flows: Implications for interpreting shearinduced mixing in the ocean", J. Phys. Oceanogr. 23, 15081522.
Klaassen, G.P. & W.R. Peltier, 1991, "The influence of stratification on secondary instability in free shear layers", J. Fluid Mech. 227, 71106.
Kolmogorov, A.N., 1941, "The local structure of turbulence in incompressible viscous fluid at very large Reynolds number", reprinted in Proc. Roy. Soc. London 434, 913 (1991).
Kraichnan, R., 1968, "Small scale structure of a scalar field convected by turbulence", J. Fluid Mech. 11, 945953.
Hill, R.J., 1978, "Models of the scalar spectrum for turbulent advection", J. Fluid Mech. 88, 541562.
Moum, J.N., 1996, "Efficiency of mixing in the main thermocline", J. Geophys. Res., 101 (C5), 12,05712,069.
Oakey, N.S. 1982: Determination of the rate of dissipation of turbulent energy from simultaneous temperature and velocity shear measurements. J. Phys. Oceanogr. 12, 256271.
Sreenivasan, K.R., 1991, "On local isotropy of passive scalars in turbulent shear flows", Proc. Roy. Soc. London 434, 165182.
Wijesekera, H., & T. Dillon, 1997, "Shannon entropy as an indicator of age for turbulent overturns in the oceanic thermocline", J. Geophys. Res. 102, 32793291.