Direct simulations are used to study turbulence and mixing in Holmboe waves. Previous results showing that mixing in Holmboe waves is comparable to that found in the better-known Kelvin–Helmholtz (KH) billows are extended to cover a range of stratification levels. Mixing efficiency is discussed in detail, as are effective diffusivities of buoyancy and momentum. Entrainment rates are compared with results from laboratory experiments. The results suggest that the ratio of the thicknesses of the shear layer and the stratified layer is a key parameter controlling mixing. With that ratio held constant, KH billows mix more rapidly than do Holmboe waves. Among Holmboe waves, mixing increases with increasing density difference, despite the fact that the transition to turbulence is delayed or prevented entirely by the stratification. Results are summarized in parameterizations of the effective viscosity and diffusivity of Holmboe waves.
Because active scalars in seawater (primarily temperature and salinity) diffuse much more slowly than does momentum, there is a tendency to form thin layers of stable density stratification associated with shear on larger vertical length scales. In the simplest case where a single thin stratified layer lies at the center of a thicker shear layer such that shear and stratification are both symmetric about some central plane, the flow can support two classes of instability, both of which may lead to turbulence and mixing. These are the Kelvin–Helmholtz (KH) instability, which consists of a stationary train of billows focused at the central plane, and the Holmboe instability, which consists of a pair of oppositely propagating wave trains, focused above and below the central plane, that interfere to form a standing wave–like structure (Holmboe 1962; Smyth et al. 1988). In this paper, we examine mixing processes in direct numerical simulations (DNSs) of Holmboe instability. The study also includes a KH case for comparison.
Historically, KH instability has been studied more thoroughly because the conditions for its growth occur in both air and water and because its exponential growth rate tends to be significantly larger (Peltier and Caulfield 2003; Smyth et al. 2001, and references therein). Conditions for Holmboe instability are common in aqueous environments, however, and have been studied in exchange flows (Zhu and Lawrence 2001), river outflows (Yoshida et al. 1998), and salt wedge intrusions (Yonemitsu et al. 1996). Because stratification is concentrated in a thin layer, the instability exists even when the bulk Richardson number in the stratified layer is much larger than the usual critical value ¼ (Miles 1961).
Smyth and Winters (2003, hereinafter SW03) suggest that Holmboe instability may generate significant levels of turbulence and mixing, not despite the slow growth rate of the primary instability but because of it. The reasoning is that the preturbulent stage of instability evolution develops sharp scalar gradients that mix rapidly without a corresponding loss of kinetic energy to dissipation. If the preturbulent phase lasts longer in a slower-growing instability, then a slow growth rate favors strong mixing. Our objective here is to test this hypothesis using a set of simulations that cover the relevant parameter space more thoroughly than did SW03. We will show that the duration of the preturbulent phase is, indeed, a crucial factor governing the amount of mixing, but also that the connection between that duration and the growth rate of the primary instability is less simple that was assumed in SW03.
The DNS dataset and the methods used to create it are described in section 2. Section 3 gives a brief overview of the DNS results. In section 4, we discuss the energy budget for stratified turbulence and various useful measures of the efficiency of mixing. We discuss mixing in terms of entrainment rates and diffusivities for both velocity and buoyancy in section 5. Conclusions are summarized in section 6 and are discussed in the context of previous studies in section 7.
a. The mathematical model
Our mathematical model employs the Boussinesq equations for velocity u, density ρ, and buoyancy b in a nonrotating physical space measured by the Cartesian coordinates x, y, and z:
The variable ρo is a constant characteristic density. The buoyancy is related to the density by b = −g(ρ − ρo)/ρo. The gravitational acceleration is denoted g, and k is the vertical unit vector. Molecular viscous effects are represented by the Laplacian operator, with kinematic viscosity ν.
The augmented pressure field Π is specified implicitly by the incompressibility condition
and the buoyancy evolves in accordance with
in which κ is a molecular diffusivity for the stratifying scalar.
We assume periodicity in the horizontal dimensions,
in which f is any solution field and the periodicity intervals Lx and Ly are constants. At the upper and lower boundaries z = ±Lz/2, we impose an impermeability condition on the vertical velocity,
and stress-free, adiabatic conditions on the horizontal velocity components u and υ and on b,
These imply a condition on Π at the upper and lower boundaries:
b. Initial conditions
The model is initialized with a parallel flow in which shear and stratification are concentrated in horizontal layers centered on the plane z = 0,
The constants ho, Δu, and Δb represent the initial thickness of the shear layer and the associated changes in velocity and buoyancy. The thickness of the stratified layer is hb = h0R−1 (Fig. 1).
The problem is nondimensionalized using the length scale h0, the velocity scale Δu, and the time scale h0/Δu. These scales, together with the various physical parameters in the problem, combine to make a set of nondimensional parameters whose values specify the problem completely. The level of stratification relative to the shear is quantified by the initial central (or bulk) Richardson number Ri0, which is defined as the value of the gradient Richardson number at z = 0 and t = 0,
The nondimensional viscosity is 1/Re0, where
is the initial Reynolds number. (Both the central Richardson and Reynolds numbers increase over time as the shear layer spreads.) The relationship between viscosity and the scalar diffusivity is expressed by the Prandtl number:
In addition to the profiles described above, the initial conditions included a two-part perturbation designed to efficiently stimulate both primary and secondary instabilities. First, a disturbance proportional to the eigenfunction of the fastest-growing mode was added. The amplitude was chosen so that the maximum vertical displacement was ah0/2, where a = 0.4 or 0.2 (see Table 1). This amplitude is large enough to efficiently stimulate the primary mode, yet small enough that linear perturbation theory provides a reasonable approximation to the initial growth (SW03). Second, a random velocity field was added in order to excite three-dimensional motions. At each point in space, the three components of the velocity increment were chosen from a list of random numbers whose probability distribution was uniform between the limits ±cΔu/2, where c = 10−2 or 10−3 (Table 1). During the first time step, the random motions were automatically made solenoidal by the pressure gradient force.
The dimensions of the computational domain are set by the wavelength of the fastest growing mode of linear theory, which we describe in the next subsection.
c. Linear stability analysis
Figure 2 shows the stability characteristics of normal modes of the velocity and density profiles (9) and (10) versus R−1, the ratio of stratified layer thickness to shear layer thickness, and J = H0Δb/Δu2, which is proportional to the buoyancy change across the stratified layer. Results to follow are given in terms of the central Richardson number Ri0 = JR; hence that parameter is plotted for reference in Fig. 2d.
Figures 2a and 2b show the real and imaginary parts of the growth rate. Significant instability is found mainly in two regions—the stationary KH modes lying approximately in the wedge Ri0 < 0.25 (Fig. 2d) and the oscillatory Holmboe modes with Ri0 > 0.25, R−1 < 0.5. The wavenumber of the dominant instability (Fig. 2c) is nearly uniform over much of the KH regime, but varies considerably for Holmboe modes.
Bullets in Fig. 2 show the relationship between the present study and SW03. SW03 varied the central Richardson number in order to obtain KH and Holmboe instabilities, while keeping the net buoyancy change across the stratified layer fixed (the two circles with crosses). That experiment produced the surprising result that the Holmboe instability (R = 3, J = 0.15, Ri0 = 0.45) mixed more effectively than did the KH instability (R = 1, J = 0.15, Ri0 = 0.15), despite the fact that the latter had a much larger growth rate (Fig. 2a). Here, we hold R fixed while varying J. At one extreme, we set J = 0.05 (Ri0 = 0.15) to produce KH instability. In contrast, our most strongly stratified case had J = 0.33 (Ri0 = 1.0). This corresponds to the maximum growth rate for Holmboe modes with R = 3. Note that this mode is strongly oscillatory, with frequency exceeding growth rate by an order of magnitude. Between these extremes are three intermediate cases of Holmboe waves, one of which is identical with the Holmboe case investigated in SW03 (R = 3, J = 0.15, Ri0 = 0.45).
d. Parameter values
In accordance with the linear stability results summarized above, Holmboe-like waves are typically observed on stratified shear layers where the scale ratio R is significantly greater than unity and the minimum gradient Richardson number exceeds ¼. For example, Yoshida et al. (1998) observed cusped waves at a river inflow where R ∼ 4 (estimated from their Fig. 5) and J = 0.56.1 In a second observation, R was much smaller, and no instability was present. As usual, the parameter choices for our DNS runs represent a compromise between realism and computational practicality, spanning the Holmboe regime of Richardson number and scale ratio at a Reynolds number that is small relative to most geophysical mixing events (Table 1).
In addition to the five main cases in which Ri0 is varied, the dataset includes four auxiliary runs designed to test sensitivity to grid resolution (run 6), Reynolds number (run 7), and initial perturbation amplitude (runs 8 and 9). (The auxiliary runs are all variants on run 3, the Ri0 = 0.45 case.) For all cases, the initial Reynolds number Re0 is 1200. These choices, together with Pr = 9 (an appropriate value for water under typical geophysical conditions) and R = 3 completely specify the equations of motion and the initial conditions.
The domain length Lx is set equal to the wavelength of the fastest-growing mode. The domain width Ly is Lx/2, approximately 3 times the spanwise wavelength of the fastest-growing three-dimensional instability of KH billows in air as described by Klaassen and Peltier (1991). The domain height Lz is set equal to Lx/2, except in cases where extensive vertical motions were anticipated (Table 1).
e. Numerical methods
The numerical code is an extension of that described by Winters et al. (2004). It uses Fourier pseudospectral discretization in all three dimensions. Time stepping is via the third-order Adams–Bashforth operator, with time step determined by a Courant–Friedrichs–Lewy stability condition. Viscous and diffusive terms are integrated exactly. Buoyancy is resolved on a fine grid with spacing equal to one-half the spacing used to resolve the other fields. Aliasing errors are reduced by applying to both grids at every time step an isotropic filter with a cosine-bell shape that decreases gradually from amplitude 1 to 0.6 over the range 0.8–1 times the Nyquist wavenumber. For further discussions of the numerical method, see Winters et al. (2004) and Smyth et al. (2005).
Spatial resolution is isotropic and is designed so that, at peak turbulence levels, the smallest scales on which significant variations exist (a few times the Batchelor scale; see Moin and Mahesh 1998; Smyth et al. 2005) are well resolved. Here, the streamwise (x) wavelength of the primary instability is resolved using 128 grid nodes for velocity and pressure and 256 nodes for buoyancy (with corresponding grid lengths in the y and z directions to maintain isotropic resolution). An auxiliary test run was made with 160 nodes in x for velocity and pressure and 320 nodes for buoyancy (run 6 on Table 1). The resolution increase made no significant difference to the mixing statistics reported here, confirming that our spatial resolution is sufficient.
3. Overview of mixing events
The Holmboe wave consists of a pair of oppositely propagating wave trains that interfere to form a standing wave. At the state of minimum potential energy, the isopycnals are not flat but rather form sharp crests above and below the stratified layer. An example is shown in Fig. 3a, which is taken from the most strongly stratified case considered here, that for which Ri0 = 1. At the opposite point in the oscillation cycle, the state of maximum potential energy, the wave is approximately sinusoidal.
When stratification is reduced, the oppositely propagating component waves interact more strongly, and hence the Holmboe wave evolves complex, vortical structures characteristic of turbulence (Figs. 3c,d; also see SW03). At the other extreme, when Ri0 < 0.25, the component waves phase lock and wrap around each other to form KH billows. Each billow contains multiple layers of unstable density gradients (Fig. 3e) that quickly become turbulent (Fig. 3f).
Each of these cases involves irreversible mixing amplified to some degree by three-dimensional secondary circulations. The dynamics of the secondary circulations are described in detail in Smyth (2006). In each case, the disturbance eventually dies away, leaving behind a stable, parallel, stratified shear layer in which both shear and stratified layer thicknesses have been permanently increased through mixing.
4. The disturbance energy budget and the efficiency of mixing
a. Theoretical development
Sheared, stratified turbulence may be viewed as the intermediary through which kinetic energy from the mean flow is transferred into potential energy and internal energy. In this subsection, we describe the energy budget for a stratified shear flow in which there are no boundary fluxes (Fig. 4).
Kinetic energy is divided into two parts describing the mean flow and the disturbance. Similarly, potential energy is subdivided into a “background” part, defined as the minimum potential energy attainable via adiabatic rearrangement of fluid parcels (Winters et al. 1995), and an “available” part, which can be converted to background potential energy by mixing. These quantities evolve according to the following equations, as illustrated schematically in Fig. 4:
The mean flow kinetic energy is defined as
where the angle brackets denote an average over the coordinate directions indicated in the subscript. The mean flow is u(z, t) = 〈u〉xy. Also, 𝒮 represents the shear production
where u′ = u − u and w′ = w are disturbance velocities, and represents rate of mean flow kinetic energy conversion to internal energy via viscous dissipation:
Energy drawn from the mean flow by the shear production 𝒮 is converted into disturbance kinetic energy, defined as
The buoyancy flux is given by
where b′ = b − b and b(z, t) = 〈b〉xy. The viscous dissipation rate of disturbance kinetic energy is
is the disturbance strain tensor.
The total potential energy,
is subdivided as 𝒫t = 𝒫a + 𝒫b. The background potential energy is given by
where br(z, t) is a three-dimensionally reordered buoyancy field defined as in Winters et al. (1995). A positive buoyancy flux converts available potential energy into disturbance kinetic energy. The mixing rate has the expression
is the rate of potential energy increase due to molecular diffusion in the absence of fluid motions.
In many ocean mixing processes, one wishes to know the fraction of the available energy that is transferred irreversibly into the potential energy field. This need has led to various attempts to define the “efficiency” of mixing. Here we will discuss three definitions, all of which are useful in different ways. Perhaps the most obvious definition of mixing efficiency results from thinking of stratified turbulence as an ideal engine in the thermodynamic sense and computing its mechanical efficiency, that is, the ratio of work done to energy input. “Stratified turbulence” refers to the range of velocity and buoyancy fluctuations possible in a stratified fluid and is commonly thought of as a mixture of gravity waves and classical turbulence. We therefore identify the energy of the stratified turbulence as 𝒦′ + 𝒫a (shaded box on Fig. 4). The mechanical efficiency is then the ratio ei(t) = ℳ/𝒮, where the subscript i indicates that the ratio is evaluated instantaneously as a function of time. Rather than being an ideal engine, however, stratified turbulence is one of those engines that is capable of returning energy to its source, like an electric vehicle that can use energy to recharge its battery. When the shear production is negative, the engine is in this recharging state, and the ratio of work done to energy input is negative. The notion of efficiency thus becomes more complex and, perhaps, less useful than one would like.
A second approach is based on the ideas of McEwan (1983). We effectively abandon the mechanical analogy by defining efficiency as the ratio of work done to energy wasted (dissipated), that is, Γi(t) = ℳ/ε′. Here Γi is also useful in estimating the buoyancy flux from observational measurements of dissipation (e.g., Moum 1996) and for discriminating between shear-driven turbulence and double-diffusive mixing (e.g., St. Laurent and Schmitt 1999). Unlike ei, Γ never becomes negative. It can, however, become arbitrarily large when dissipation is small: Γi is related to ei by
The first factor on the right-hand side is nearly equal to ei when ei ≪ 1, as is usually the case. This familiar algebraic relation is supplemented by a factor describing the storage of energy by stratified turbulence, which will tend to increase Γi when turbulence is growing and decrease it when turbulence is decaying.
An elegant alternative approach was taken by Caulfield and Peltier (2000), who avoided partitioning the kinetic energy into a mean part and a disturbance part, effectively expanding the “engine” of stratified turbulence to include the mean flow. This engine has no energy source; it has, instead, a finite reservoir of energy that can lose (but never gain) energy via dissipation and mixing. Caulfield and Peltier defined the efficiency of this process as the ratio of mixing to total energy expended:
This approach is useful in the context of large-scale energy budgets (e.g., Munk and Wunsch 1998). Suppose, for example, that a wind event supplies the ocean current system with an increment of kinetic energy. Eventually, we expect that this energy increment will be shed due to instability and that a fraction of it will go into mixing the ocean: ℰi quantifies that fraction.
While instantaneous mixing efficiencies can give useful insights into the mechanics of mixing, variants that describe the cumulative mixing over a complete event may be of greater practical use:
Here Γc is simply related to ec because the storage term integrates to zero:
In high-Reynolds-number flows where is negligible, ℰc ≈ ec.
We now apply these ideas toward understanding the physics of mixing in our sequence of simulated mixing events. We begin with run 3, the Holmboe case with Ri0 = 0.45. This case is similar to the Holmboe wave considered in SW03, differing only in the details of the initial noise fields. After this, we will look at the more strongly stratified Holmboe wave (run 1) and the weakly stratified KH wave (run 5). We next discuss comparative results for all cases and close the section with a look at the effects of the initial noise amplitude.
b. A Holmboe wave in moderate stratification
For run 3 (Figs. 3c,d), the disturbance kinetic and available potential energies exhibit the oscillatory growth characteristic of Holmboe instability (Fig. 5a). Part of this oscillation is due to a wavelike exchange of energy between 𝒦′ and 𝒫a via the buoyancy flux. This aspect of the oscillation vanishes when the two energy reservoirs are added to create the total energy of the stratified turbulence. The remaining oscillation is due to a periodic exchange between the disturbance and mean flow kinetic energies. The background potential energy grows monotonically, and the most rapid growth is found just after the peak in 𝒦′ + 𝒫a, when the latter is declining rapidly.
Figure 5b shows the energy conversions that lead to mixing. These curves have been low-pass filtered to attenuate the strong oscillations in the shear production and hence in the storage rate. The dominant balance is between the shear production (dashed curve) and the dissipation (dotted curve). Both growth and decay of the turbulence and irreversible mixing are associated with imperfections in this balance. The imbalance has the character of a time lag, with production peaking earlier than dissipation.
The time lag between production and dissipation results from the downscale energy cascade that characterizes three-dimensional flow. Early in the flow evolution, the flow is dominated by large-scale structures that drive a vigorous flux of energy from the mean flow to the disturbance while dissipating energy at a relatively low rate. Some of the resulting excess energy goes into irreversible mixing, but the greater part is stored in the turbulence. Over time, turbulent energy cascades to smaller scales so that by the latter half of the mixing event, small-scale motions are more prevalent. This leads to an excess of dissipation over production (cf. Smyth and Moum 2000), which depletes the energy of the turbulence. As in the growth phase, a portion of this energy imbalance is converted irreversibly into background potential energy.
The evolution of Γi (thick curve on Fig. 5c) is typical of mixing events that originate from instabilities; Γi decreases from a relatively high value2 early in the mixing event to a value close to 1/Pr, the value for monochromatic gravity waves (Staquet and Bouruet-Aubertot 2001), after turbulence has subsided. Early in the event, the ratio ei is smaller, indicating that, with the oscillations of the Holmboe wave averaged out, somewhat less than 20% of the energy gained from the mean flow via shear production goes into raising the background potential energy 𝒫b. Equation (29) now suggests a way to understand variations in Γi. Early in the event, some of the energy from shear production that does not go into 𝒫b also does not go into dissipation, where it would reduce the value of Γi, but instead is stored as 𝒦′ + 𝒫a. In the decay phase, Γi remains close to the wave value 1/Pr, and it is ei that is elevated, as energy for ℳ comes not from shear production, which would reduce ei, but rather from the release of energy stored in 𝒦′ + 𝒫a. In the final stage of decay, shear production drops to nearly zero while weak mixing persists due to the continuing release of energy stored in the turbulence, causing ei to increase sharply.
The value of ℰi shows that, of the total energy lost to the mean flow at any given time, about 5%–10% goes into raising 𝒫b. This relatively small percentage indicates the importance of mean flow dissipation , which absorbs a significant fraction of mean flow kinetic energy in this low Reynolds number mixing event.
c. A Holmboe wave in strong stratification
The Holmboe wave at J = 0.33, Ri0 = 1 (run 1, Figs. 3a,b) is strongly oscillatory, as predicted by linear theory (Fig. 6). The oscillation dominates the evolution of 𝒦′ and 𝒫a to the extent that the two are indistinguishable on the scale of Fig. 6a. The oscillation is much smaller in the disturbance energy 𝒦′ + 𝒫a, indicating that the oscillatory energy transfer is almost entirely between 𝒦′ and 𝒫a via the buoyancy flux. The maximum value of 𝒦′ + 𝒫a is just above 3 × 10−3, very similar to the maximum energy attained in the previous case (Ri0 = 0.45). In contrast, the background potential energy grows to values much larger than those found at Ri0 = 0.45, a difference that is discussed in detail below.
Before plotting the energy transfers, we filter out the fast oscillations. This reveals a second oscillation of about 1/7 the frequency. The disturbance energy 𝒦′ + 𝒫a peaks at about t = 220. Before this, the shear production is around 4 × 10−5, which is quite similar to the range of values in the early part of run 3 (cf. Fig. 5b). Unlike run 3, however, this run does not exhibit a subsequent period of rapid disturbance growth associated with the turbulent breakdown of the Holmboe wave. Instead, the dissipation rate ɛ′ grows gradually until shortly after t = 220 when it becomes approximately equal to the shear production and the wave begins to decay. The mixing rate ℳ remains near 10−5 until after t = 300. This is in marked contrast to the previous case Ri0 = 0.45, in which ℳ decays rapidly after t = 180. Eventually, ℳ decreases to zero as in the previous case.
The evolution of Γi is similar to the previous case: it is undefined near t = 0, but settles down to a value near 0.4 for most of the growth phase. Then Γi decreases gradually to a value near 1/Pr. The ratio ei is just above ⅕ throughout most of the growth phase. It then increases to become greater than Γi for a brief period when the disturbance energy is decreasing, before decreasing to a value not much different from 1/Pr late in the run. Through most of the growth phase ℰi is near 10%, and then it decreases to zero.
The main result that we find in comparing this strongly stratified Holmboe wave with the more weakly stratified case discussed previously is that the strongly stratified wave accomplishes more mixing (cf. thick, solid curves in Figs. 5a and 6a). This strong mixing happens not despite, but because of, the fact that the wave never becomes turbulent. The growth phase, during which mixing is relatively efficient, is sustained for about 2 times as long because the wave fails to break, and the ultimate result is a greater irreversible increase in potential energy.
d. A KH billow
The KH instability at J = 0.05, Ri0 = 0.15 (run 5, Figs. 3e,f) reaches energy levels higher than in any of the Holmboe cases considered here (Fig. 7). The linear KH mode is not oscillatory, but at finite amplitude the disturbance exhibits slow nutations characteristic of an elliptical vortex in a shear flow. The shear production rate (Fig. 7b) reaches a maximum an order of magnitude larger than the fastest-growing Holmboe mode at Ri0 = 1, then oscillates briefly before settling down to a monotonic decrease. As in the Holmboe cases, the dissipation rate peaks later in the event, and this peak coincides with the release of energy stored during the earlier, nondissipative phase. As in all cases, a slow but steady conversion feeds energy irreversibly into the background potential energy reservoir.
Because energy transfers were not low-pass filtered, the ratio ei exhibits complex behavior, including singularities and negative values associated with sign changes in 𝒮. Once these oscillations cease, however, we see an increasing trend in ei that contrasts with a decreasing trend in Γi, as seen previously in the Holmboe wave cases. This again illustrates the effect of the storage and subsequent release of energy in a time-limited mixing event. Again Γi is large in the early stages, reaching peaks near 0.6; ℰi evolves much as Γi does, which is not surprising since ε′ ≫ for this case and hence ℰi ≈ Γi/ (1 + Γi). In general ℰi is larger than in the Holmboe waves, reaching a peak value of 0.3.
The cumulative mixing ∫t0 ℳ dt′ is greatest by far in the KH wave case (Fig. 8). Beside reaching the largest final value, it grows most rapidly in the early phases. This is counter to the result of SW03, who showed that a KH billow with R = 1 mixes less than a Holmboe wave with the same buoyancy difference Δb but R = 3. Here R = 3 in all cases and Δb is varied, so there is no discrepancy in the conclusions.
Among the Holmboe wave cases, the growth of ∫t0 ℳ dt′ begins earliest at Ri0 = 0.7, but the ultimate amount of mixing increases monotonically with Ri0, yielding its maximum value for the nearly laminar wave at Ri0 = 1.
The instantaneous mixing efficiency Γi of the Holmboe waves (Fig. 9) shows a clear dependence on Ri0. In the early phase (roughly t = 0–50), Γi is undefined because ε′ ≈ 0. Following this comes a period of elevated Γi, after which Γi decreases and ultimately approaches 1/Pr. This interval of high Γi corresponds to the preturbulent phase and is longest in the most strongly stratified case (where transition to turbulence does not occur). The length of the preturbulent phase decreases monotonically with decreasing Ri0, disappearing entirely in the most weakly stratified case Ri0 = 0.3, where the wave breaks rapidly.
The cumulative mixing efficiency, Γc = ∫tf0 ℳ dt/∫tf0 ε′dt, varied considerably among the different cases (Fig. 10). The KH case showed the highest value. Among the Holmboe cases, the highest Γc occurred at Ri0 = 1. The open circle on Fig. 10 indicates Γc for the high-Reynolds-number case with Re0 = 2400, Ri0 = 0.45. The similarity between this case and its counterpart with Re0 = 1200 suggests that dependence on the Reynolds number is weak.
The net fraction of mean flow energy that went into irreversible potential energy increase ℰc is generally smaller than Γc but otherwise behaves similarly: it is greatest for the weakly stratified KH wave but increases monotonically with stratification for the Holmboe wave cases. In contrast to Γc, ℰc is sensitive to the Reynolds number. In the single high-Reynolds-number case, ℰc increased significantly owing to the reduced role of mean flow dissipation.
f. Dependence on the initial disturbance amplitude
The stages in a mixing event depend to some degree on the characteristics of the initial disturbance. The cases shown in Fig. 11 are identical except that the initial random noise field was reduced by an order of magnitude in the latter. The results are identical until about t = 100, at which point three-dimensional motions become significant in the run with larger initial noise. While the mixing rate ℳ is initially unaffected, increased dissipation causes Γi to decrease. More important, the increased dissipation causes the turbulence to decay so that ℳ begins to decrease at about t = 170. In contrast, the event with lower initial noise does not develop strong dissipation until t ∼ 200. The cycle of increased dissipation, reduced Γi, and finally reduced mixing is the same except for being displaced in time. This apparently minor difference has a major impact on the net amount of mixing accomplished. The net mixing was increased by nearly 20% (Γc increased from 0.17 to 0.20) because the highly efficient, two-dimensional stage lasted longer when the initial noise was reduced.
Next, we investigate dependence on the initial amplitude of the primary instability. Runs 8 and 9, shown in Fig. 12, are identical except that the initial eigenfunction amplitude was halved for run 9. This has the effect of delaying the growth of the two-dimensional Holmboe wave. As expected, mixing is weaker at early times (t < 100). Dissipation is also weaker during this time, though. For t = 50–150 (after the initial phase in which it is undefined), Γi is essentially equal for the two cases. But, the rapidly growing two-dimensional instability in run 8 becomes turbulent and collapses around t = 200, whereas the wave with the smaller initial amplitude does not break until t = 250. Once again, the result of a weaker initial perturbation is that the highly efficient two-dimensional mixing stage is prolonged. In this case, the net amount of mixing is increased by about 10% (Γc increased from 0.20 to 0.22).
5. Entrainment and diffusivity
In all cases, both the shear layer and the stratified layer thickened over time due to the combination of molecular diffusion, turbulent mixing, and wave motions (Fig. 13). Layer thicknesses were calculated as
where f(z,t) represents the mean profile of either velocity or buoyancy and Δf is its net change. Both hu and hb exhibited an initial period of roughly linear growth modulated by oscillations. This growth rate was used to construct a dimensionless entrainment rate
For the Holmboe cases, Eu exceeds Eb (Fig. 14). The difference is slight at Ri0 = 0.3 but becomes more pronounced with increasing stratification. In the KH case, Eb was slightly larger than Eu, and both were larger than the largest values found in the Holmboe cases. In the Holmboe cases, the entrainment rate describing the thickening of the shear layer shows no distinct dependence on stratification, being within 10% of 0.005 in all cases. In contrast, the thickening of the stratified layer is slower and more dependent on stratification. It is clear that entrainment is strongly decreased by stable stratification in all cases, as the comparable entrainment rate for an unstratified shear layer is O(10−1) (e.g., Rogers and Moser 1994).
The diffusivity Kirrb describing irreversible scalar mixing is computed instantaneously based on the Winters and D’Asaro (1996) formalism. At each time step, the buoyancy field is resorted to form the state of minimum potential energy 𝒫b as described in section 4a. The evolution of 𝒫b is then used to define a spatially uniform (but temporally varying) diffusivity. A detailed description of the method may be found in Smyth et al. (2005).
In each case, this diffusivity rises to a maximum then decays as the disturbance subsides (Fig. 15). The peak diffusivity decreases with increasing stratification, and is much larger in the KH case than in the Holmboe cases. This dependence on stratification is partially countered by the fact that the interval during which Kirrb is significantly greater than the molecular level lasts longer in the more strongly stratified cases. Also, the fluxes may be larger in the more strongly stratified cases because the gradients are larger.
Unfortunately, we know of no way to compute an instantaneous momentum diffusivity that describes irreversible mixing only. Instead, we compute characteristic diffusivities for each event based on the kinetic energy budget. The eddy viscosity is defined via the shear production rate:
where the time integrals span the entire mixing event. Similarly, we compute a net mass diffusivity based on the buoyancy flux:
As is true of the entrainment rate, the net diffusivities for the KH case are much larger than those for the Holmboe cases (Fig. 16). The diffusivities are approximately independent of Ri0 for the Holmboe cases. Although there is a hint of systematic decrease in Kbuoy, the present data are insufficient to assess its significance. The disturbance Prandtl number Kmom/Kbuoy is near unity for the KH case and is about 3 for all Holmboe cases.
The net diffusivities for the Holmboe cases are low (i.e., comparable to molecular values) as a result of the low Reynolds numbers attained in these DNS experiments. Based on scaling arguments, one would expect the eddy diffusivities to increase in direct proportion to the Reynolds number. In an auxiliary simulation with Re = 2400 (run 7, plus sign and asterisk on Fig. 16), the nondimensional diffusivities essentially equaled their values in the corresponding Re = 1200 run (i.e., after doubling Re, the diffusivities increased by only a few tens of percent). This means that, in dimensional terms, the diffusivities scale with Reynolds number, as expected. Based on the mean values for the Holmboe cases, we suggest that Holmboe wave diffusivities may be predicted as Kmom = 3.6 × 10−4h0Δu and Kbuoy = 1.2 × 10−4h0Δu. This prediction is of limited practical value, as the initial layer depth is difficult to specify in practice. In cases studied here, however, the shear layer thickened by a factor of 2.0 ± 0.5, including the run at higher Reynolds number.
The shear layer thickness, observed at any given time, would therefore be between 1 and 2.0 ± 0.5 times h0. One could assume that h0 is the observed shear layer thickness h divided by 1.5, and is correct to within a factor of 2. Based on this, we propose the parameterization
where h and Δu are the instantaneous thickness of the shear layer and velocity change across it. This parameterization should be used with caution when Ri0 ≫ 1 or Re0 ≫ 103.
We have investigated various aspects of mixing in a set of simulations consisting of four Holmboe waves at various levels of stratification and a single KH wave, along with auxiliary Holmboe wave simulations in which the Reynolds number and the initial noise levels were varied. The main conclusions are as follows.
Both KH billows and Holmboe waves follow a life cycle that is probably characteristic of time-limited, shear-driven mixing events in general. First, shear production causes a transfer of energy from the mean flow to the disturbance energy reservoir, 𝒦′ + 𝒫a. When the instability grows to asufficiently large amplitude, secondary circulations grow and, in most cases, trigger the onset of turbulence and the breakdown of the primary instability. Beyond this point, shear production decreases and energy stored in the disturbance decays. Disturbance energy cascades to small scales over time so that its decrease in the later phase is mostly due to viscous dissipation. At each stage of the event, a small fraction of the energy transfer goes irreversibly into raising the background potential energy.
We have examined three ratios that quantify the efficiency of mixing in different ways. Two of these become large at some point in the event—Γi = ℳ/ε′ is large in the growth phase mainly because its denominator is small. The reverse is true of ei = ℳ/𝒮, which is large later in the decay phase because its denominator is small: ℰi behaves very differently, increasing to a peak value when mixing is strongest, then decaying again. These changes are associated with storage and subsequent release of energy by the disturbance. Here ℰi is relatively small, peaking at 0.20 for the KH billow and less than 0.10 for the Holmboe waves.
When stratification is strong, the transition to turbulence occurs late or not at all. As a result, the two-dimensional Holmboe wave train, along with its characteristic large values of Γi, persists for a relatively long time.
The high values of Γi in the preturbulent phase do not indicate especially strong mixing, as those values are high primarily because ε′ is small. The latter fact, however, indicates that turbulent breakdown is not occurring. The longer this state persists, the greater the net mixing. As a result, the net mixing due to Holmboe waves is greatest when the laminar phase lasts longest, that is, in the strongly stratified cases.
Similarly, the cumulative mixing efficiency of Holmboe waves, by any measure, is greatest in the most strongly stratified cases. (If we continued the sequence of simulations to still higher Ri0, it is possible that this trend would eventually reverse.)
The present results pertain to relatively low Reynolds number events; however, Γi appears to be quite insensitive to Re0. Smyth et al. (2001) have demonstrated similar Reynolds number independence of Γi in KH billows. In contrast, ℰi involves mean flow dissipation and is therefore distinctly larger in the higher Reynolds number case.
Entrainment rates (scaled by Δu) for the KH billow are (2–3) × 10−2. For the Holmboe waves, momentum entrainment rates are (4–5) × 10−3 for all Ri0, while the scalar entrainment rate is smaller and decreases with increasing Ri0.
The instantaneous scalar diffusivity of the KH billow reaches a peak value two orders of magnitude above the molecular value. Among the Holmboe waves, the maximum value is smaller and decreases at stronger stratification, but the interval over which large values persist increases.
Net diffusivites are large for the KH case, small and fairly uniform among the Holmboe cases. Dimensional diffusivities appear to be proportional to the Reynolds number. Our results are summarized in (39), which may be used to parameterize eddy viscosity and diffusivity in Holmboe waves in this parameter range. A salient feature of (39) is that the effective Prandtl number for Holmboe waves is approximately 3. The data hint that the scalar diffusivity decreases at higher Ri0; if so, the effective Prandtl number would increase further in that parameter regime. Further experiments are needed to test this possibility.
Mixing events that begin in an environment relatively free of preexisting turbulence mix more effectively. This is because the onset of turbulence (which mixes vigorously but inefficiently and thus annihilates the wave) is delayed.
In studies of turbulent Kelvin–Helmholtz billows, Smyth et al. (2001) showed that Γi is highest early in the mixing event, when the flow is not yet turbulent. This result, together with the alignment statistics of KH billows (Smyth 1999), suggested the following interpretation. In the preturbulent phase, the strain field evolves slowly enough that scalar gradients are able to approach the optimal alignment for compression, leading to rapid diffusion, even though the velocity field is only weakly dissipative. Thus, any mixing event that involves a lengthy preturbulent phase is likely to achieve significant mixing with a relatively small expenditure of kinetic energy.
This view was challenged by the results of SW03, who found elevated mixing efficiency in a preturbulent Holmboe wave, despite the fact that the preturbulent strain field was oscillatory. Values of Γi in preturbulent Holmboe waves are generally not as large as those found in KH billows, however [Figs. 5c, 6c, 7c; Smyth et al. (2001)]. This suggests that preturbulent Holmboe waves are intermediate between preturbulent KH billows and fully developed turbulence, both in terms of persistent strain and in terms of efficient mixing. Thus, the hypothesis that high mixing efficiency in preturbulent flows is caused by persistent strain remains viable.
SW03’s suggestion that a low linear growth rate may favor effective mixing is not supported by the present results; in fact, the overall level of mixing correlates positively with the linear growth rate of the primary instability. Our results do, however, support SW03’s suggestion that a long preturbulent phase favors effective mixing. SW03 assumed that a slower-growing primary instability would necessarily take longer to become turbulent. We have seen that this is not so; the strongly stratified Holmboe cases become turbulent later than the weakly stratified cases despite their higher growth rates. The ultimate effect of this delayed transition to turbulence is increased net mixing. In the extreme case Ri0 = 1, the primary growth rate is a maximum but the transition to turbulence is suppressed entirely. The result is a very long preturbulent phase, and the greatest irreversible potential energy gain of all the Holmboe waves.
To compare mixing in KH and Holmboe waves, SW03 held the buoyancy difference across the stratified layer constant and varied its thickness; specifically, they set R = 3 to obtain a Holmboe wave and R = 1 for a KH wave. They found that levels of potential energy gain were similar in the two cases; in fact, the Holmboe wave generated a somewhat greater level of mixing. Taken by itself, this result would appear to suggest that the intrinsic mechanisms of the KH and Holmboe waves are similarly effective at generating mixing. The laboratory experiments of Strang and Fernando (2001) contradict this conclusion, showing that mixing is significantly stronger, by various measures, in the KH case.
The present results resolve this discrepancy. In contrast to SW03, we have held R constant at 3 while varying the buoyancy change across the stratified layer. (This distinction is illustrated in Fig. 2.) The present approach is a better approximation to that taken in laboratory experiments (e.g., Strang and Fernando 2001). The KH billow with R = 3 turns out to mix much more effectively than Holmboe waves with the same R. Taken together, these results show that the scale ratio R is a (and perhaps the) crucial factor in determining mixing efficiency in stratified shear layers. This is not hard to understand when one notes that, in flows with large R, scalar gradients are concentrated in the regions of strongest shear.
The effective Prandtl number for Holmboe waves is approximately 3 in the parameter range 0.45 < Ri0 ≤ 1, and may trend toward higher values at higher Ri0. This property is potentially important for the longevity of Holmboe waves (and their associated mixing) since it tends to maintain the sharpness of the density interface relative to the shear layer. As noted by Alexakis (2005), the regime Ri0 ≫ 1 may be complicated by higher harmonics of the fundamental Holmboe instability.
We gratefully acknowledge helpful discussions with Joe Fernando and production support from Rita Brown. The work was funded by the National Science Foundation under Grant OCE0221057.
Corresponding author address: Dr. William Smyth, College of Oceanic and Atmospheric Sciences, Oregon State University, Corvallis, OR 97330. Email: firstname.lastname@example.org
SW03 found Γi values that were 2 times as large as those reported here for the same case. This was due to a coding error that resulted in ε′ being too small by a factor of 2; hence, Γi is too large by the same factor. The error appears only in Figs. 12c–f of SW03, and it does not affect the conclusions.