Abstract

Motivated by the tendency of high-Prandtl-number fluids to form sharp density interfaces, the authors investigate the evolution of Holmboe waves in a stratified shear flow through direct numerical simulation. Like their better-known cousins, Kelvin–Helmholtz waves, Holmboe waves lead the flow to a turbulent state in which rapid irreversible mixing takes place. In both cases, significant mixing also takes place prior to the transition to turbulence. Although Holmboe waves grow more slowly than Kelvin–Helmholtz waves, the net amount of mixing is comparable. It is concluded that Holmboe instability represents a potentially important mechanism for mixing in the ocean.

1. Introduction

Turbulence in the ocean is often governed by a competition between vertical shear, which promotes instability, and statically stable density stratification, which acts to stabilize the flow. A simple, oceanographically relevant model for this physical regime is the stratified shear layer. In this model, flow evolves from initial conditions in which two homogeneous water masses of different densities are separated by a horizontal, sheared interface (as may occur, for example, in a river outflow). The lower layer is presumed to have the higher density so that the system is statically stable. If the shear across the interface is strong enough to overcome the stabilizing effect of the stratification, instability leads to the development of wavelike structures that may break and generate turbulence. In the absence of forcing, turbulence eventually decays and the mean flow relaxes to a stable, parallel state. Mixing and potential energy gain resulting from such events represent an important facet of ocean dynamics that is only partially understood.

Previous studies of turbulence evolving in stratified shear layers have usually assumed that the vertical changes in velocity and density change occur on similar length scales. In this case, sufficiently strong shear results in the well-known Kelvin–Helmholtz (KH) instability, which leads to the growth of periodic arrays of billows (e.g., Thorpe 1987). However, the low molecular diffusivities of heat and (especially) salt in the ocean suggest that this assumption may not always be valid; that is, there may be a tendency for density to change more abruptly than velocity (Fig. 1). If R, the ratio of length scales over which velocity and density change, is greater than about 2.4, strong stratification no longer acts to stabilize the flow but rather causes Kelvin–Helmholtz instability to be supplanted by Holmboe instability, whose finite-amplitude expression is an oscillatory, standing wavelike structure (e.g., Smyth et al. 1988; Smyth and Peltier 1989). Holmboe waves are difficult to observe unambiguously due to their oscillatory nature. However, the conditions for their growth are known to occur in exchange flows (Zhu and Lawrence 2001), salt wedge intrusions (Yonemitsu et al. 1996), and river outflows (Yoshida et al. 1998). Low Reynolds number mixing events such as those that dominate in the main thermocline (Smyth et al. 2001) are also expected to be influenced by Holmboe-like dynamics.

Fig. 1.

Profiles of streamwise (a) velocity and (b) density that define a stratified shear layer [(11), (12)]. In this case, the scale ratio R = 3.

Fig. 1.

Profiles of streamwise (a) velocity and (b) density that define a stratified shear layer [(11), (12)]. In this case, the scale ratio R = 3.

Our purpose is therefore to investigate the nature of turbulence that results from Holmboe waves in a strongly stratified shear layer with a sharp density interface. This flow geometry has largely been neglected as a potential source of turbulence in the ocean (Thorpe 1987; Gregg 1987), in part because the exponential growth rate of Holmboe instability is generally much smaller than that of Kelvin–Helmholtz instability. Here, we will demonstrate that this neglect is unjustified, that in fact Holmboe instability can generate strong turbulence and can lead to mixing comparable with Kelvin–Helmholtz instability.

Kelvin–Helmholtz billows have been studied extensively, both in the laboratory (e.g., Thorpe 1973, 1985; Koop and Browand 1979) and theoretically. Theoretical studies have involved linear stability analyses of the primary instability (e.g., Hazel 1972), nonlinear simulations of the primary instability in two dimensions (e.g., Patnaik et al. 1976; Klaassen and Peltier 1985a), and secondary stability analyses of the finite-amplitude, two-dimensional billow (Klaassen and Peltier 1985b, 1991; Caulfield and Peltier 2000). The transition process has been simulated in three dimensions by Caulfield and Peltier (1994) and Cortesi et al. (1998). The full KH life cycle, including the growth and decay of turbulence, has been simulated more recently as computer capacities have become equal to the task (e.g., Scinocca 1995; Smyth 1999; Cortesi et al. 1999; Smyth and Moum 2000b; Caulfield and Peltier 2000; Staquet 2000; Smyth et al. 2001). Observations of KH-like billows in the ocean have been documented in several studies (Farmer and Smith 1980; Hebert et al. 1992; Seim and Gregg 1994; DeSilva et al. 1996; Farmer and Armi 1999; Moum et al. 2002, manuscript submitted to J. Phys. Oceanogr.). Smyth et al. (2001) have shown that turbulence developing in direct numerical simulations of KH billows is statistically consistent with turbulent patches in the thermocline as revealed by in situ microstructure measurements.

Studies of the Holmboe wave have been much less numerous. The primary instability was first described theoretically by Holmboe (1962) and was investigated subsequently by Hazel (1972) and Smyth and Peltier (1989, 1990). A mechanistic description of the primary instability has been provided by Baines and Mitsudera (1994) (also see Caulfield 1994). Laboratory investigations have been carried out by Browand and Winant (1973), Lawrence et al. (1991), Pouliquen et al. (1994), Lawrence et al. (1998), Zhu and Lawrence (2001), Strang and Fernando (2001), and Hogg and Ivey (2002, hereinafter HI). The nonlinear growth of the Holmboe wave has been simulated in two dimensions by Smyth et al. (1988), Smyth and Peltier (1991), and Lawrence et al. (1998). Secondary stability analysis was carried out in a small range of parameter values straddling the KH–Holmboe transition by Smyth and Peltier (1991). The three-dimensional development of the Holmboe wave is computed here for the first time.

Our objective is to perform a direct comparison of turbulence evolution in Holmboe and KH waves. To do this, we consider two very similar stratified shear layers. The only difference between the two cases is the relative scales over which density and velocity change. In the first case, velocity and density vary over similar vertical distances and turbulence evolves through KH instability. The second case is identical except that the density change is concentrated in a thinner layer (as in Fig. 1) so that the flow exhibits Holmboe instability. These two flows are investigated by means of direct numerical simulation (DNS). In each case, we characterize turbulence and mixing in terms of irreversible potential energy gain and mixing efficiency. The simulated Holmboe waves, like KH waves, exhibit high mixing efficiency prior to the transition to turbulence (Winters et al. 1995; Smyth and Moum 2000b; Caulfield and Peltier 2000; Staquet 2000; Smyth et al. 2001). Because the Holmboe waves grow slowly, this highly efficient preturbulent mixing phase lasts for a relatively long time. In both KH and Holmboe cases, the subsequent turbulent phase is characterized by mixing efficiency near the standard value of 0.2. The net potential energy increase is actually greater in the Holmboe case.

In section 2, we describe the mathematical model upon which our calculations are based. Linear stability analyses are used to guide the choice of parameter values for the nonlinear simulations. Our main results are given in section 3. We describe the temporal evolution and spatial structure of turbulent KH and Holmboe waves, and seek to relate the transition to turbulence to existing understanding of the secondary instabilities of the two-dimensional waves (sections 3a and 3b). We then compare the two flows in terms of irreversible potential energy gain (section 3c). The main conclusions are summarized in section 4, and directions for future work are outlined in section 5.

2. Methodology

a. The mathematical model

Our mathematical model employs the Boussinesq equations for velocity u, density ρ, and pressure p in a nonrotating physical space measured by the Cartesian coordinates x, y, and z:

 
formula

Subscripts following commas denote partial differentiation. The variable ρ0 is a constant characteristic density, here set equal to 1027 kg m−3. The gravitational acceleration g has the value 9.81 m s−2, and k is the vertical unit vector. Viscous effects are represented by the usual Laplacian operator with kinematic viscosity ν = 0.995 × 10−6 m2 s−1.

The augmented pressure field Π is specified implicitly by the incompressibility condition

 
∇·u = 0,
(3)

and the density evolves in accordance with

 
ρ,t = −∇·ρu + κ2ρ,
(4)

in which κ is a molecular diffusivity for density. Here, we use the value 1.06 × 10−7 m2 s−1 for κ in order to model thermal stratification in the ocean.

We assume periodicity in the horizontal dimensions:

 
f(x + Lx, y, z) = f(x, y + Ly, z) = f(x, y, z),
(5)

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:

 
formula

and stress-free, adiabatic conditions on the horizontal velocity components u and υ and on ρ:

 
formula

These imply a condition on Π at the upper and lower boundaries:

 
formula

b. Numerical solution methods

Spatial derivatives are computed using full Fourier transforms in both horizontal directions and half-range sine and cosine transforms in the vertical, as required by the boundary conditions. The time evolution of the viscous and diffusive terms in (1) and (4) is evaluated exactly in Fourier space. (This is made possible by our choice of Fourier discretization on equally spaced collocation points.) The remaining terms are stepped forward in time using a third-order Adams–Bashforth method.

Because the dynamics of interest here depend on the low diffusivity of seawater, it is critical that we maintain adequate spatial resolution in the scalar field. In unstratified flows, the resolution requirement is determined by the Kolmogorov scale LK = (ν/ε)1/4, where ε is the kinetic energy dissipation rate. Grid spacing of 4LK is generally sufficient: Moin and Mahesh (1998) list several examples of successful simulations with grid spacing ranging from 3.75LK to 14.3LK. (Flows near solid boundaries are exceptions, as finer resolution is needed in the wall-normal direction.) In stratified flow with Prandtl number (Pr ≡ ν/κ) in excess of unity, the density field requires especially fine resolution, and the appropriate scale is the Batchelor scale, LB = LK/Pr1/2 (Batchelor 1959). Previous studies of stratified shear layers with high Prandtl number have used grid spacing ≤2.5LB (e.g., Smyth 1999; Smyth and Moum 2000a,b; Smyth et al. 2001). For the present work, we have been especially cautious in this regard. The minimum value of LB, which occurs at the center of the domain just after the transition to turbulence, is about one half the grid spacing; that is, the grid spacing never exceeds 2LB. (In addition, we use fully spectral discretizations as opposed to the spectral–finite difference methods used in the studies cited above.)

Our grid increments were set equal at

 
Δx = Δy = Δz = Lx/255
(9)

and the array sizes were set to

 
{Nx, Ny, Nz} = {256, 128, 128}.
(10)

Later in this paper (section 3b), spatial resolution will be checked directly using wavenumber spectra of the density gradient field.

The code is parallelized using Message Passing Interface (MPI) directives. The simulations described here were run using 32 processors on a Cray T3E.

c. Initial conditions

The model is initialized with a parallel flow in which shear and stratification are concentrated in a horizontal layer surrounding the plane z = 0:

 
formula

The constants h0, Δu, and Δρ represent the initial thickness of the shear layer and the associated changes in velocity and density; R is the ratio of shear layer thickness to stratified layer thickness (Fig. 1).

In order to obtain a turbulent flow efficiently, we add to the initial mean profiles a perturbation field designed to stimulate both two-dimensional primary and three-dimensional secondary instabilities. The horizontal velocity components are prescribed explicitly as described below; the vertical component is then obtained through (3). Primary instability is stimulated by a velocity perturbation whose streamwise component is

 
formula

Because this initial perturbation is weak enough to obey linear physics, its precise form has little effect on the statistical quantities of interest here. Random perturbations are added to the velocity and density fields to seed secondary instabilities. The density fluctuations have the form

 
formula

where rρ is a random deviate distributed uniformly between −1 and 1. The vertical dependence is chosen to enforce the limits |ρρ0| ≤Δρ/2. The random horizontal velocity perturbations have the form

 
formula

where −1 < ru < 1 and similarly for υrand.

d. Parameter values

The constants h0 and Δu and Δρ and R that define the initial mean flow can be combined with the fluid parameters ν and κ and the geophysical parameter g to form four dimensionless groups whose values determine the stability of the flow at t = 0. In addition to R, these are

 
formula

The initial macroscale Reynolds number Re0 expresses the relative importance of viscous effects. In the present simulations, Re0 is set to 1200, large enough that the initial instability is nearly inviscid. The Prandtl number Pr is set to 9, the appropriate value for thermal stratification in water at 12°C. The centerline Richardson number Ri0 is equal to the initial gradient Richardson number Ri(z) = −gρ̃,z/ρ0ũ2,z evaluated at z = 0, and it thus quantifies the relative importance of shear and stratification within the stratified layer. As the flow evolves, both the centerline Richardson and Reynolds numbers increase in proportion to the increasing thickness of the shear layer (Smyth and Moum 2000b).

There are two common choices for the length, velocity, and density scales used to describe stratified shear layers; therefore, care must be taken when comparing different studies. Theoretical studies (e.g., Smyth and Peltier 1989; Cortesi et al. 1998) have often used h0/2, Δu/2, and Δρ/2, because this choice simplifies some of the mathematical expressions. The choice of h0, Δu, and Δρ has been made more commonly in the observational and experimental literature (e.g., Seim and Gregg 1994; Zhu and Lawrence 2001; HI). Because our DNS work is done in close coordination with observations (e.g., Smyth et al. 2001), we prefer the latter convention.

The potential for instability may be guaged a priori by applying the Miles–Howard criterion (Miles 1961; Howard 1961), which states that instability is only possible if Ri(z) < 1/4 for some z. (This condition for instability is necessary but not sufficient.) The vertical dependence of Ri depends on the value of the scale ratio R (Fig. 2). When R < 2, Ri(z) has a single minimum, Ri(0) = Ri0, and increases without bound at large |z|. Ri0 < 1/4 is therefore a necessary condition for instability in this class of flows. When 2 < R < 2, z = 0 is a local maximum of Ri(z), and is flanked by minima located in the upper and lower halves of the shear layer. No unstable modes have been found in this regime to our knowledge, although their existence is not prohibited by the Miles–Howard criterion. When R > 2, Ri(z) vanishes far above and below the shear layer, and instability is therefore allowed for any Ri0. In fact, instability at Ri0 > 1/4 requires R > 2.4 (Smyth and Peltier 1989).

Fig. 2.

Profiles of the gradient Richardson number for (11) and (12) with Ri0/RJ = 0.15. Labels indicate the value of R for each curve. The thicker curves represent the cases to be investigated through DNS. The vertical dashed line indicates Ri = 1/4

Fig. 2.

Profiles of the gradient Richardson number for (11) and (12) with Ri0/RJ = 0.15. Labels indicate the value of R for each curve. The thicker curves represent the cases to be investigated through DNS. The vertical dashed line indicates Ri = 1/4

To motivate the remaining parameter choices, we now explicitly compute the linear stability characteristics of the stratified shear flow described by (11)–(16) with Re0 = 1200 and Pr = 9. Equations (1)–(4) are linearized about the parallel state (11) and (12). Aside from small effects due to viscosity and diffusion, (11) and (12) represent a steady solution, and so perturbations can be written in the normal mode form:

 
ϕ(x, y, z, t) = ϕ̂(z)eσt+i(αx+βy).
(17)

Here, ϕ may represent any perturbation field, α and β are streamwise and spanwise wavenumbers nondimensionalized by 2/h0, and σ is a (possibly complex) exponential growth rate nondimensionalized by Δu/h0. An extension of Squire's theorem (e.g., Smyth and Peltier 1990) shows that the dominant unstable modes of this flow are two-dimensional in the region of interest to us, and so we assume that β = 0 for these calculations. Discretizing the vertical dependence (using third-order compact derivatives) yields a matrix eigenvalue problem whose solution furnishes the growth rates and structure functions of all normal modes.

We now consider a class of flows that are identical except for the value of the scale ratio R. The parameters Δρ, Δu, and h0 are fixed so that

 
formula

The ratio Ri0/R represents a bulk Richardson number and is often denoted “J” (e.g., HI). The value 0.15 was chosen for convenience in the nonlinear simulations to follow. In this class of flows, the stratification profile varies from weak stratification spread over a thick layer (at small R and Ri0) to strong stratification concentrated in a thin layer (at large R and Ri0), but the net density change does not vary.

Figure 3 shows the growth rate (the real part of σ) as a function of α for a range of R values. Below R = 1.6 (Ri0 = 0.25), a band of unstable wavenumbers is evident. These are stationary modes (real phase speed is zero) and are an extension of the KH instability. The flow is stable to perturbations of the form (17) for 1.6 < R < 2.4. For R > 2.4, we find a domain of instability that does not exist in the R = 1 case. These are the oscillatory modes first described by Holmboe (1962). They occur in pairs having equal growth rate and equal but oppositely directed phase speed. At finite amplitude, these modes produce the standing-wave-like structures now known as Holmboe waves.

Fig. 3.

Nondimensional growth rate vs streamwise wavenumber and centerline Richardson number for linear, normal mode instabilities of a stratified shear layer. Ri0/RJ = 0.15, Re0 = 1200, and Pr = 9. The domain depth Lz is 4h0. Waviness in the stability boundary is due to the difficulty of resolving modes with vanishingly small growth rates. The labels “1” and “2” indicate the two cases to be investigated via DNS

Fig. 3.

Nondimensional growth rate vs streamwise wavenumber and centerline Richardson number for linear, normal mode instabilities of a stratified shear layer. Ri0/RJ = 0.15, Re0 = 1200, and Pr = 9. The domain depth Lz is 4h0. Waviness in the stability boundary is due to the difficulty of resolving modes with vanishingly small growth rates. The labels “1” and “2” indicate the two cases to be investigated via DNS

In this paper, we will examine direct numerical simulations of two flows chosen to represent the KH and Holmboe regimes (denoted as “1” and “2,” respectively, on Fig. 3). The only difference between the two flows is that the density changes over the same length scale as the velocity (i.e., R = 1) in the KH case but changes over a shorter length (R = 3) in the Holmboe case. Because of this, the density gradient maximum at z = 0 is stronger by a factor of 3 in the Holmboe case, and Ri0 is correspondingly larger. The nondimensional exponential growth rates of the KH and Holmboe modes are 0.069 and 0.022, respectively. Thus, by compressing the stratified layer, we have reduced the growth rate by more than a factor of 3. One might reasonably expect that this stabilization of the original flow would reduce the strength of the resulting turbulence, or even prevent its emergence entirely. Our DNS experiments will test this prediction.

The domain length Lx was chosen so as to accommodate one wavelength of the primary instability since pairing instability is not present at the high levels of stratification used here. [The pairing instability of KH billows is discussed by Klaassen and Peltier (1989). The fact that Holmboe waves do not pair in this region of parameter space was confirmed using an auxiliary sequence of two-dimensional simulations.] The length scale h0 was set to 0.1795 m. The (dimensional) wavelength was chosen as 2π/α × h0/2 with α = 0.35 (cf. Fig. 3), so that Lx = 1.6116 m. Here Lz was chosen as Lx/2, sufficient to prevent the upper and lower boundaries from influencing the flow evolution until after the turbulent mixing phase was complete. The domain width Ly was chosen to be Lx/2. This aspect ratio is justified in light of the small spanwise wavelength of the dominant secondary instabilities that lead the flow to a turbulent state (Klaassen and Peltier 1991; Smyth and Peltier 1991; Caulfield and Peltier 2000). Sensitivity experiments in which Ly was varied yielded no difference in the turbulence statistics of interest here, although scale selectivity in the spanwise direction was weaker than expected on the basis of linear theory. For both runs, the shear timescale h0u was set to 28.28 s to give dimensional values relevant to the ocean and to facilitate comparison with previous studies. The model timestep was set to Δt = 0.01 h0u.

3. Results

We begin this discussion by describing the basic patterns of temporal and spatial variability occurring in our two simulations. This description will be set mainly in the context of scales predicted previously via linear stability analysis. We will then quantify mixing in the two cases in terms of irreversible potential energy gain.

a. Temporal evolution

In each simulation, the kinetic energy associated with the random component of the initial flow decayed rapidly during the first few time steps (Fig. 4). This was followed by a longer period of two-dimensional growth as predicted by linear stability theory. Streamwise and vertical kinetic energies grew during this phase, while spanwise motions continued to decay. As expected, the Holmboe mode grew more slowly than the KH mode, and its growth was modulated by a pronounced oscillation. The short lines near the streamwise kinetic energy curves on Fig. 4 indicate exponential growth rates predicted using the linear stability analyses of the previous section. Both modes grew at rates close to those predicted by linear theory. (We do not expect exact correspondence in this case because the initial perturbations are not eigenfunctions and because the mean flow spreads somewhat in time due to viscosity and mass diffusion.)

Fig. 4.

Kinetic energy evolution for (a) case 1 and (b) case 2. Thick solid lines: kinetic energy of streamwise velocity perturbations. Dashed (dash–dotted) lines represent kinetic energy of spanwise (vertical) velocity components. Short, thin lines located near t = 0 indicate the exponential growth rate of the normal mode instability from linear theory. Steeper lines located near the spanwise kinetic energy curve indicate representative values of the growth rate of secondary instability from (a) KP (their Fig. 6; Ri = 0.16) and (b) SP (their Fig. 11a). Slopes correspond to 2σ, the growth rate for kinetic energy

Fig. 4.

Kinetic energy evolution for (a) case 1 and (b) case 2. Thick solid lines: kinetic energy of streamwise velocity perturbations. Dashed (dash–dotted) lines represent kinetic energy of spanwise (vertical) velocity components. Short, thin lines located near t = 0 indicate the exponential growth rate of the normal mode instability from linear theory. Steeper lines located near the spanwise kinetic energy curve indicate representative values of the growth rate of secondary instability from (a) KP (their Fig. 6; Ri = 0.16) and (b) SP (their Fig. 11a). Slopes correspond to 2σ, the growth rate for kinetic energy

Once the spanwise kinetic energy began to grow, its growth was considerably more rapid than that of the streamwise and vertical components. This is a manifestation of three-dimensional, secondary instabilities of the finite amplitude, two-dimensional waves. Linear, secondary instabilities of KH billows have been computed for the case Re0 = 1200, Pr = 1, Ri0 = 0.16 by Klaassen and Peltier (1991, hereinafter KP). [A higher-resolution study has been conducted recently by Caulfield and Peltier (2000), but the results are limited to Ri0 ≤ 0.1 and are thus of less relevance here.] Although the parameter values of KP are slightly different from ours, the growth rate predicted from their analyses compares well with the growth of the spanwise kinetic energy seen in Fig. 4a.

Secondary stability analysis for the Holmboe case simulated here has not yet been carried out. The closest existing results are those of Smyth and Peltier (1991, hereafter SP), which focused on much smaller values of the centerline Richardson number (Ri0 ≤ 0.28). That study predicted that growth rates of secondary instabilities would greatly exceed those of primary instabilities, and that result is also evident in Fig. 4b. In quantitative terms, though, the growth rate observed here (about 0.045 in nondimensional units) is far less than that predicted by SP (∼0.12). This quantitative discrepancy is perhaps not surprising given the significant difference in parameter values.

In each case, the spanwise kinetic energy saturates at a value nearly equal to the vertical kinetic energy. There follows an extended period in which the three components decay in parallel, with the spanwise and vertical components nearly equal and the streamwise component almost an order of magnitude larger.

b. Spatial structure

We now describe the evolving spatial structure of our simulated flows, beginning with the KH case. During the two-dimensional growth phase, the density field (Fig. 5a) reveals the standard “cat's eye” structure of the KH billow. The centerline Richardson number (Ri0 = 0.15) is high relative to previous KH simulations (e.g., Caulfield and Peltier 2000; Smyth and Moum 2000b), and as a result the billow is strongly elliptical. At the state of maximum kinetic energy (Fig. 5b), the density field is dominated at large scales by the KH billow; at small scales by disordered, three-dimensional motions concentrated in the vortex core. This chaotic, strongly dissipative flow state is referred to loosely here as “turbulence,” although the Reynolds number is barely large enough to allow such a state to exist. Although the low Reynolds numbers used here are necessary in order to maintain the required spatial resolution in the density field for these high Prandtl number flows, low Reynolds number events such as this may actually account for much of the mixing occurring in the ocean thermocline (Smyth et al. 2001).

Fig. 5.

(a), (b) Density structure of the KH billow as it approaches maximum amplitude, (c) as it evolves into a Holmboe wave, and (d) after one oscillation period of the Holmboe wave. The range of shading from dark to light corresponds to the range of density from high to low, but the correspondence is not exact because additional lighting effects are used to highlight three-dimensionality. The highest and lowest densities are rendered transparent, as is all fluid in the upper right-hand quarter of the billow in (a), (b), and (c)

Fig. 5.

(a), (b) Density structure of the KH billow as it approaches maximum amplitude, (c) as it evolves into a Holmboe wave, and (d) after one oscillation period of the Holmboe wave. The range of shading from dark to light corresponds to the range of density from high to low, but the correspondence is not exact because additional lighting effects are used to highlight three-dimensionality. The highest and lowest densities are rendered transparent, as is all fluid in the upper right-hand quarter of the billow in (a), (b), and (c)

Shortly after the state of maximum kinetic energy is reached, the KH wave separates into a pair of oppositely propagating vortices resembling the Holmboe wave (Figs. 5c,d; cf. Fig. 7). This wave persisted throughout the remainder of the simulation, although it never became as well-defined as the primary Holmboe wave that appeared in the explicit “Holmboe” simulation. This transition has been observed previously in laboratory experiments (Browand and Winant 1973; HI). Because the Prandtl number exceeds unity, the shear layer tends to spread more rapidly than does the stratified layer, with the result that the mean flow tends to evolve to a state susceptible to a Holmboe-like instability. [In the absence of instability, the scale ratio would approach the limiting value R = Pr (Smyth et al. 1988), which is one reason we chose {R = 3, Pr = 9} to characterize the Holmboe regime.] Close examination of the mean profiles (not shown) reveals that the thin stratified layer is actually divided into two sublayers separated by a central mixed layer. Caulfield (1994) has shown theoretically that this more complex flow geometry is, in fact, more unstable to Holmboe-like modes than is Holmboe's original, two-layer density distribution. A related phenomenon was noticed in simulations of more weakly stratified flow by Smyth and Moum (2000b): a KH billow in its final stages of decay spontaneously produced a secondary KH billow as the result of nonlinear modification of the gradient Richardson number profile.

Fig. 7.

Isopycnal representation of two-dimensional cross sections through the Holmboe wave at selected times. Isopycnal levels are {−0.95, −0.8, −0.5, −0.2, 0.2, 0.5, 0.8, 0.95} × Δρ/2

Fig. 7.

Isopycnal representation of two-dimensional cross sections through the Holmboe wave at selected times. Isopycnal levels are {−0.95, −0.8, −0.5, −0.2, 0.2, 0.5, 0.8, 0.95} × Δρ/2

The complex structure of the mature KH billow (Fig. 5b) is the finite amplitude expression of three-dimensional secondary instability of the two-dimensional primary billow as computed by KP. We now compare that structure with the predictions of KP, focusing first on variability in the y direction and then on structure in the xz plane.

The spanwise density gradient spectrum (Fig. 6a) shows that the dominant nondimensional wavenumber at t = 1979 s is consistent with the growth rate maximum predicted by KP for the fastest-growing mode. However, the degree of scale selectivity predicted by linear theory is not evident here; density fluctuations cover a wide band of spanwise wavenumbers. This suggests both upscale and downscale cascades due to the finite amplitude of this fully nonlinear secondary flow.

Fig. 6.

(a) Power spectrum of the nondimensional spanwise density gradient, ρy/ρ0, vs nondimensional spanwise wavenumber β, computed at t = 1979 s. The spectrum is averaged over {Lx/4 < x < 3Lx/4; − Lz/4 < z < Lz/4}; the spanwise wavenumber is nondimensionalized by h0/2 to give β. Inset in (a) shows the exponential growth rate spectrum of the fastest-growing secondary instability, replotted from KP (their Fig. 20). (b) Contours of K3D, the kinetic energy associated with the secondary instability. The kinetic energy corresponds to velocity fluctuations about the spanwise-averaged velocity, and is itself spanwise averaged. Contour values are {0.25, 0.50, 0.75, 1.00} × 3.2 × 10−9 J kg−1. The background is an image plot of the spanwise-averaged density field

Fig. 6.

(a) Power spectrum of the nondimensional spanwise density gradient, ρy/ρ0, vs nondimensional spanwise wavenumber β, computed at t = 1979 s. The spectrum is averaged over {Lx/4 < x < 3Lx/4; − Lz/4 < z < Lz/4}; the spanwise wavenumber is nondimensionalized by h0/2 to give β. Inset in (a) shows the exponential growth rate spectrum of the fastest-growing secondary instability, replotted from KP (their Fig. 20). (b) Contours of K3D, the kinetic energy associated with the secondary instability. The kinetic energy corresponds to velocity fluctuations about the spanwise-averaged velocity, and is itself spanwise averaged. Contour values are {0.25, 0.50, 0.75, 1.00} × 3.2 × 10−9 J kg−1. The background is an image plot of the spanwise-averaged density field

Next, we compare the xz structure of the simulated wave at t = 1979 s with that of the fastest-growing eigenmode of KP. We isolate motions associated with secondary instability by computing the kinetic energy associated with velocity fluctuations about the spanwise-averaged state:

 
formula

K3D is concentrated in the convectively unstable regions on the upper-left and lower-right edges of the billow core (Fig. 6b). This structure corresponds well with the linear eigenfunctions of KP (their Fig. 21).

We now investigate the spatial structure of the Holmboe wave. Figure 7 shows two-dimensional cross-sections through the wave at selected times covering approximately one oscillation period. At this point in its evolution, the wave is nearly two-dimensional and is similar to those simulated previously by Smyth et al. (1988), SP, and Lawrence et al. (1998). The standing-wave-like structure of the Holmboe wave may be understood intuitively if one thinks of the thin layer of intense stratification around z = 0 as a flexible barrier, at which vertical motions are strongly inhibited. Above and below this barrier are sheared regions that are relatively free of stratification, so we might expect KH-like instability in these regions. If the barrier were rigid (i.e., if Ri0 were infinite), there would be no such instability because neither half of the velocity profile contains an inflection point. Because the stratification is finite, however, coupling between the upper and lower halves of the shear layer is allowed, and this coupling induces a positive feedback that causes KH-like billows to grow in the upper and lower layers [see Baines and Mitsudera (1994) for a more complete description of this mechanism].

Figure 8 illustrates the three-dimensional geometry of the mature Holmboe wave at various points in its evolution. At t = 3817 s (Fig. 8a) the spanwise kinetic energy is still more than an order of magnitude smaller than the streamwise and vertical components (Fig. 4), but spanwise structure is already evident. Inspection suggests that the dominant spanwise length scale is Ly/4. At t = 4383 s (Fig. 8b), spanwise kinetic energy is not much smaller than the other two components, and the three-dimensionality of the density field is correspondingly more fully developed. By t = 4949 s (Fig. 8c), all three components of the kinetic energy have reached their maxima. The billows now exhibit a highly complex structure with variations over a wide range of length scales. As the large billows pass each other, strong overturning is generated in the intervening fluid such that the central density interface itself is overturned. This generates overturns that are highly localized in x and z (arrow on Fig. 8c).

Fig. 8.

Density structure of the Holmboe wave through approximately two oscillation cycles. The range of shading from dark to light corresponds to the range of density from high to low, but the correspondence is not exact because additional lighting effects are used to highlight three-dimensionality. The highest and lowest densities are rendered transparent, as are sections of the upper billow in (b) and (c)

Fig. 8.

Density structure of the Holmboe wave through approximately two oscillation cycles. The range of shading from dark to light corresponds to the range of density from high to low, but the correspondence is not exact because additional lighting effects are used to highlight three-dimensionality. The highest and lowest densities are rendered transparent, as are sections of the upper billow in (b) and (c)

Figure 8c represents the last oscillation cycle visible in the kinetic energy evolution (Fig. 4) before the flow enters the turbulence decay phase. At later times (Fig. 8d), the oscillatory character of the large-scale wave is still evident in the density field, but the billows are severely distorted and there is very little displacement of the central density interface. Because of the latter, changes in available potential energy (and thus in kinetic energy) associated with the oscillation are small. There is, however, considerable activity in the upper and lower parts of the shear layer, and stratified fluid is ejected far above and below the original stratified layer. These ejections often take the form of loops, one of which is visible near the upper right-hand corner of Fig. 8d. Closer inspection (Fig. 9) shows that these structures are similar to the hairpin vortices that appear in boundary layers (e.g., Head and Bandyopadhyay 1981). Similar structures have been noticed in laboratory experiments on Holmboe waves (G. Lawrence 2001, personal communication). They form when localized vertical motions (e.g., the overturn shown by the arrow on Fig. 8c) displace isopycnals from their equilibrium depths, and the displaced fluid is differentially advected by the vertically sheared streamwise flow. This differential advection of the upper and lower parts of the loop tends to amplify vorticity in the loop via vortex stretching. [See Lin and Corcos (1984) for a more detailed discussion.]

Fig. 9.

Vortex lines associated with the loop structure seen in Fig. 8c (also Fig. 10b). Domain boundaries are shaded according to the vorticity component normal to each boundary, with dark (light) indicating clockwise (counterclockwise) rotation

Fig. 9.

Vortex lines associated with the loop structure seen in Fig. 8c (also Fig. 10b). Domain boundaries are shaded according to the vorticity component normal to each boundary, with dark (light) indicating clockwise (counterclockwise) rotation

The process of loop formation is illustrated in Fig. 10. Figure 10b represents the same stage of flow evolution as Fig. 8c (also Fig. 9), but more layers of low density are rendered transparent in order to expose the localized overturn. Figure 10a shows a similar view taken at a slightly earlier time. Between Figs. 10a and 10b the upper part of the overturn has not moved significantly in the x direction (because its vertical location is near z = 0), but the lower billow has been advected to the left so that it straddles x = 0, and part of it has reentered near the lower right-hand corner of the domain (cf. Fig. 9). The straining motion of the background flow has extruded the overturn into a loop whose apex is shown by the thick arrows. The thinner arrow on Fig. 10b indicates a region near the base of the loop where overturns in the yz plane reveal the counterrotating structures visible in the vorticity field (Fig. 9).

Fig. 10.

Density structure of the Holmboe wave at two times, illustrating loop formation. Shading is as in Fig. 5. Arrows indicate different regions of the loop structure discussed in the text

Fig. 10.

Density structure of the Holmboe wave at two times, illustrating loop formation. Shading is as in Fig. 5. Arrows indicate different regions of the loop structure discussed in the text

A simple model for a one-dimensional vortex maintained against viscosity by stretching is the Burgers vortex (Batchelor 1967), whose cross section is a Gaussian function and whose diameter is 4ν/S, where S is an imposed strain field that stretches the vortex along its axis. The Burgers vortex has been used previously as a model for one-dimensional vortices that appear in homogeneous turbulence (Jimenez 1992; Andreotti 1997; Jimenez and Wray 1998). In the present context, S may be approximated by the strain due to the horizontally averaged streamwise flow u(z), that is, S = 0.5 du/dz. We approximate du/dz as the velocity difference 0.5Δu ≈ 3 × 10−3 m s−1 divided by a typical height 0.3 m, so that S ≈ 0.01 s−1. With ν ≈ 10−6 m2 s−1, the diameter of a Burgers vortex becomes 0.06 m. This is a reasonable approximation to the diameters of the loops seen in Figs. 8d and 10. The spatial structure of the loops seen here is therefore consistent with the notion that they are shaped by a competition between straining by the mean shear and viscosity. Mass diffusion and baroclinic torques due to convective overturning are also likely to be important factors.

The loop shown in Fig. 10 is distinct from that seen previously in Fig. 8d. The former dissipated shortly after t = 5000 s without extending very far in the vertical; the latter formed later in a higher region of the flow and transported stratified fluid far above its original depth. In even more extreme cases, these ejections reached the upper and lower boundaries. A loop structure also appeared in the later stages of the KH evolution (Fig. 5d), but did not eject fluid nearly as far from the shear layer as did the loops arising from the Holmboe wave.

We now describe the spatial structure of the mature Holmboe wave with reference to the secondary stability analyses of SP. Once again, those analyses were performed at a significantly smaller level of stratification than the present simulations, so we do not expect quantitative correspondence.

Smyth and Peltier (1991) found that the structure of the dominant secondary instability varied significantly with the phase of the Holmboe wave (e.g., their Fig. 10). At some points, the dominant mode was concentrated around highly localized overturns on the central density interface; at other points, instability was stronger in the upper and lower billows. At any given instant, the growth rate of secondary instability was large enough compared with the oscillation frequency of the wave that the characterization of the instability as a spatially self-similar, exponentially growing normal mode was self-consistent. Nevertheless, the variability of the mode shape indicates that the secondary disturbance existing at any given time contains contributions from modes with different shapes that emerged at previous times. Smyth and Peltier (1991) documented the dependence of growth rate on spanwise wavenumber for only a single mode, namely, the highly localized mode that appeared when the central density interface was overturned. (This mode exhibited the largest exponential growth rate.) This dependence gave the impression of an “ultraviolet catastrophe”; that is, the growth rate increased monotonically with increasing wavenumber up to the highest wavenumber resolved in the numerical calculation. Such a mode is thought to generate a direct transfer of energy into the dissipation range of the wavenumber spectrum. It was later suggested by Smyth and Peltier (1994) that this behavior is extremely sensitive to the time dependence of the background state: modes with small length scales could be effectively decorrelated by even slow changes in the background flow. If this suggestion is correct, the ultraviolet catastrophe will appear only on background flows that, unlike the Holmboe wave, are very nearly stationary.

Figure 11 shows the structure of the Holmboe wave at approximately the same phase in three successive oscillation cycles approaching the maximum amplitude state. The left-hand column shows the spanwise density gradient spectrum; the inset in Fig. 11c shows SP's growth rate curve for the fastest-growing secondary instability. As the wave grows, the spectrum broadens, indicating both upscale and downscale cascades due to nonlinear mode–mode interactions. On the basis of linear theory, we would expect the shape of the growth rate curve to be reflected in the spectrum of the simulated flow, at least in the early stages of secondary instability growth. This expectation is not fulfilled. The growth rate curve suggests that spanwise variability will be found mainly around β = 4–5 and possibly higher, whereas most variability in the simulated fields occurs at lower β. At the earliest times (when linear theory should be most relevant) the strongest fluctuations occur at β < 3. It could be that this is partly a function of the initialization; although the initial noise field is white, viscous dissipation of that noise in the early part of the simulations occurs preferentially at the smallest scales so that those scales begin the three-dimensional growth phase at lower amplitude. Nevertheless, we conclude that the growth rate spectrum of SP gives an unsatisfactory description of the present results, even in qualitative terms. In particular, there is no evidence here for ultraviolet catastrophe in Holmboe waves.

Fig. 11.

Three-dimensional aspects of the Holmboe wave. (a), (c), (e) Power spectrum of the spanwise gradient of the nondimensional density, ρy/ρ0, vs nondimensional spanwise wavenumber β, computed at t = 4666 s. The spectrum is averaged over all x and z; the spanwise wavenumber is nondimensionalized by h0/2 to give β. Inset in (c) shows the exponential growth rate spectrum of secondary instability for the case Ri0 = 0.28; replotted from SP, Fig. 9a, t = 120. (b), (d), (f) Contours of K3D, the kinetic energy associated with the secondary instability (19). Contour values are {0.25, 0.50, 0.75, 1.00} × Ksc, where Ksc = 3.7 × 10−9 for (b), 4.3 × 10−8 for (d), and 3.4 × 10−7 J kg−1 for (f). The background is an image plot of the spanwise-averaged density field

Fig. 11.

Three-dimensional aspects of the Holmboe wave. (a), (c), (e) Power spectrum of the spanwise gradient of the nondimensional density, ρy/ρ0, vs nondimensional spanwise wavenumber β, computed at t = 4666 s. The spectrum is averaged over all x and z; the spanwise wavenumber is nondimensionalized by h0/2 to give β. Inset in (c) shows the exponential growth rate spectrum of secondary instability for the case Ri0 = 0.28; replotted from SP, Fig. 9a, t = 120. (b), (d), (f) Contours of K3D, the kinetic energy associated with the secondary instability (19). Contour values are {0.25, 0.50, 0.75, 1.00} × Ksc, where Ksc = 3.7 × 10−9 for (b), 4.3 × 10−8 for (d), and 3.4 × 10−7 J kg−1 for (f). The background is an image plot of the spanwise-averaged density field

The xz structure of the three-dimensional motions is shown in Figs. 11b,d,f via contour diagrams of the kinetic energy K3D superimposed on the spanwise-averaged density. The xz dependence seen here is qualititively consistent with the predictions of SP. Early in the three-dimensional growth phase (Fig. 11b), secondary circulations are concentrated in the lower billow of the Holmboe wave. This asymmetry reflects some slight, random asymmetry in the particular noise field introduced into this flow realization at t = 0. By t = 4807 s (Fig. 11f), K3D has grown by two orders of magnitude and has become much more symmetric. For the most part K3D is concentrated in the billows of the Holmboe wave. Also evident, however, is a highly localized disturbance (slightly left of center in Fig. 11f) associated with the small overturn in Fig. 8c. This structure evidently corresponds to the intensely localized convective instability found by SP (their Fig. 10, e.g., t = 120, 160).

The density gradient spectrum shown in Fig. 11e provides a particularly stringent test of the spatial resolution of the density field (cf. discussion in sec. 2b), as it corresponds to an advanced state of turbulence development. The fact that the spectrum drops smoothly over several orders of magnitude at the highest wavenumbers shows that the density field is well resolved.

c. Mixing

Mixing is usefully quantified in terms of an irreversible increase of potential energy (Winters et al. 1995; Winters and D'Asaro 1996). The computation of the irreversible potential energy increase begins with the identification of the minimum potential energy state of a given density distribution. In this state, density varies only in the vertical, and its vertical gradient is negative semidefinite. The associated (volume averaged) potential energy is defined, up to an additive constant, by

 
formula

in which ρ∗(z) is the density profile defining the equilibrium state. The total (volume averaged) potential energy is

 
formula

where the integral is over the computational volume. The available potential energy is just the difference: Pa = PtPb.

A fluid volume obeying the boundary conditions (5)–(8), at rest in the equilibrium state, gains potential energy by molecular diffusion at the (volume averaged) rate

 
formula

in which ρ|topbottom denotes the net difference in horizontally averaged density between the upper and lower boundaries (Winters et al. 1995). The energy source for this increase in background potential energy is the internal energy of the rest state.

Total and background potential energies for the KH and Holmboe cases are shown in Figs. 12a and 12b. Note that the growth of Pb is monotonic. In contrast, Pt fluctuates in response to large-scale wavelike motions, especially in the Holmboe case (Fig. 12b). In each case, the total potential energy growth is more rapid than the irreversible component in the early part of the simulation. However, Pb grows more rapidly once the flow has become fully turbulent, and Pt and Pb become equal once turbulent motions have decayed. Thin, solid curves on Figs. 12a and 12b show the time integral of Φ. This represents a substantial contribution to the potential energy increase in these low Reynolds number flows. By the time Pt and Pb become nearly equal, the rate of potential energy increase is only negligibly larger than Φ.

Fig. 12.

Mixing in KH and Holmboe waves. (a), (b) Total potential energy Pt (solid curves); background potential energy Pb (dashed curves), potential energy increase in the absence of fluid motions t0 Φ(t′) dt′ (thin, solid curves). (c), (d) Viscous dissipation rate ε (solid); rate of irreversible potential energy creation by fluid motions (dashed). (e), (f) Instantaneous flux coefficient. Horizontal lines indicate Γi = 0.2. Some curves are jagged because data were not saved frequently enough to thoroughly resolve the fast oscillations

Fig. 12.

Mixing in KH and Holmboe waves. (a), (b) Total potential energy Pt (solid curves); background potential energy Pb (dashed curves), potential energy increase in the absence of fluid motions t0 Φ(t′) dt′ (thin, solid curves). (c), (d) Viscous dissipation rate ε (solid); rate of irreversible potential energy creation by fluid motions (dashed). (e), (f) Instantaneous flux coefficient. Horizontal lines indicate Γi = 0.2. Some curves are jagged because data were not saved frequently enough to thoroughly resolve the fast oscillations

To calculate the mechanical efficiency of mixing, we would relate the rate of irreversible potential energy increase to the rate at which kinetic energy is supplied to the turbulence by the mean flow. Unfortunately, there is no way to distinguish latter from the reversible transfer of kinetic energy between the mean flow and wavelike, nonturbulent motions. Instead, we represent mixing “efficiency” by the flux coefficient Γi, which relates the rate of irreversible potential energy increase to ε, the rate at which kinetic energy is lost to viscous dissipation (Moum 1996), namely,

 
formula

The dissipation rate is given explicitly by ε = ν Σ3i,j=1 (ui,j + uj,i)2/2 and angle brackets denote a volume average over the computational domain. In this computation, Φ is subtracted from dPb/dt in order to isolate the rate of irreversible potential energy increase associated with fluid motions. The subscript i on Γi indicates the instantaneous value, as opposed to the ratio of net potential energy gain to net dissipation (Smyth et al. 2001).

In both simulations, we observe an initial phase in which dPb/dt − Φ and 〈ε〉 are relatively small and nearly equal (Figs. 12c,d). This is the period of laminar, two-dimensional growth. The transition to turbulence is marked by a dramatic increase in 〈ε〉, accompanied by a much smaller increase in dPb/dt − Φ. As turbulence decays, the two quantities decrease approximately in proportion to one another.

The approximate equality of dPb/dt − Φ and 〈ε〉 during the two-dimensional growth phase is expressed in values of the flux coefficient near unity. A central result of this work is that Holmboe waves, despite their dramatically different physical character, exhibit values of Γi near unity in the preturbulent phase just as do KH waves (Figs. 12e,f). In fact, values of Γi are slightly larger in the Holmboe case. Because of its low growth rate, the phase of highly efficient preturbulent mixing lasts much longer in the Holmboe case than in the KH case.

Close inspection shows that peaks in Γi occur just as the upper and lower vortices pass. We suspect that this is because the close interaction of the vortices creates exceptionally sharp density gradients in the intervening fluid (e.g., Figs. 7a,c). The localized instability seen in Figs. 8b and 8c may represent a significant source of mixing at this stage in the oscillation cycle.

With the transition to turbulence, Γi decreases to near 0.2 in both flows. In the KH case, this decrease is only temporary; Γi begins to increase again at t = 2800 s and reaches 0.5 by the end of the simulation. By this time, Γi has become a ratio of small numbers, and the large values may be due in part to roundoff error. Another reason for the increase, however, is the emergence of the secondary Holmboe wave (Figs. 8c,d).

Comparison of Figs. 12a and 12b reveals a second significant result of this study. Linear theory suggests that Holmboe waves will be less effective at mixing the fluid than KH waves because the instability they grow from is relatively weak. Here, we find that the Holmboe wave does grow more slowly than the KH wave, but that the net amount of mixing (i.e., irreversible potential energy gain) is quite comparable. In fact, by the time turbulence has decayed (i.e., when Pt and Pb become equal and dPb/dt becomes indistinguishable from Φ), the Holmboe wave has accomplished somewhat more potential energy increase than has the KH wave. By virtue of its rapid growth, the KH instability quickly consumes the kinetic energy stored in the mean shear. The Holmboe wave, in contrast, grows more slowly but as a result lasts longer, and ultimately accomplishes a comparable amount of mixing.

This comparison becomes more vivid when potential energy gain due to fluid motion is isolated by subtracting out the time integral of Φ (Fig. 13). Although the KH wave mixes the fluid more rapidly, the net potential energy gain attributable to the Holmboe wave is about twice as large as that due to the KH wave. In each case, slightly less than one half of the potential energy gain due to fluid motions occurred before the transition to turbulence.

Fig. 13.

Potential energy gain due to fluid motion in KH and Holmboe waves. Thick sections of curves represent the preturbulent phase. This phase is taken to end the first time that Γi drops below 0.3 (cf. Fig. 12c)

Fig. 13.

Potential energy gain due to fluid motion in KH and Holmboe waves. Thick sections of curves represent the preturbulent phase. This phase is taken to end the first time that Γi drops below 0.3 (cf. Fig. 12c)

Net changes in mean profiles for the KH and Holmboe cases are shown in Fig. 14. Initial velocity profiles (thin curves) are the same in both cases, and in both cases the shear layer spreads considerably during the simulation. Spreading is greater in the Holmboe case, due in part to the longer duration of the run. Initial density stratification near z = 0 is greater in the Holmboe case than in the KH case and remains greater as the flows evolve. However, density mixing away from z = 0 is more similar in the two cases; both Holmboe and KH instabilities alter the density profile significantly out to z ∼ ±0.2 m. This is a result of strong plume ejection events in the Holmboe case, which are in turn associated with the tendency toward reduced Richardson numbers on the edges of the shear layer as illustrated in Fig. 2. Even in the KH case, the shear layer spreads more rapidly than does the density interface, accounting for the Holmboe-like behavior seen in the late stages of that simulation (Figs. 5c,d). However, this flow retains a significant feature of turbulent KH waves that appears not to develop the Holmboe case: a thin region of well-mixed fluid near z = 0. This tendency for enhanced mixing of mass near the centerline has been noted previously by (e.g.) Scinocca (1995) and Caulfield and Peltier (2000), and was found by Smyth and Moum (2000b) to lead to the generation of new KH billows long after turbulence from the original event had decayed.

Fig. 14.

Changes in mean velocity and density profiles due to KH and Holmboe instabilities. Thick (thin) curves indicate final (original) profiles

Fig. 14.

Changes in mean velocity and density profiles due to KH and Holmboe instabilities. Thick (thin) curves indicate final (original) profiles

4. Conclusions

We have compared the evolution of turbulence in two stratified shear layers that differ only in the thickness of the density interface. (Note that this differs from the approach usually taken in laboratory experiments, in which the bulk Richardson number is varied by varying Δρ; e.g., HI.) One case yields KH instability, whose mixing characteristics are relatively well understood. The second case yields Holmboe waves, which have not been simulated in three dimensions before now. The KH billow simulated here became three-dimensional via a mechanism very similar to the secondary instability predicted by KP. After the KH billow had reached large amplitude, it evolved into a Holmboe-like standing wave. This transition is expected to be common in high Prandtl number fluids, a result that emphasizes the potential importance of Holmboe waves in the ocean.

The Holmboe wave was found to develop turbulence despite the low growth rate of the primary instability. Secondary stability analysis of Holmboe waves at lower Ri0 (SP) predicted some features of the transition to turbulence in the Holmboe wave simulated here, but significant differences are also apparent. Loop structures in the density field associated with hairpin-like vortices are a conspicuous feature of turbulent Holmboe waves. These structures are initiated by secondary instabilities (in one case this resembled the localized convective instability described by SP) and grow to large amplitude via vortex stretching.

Our main objective has been to compare the Holmboe and KH waves in terms of irreversible potential energy gain. This comparison has led us to two significant conclusions.

  1. Despite their very different spatial structure, Holmboe waves exhibit highly efficient mixing in the laminar growth phase that precedes the transition to turbulence just as do KH waves. This tends to confirm our expectation that high mixing efficiency in the preturbulent regime is characteristic of a wide range of unstable shear flows. High efficiency appears to be a result of the coherence of the strain field in laminar flow as compared with turbulent flow. Turbulent strain fields tend to be stronger, but they are also highly chaotic and therefore unable to advect the density gradient into the optimal orientation for mixing (Smyth 1999). At the same time, the high strain rate means that kinetic energy is dissipated rapidly by viscosity. Laminar strain fields, though weaker, are relatively steady and are therefore able to compress density gradients effectively while minimizing viscous dissipation.

  2. Although they result from a relatively “weak” instability (in terms of linear growth rate), Holmboe waves mix the fluid at least as effectively as do KH waves. This suggests that the traditional focus on the fastest-growing modes of linear theory is not always appropriate. Instability in stratified shear flows is the mechanism by which an unstable mean flow adjusts to a stable state. This adjustment process invariably involves some degree of mixing, as a part of the kinetic energy shed by the mean flow is converted irreversibly to potential energy. The remaining kinetic energy is lost to viscosity (or, in other situations, radiated as waves). The net amount of mixing depends not on the vigor of the linear instability, but rather on the fraction of mean flow kinetic energy that winds up as potential energy. A slowly growing instability may spend more time in the highly efficient preturbulent state and thus, given sufficient time, mix the fluid more completely than does an instability that extracts mean kinetic energy rapidly but wastes most of it. In other words, Holmboe instability may mix effectively not despite its low growth rate but because of it.

5. Future work

The work described here must be extended in several ways to attain a comprehensive understanding of Holmboe instability and its role in ocean mixing. The secondary stability analyses of SP must be extended considerably if we are to understand the transition process in Holmboe waves in terms of linear modes. In addition, nonlinear studies are limited in several respects due to finite computational resources, as we now discuss.

Our simulated flows have developed significant structure at the largest scales permitted by the boundary conditions, suggesting the possibility of undesirable boundary effects. In the streamwise direction, this is unlikely to be a problem: upscale energy cascade via KH vortex pairing happens only when stratification is much weaker than that used here (e.g., Klaassen and Peltier 1989; Smyth and Moum 2000b), and two-dimensional tests indicate that Holmboe waves do not pair (at least not in this region of parameter space). In the spanwise direction, an extended domain would be preferable. This has not been possible with the resources presently available. However, auxiliary simulations run with Ly halved yielded essentially the same results seen here: transition to turbulence in both KH and Holmboe cases, hairpin vortices, efficient mixing in the preturbulent phase and net potential energy gain greater in the Holmboe case. In the vertical, some ejection events occurring late in the Holmboe simulation encountered the boundaries at z = ±Lz/2, but it seems unlikely that these interactions influenced the rest of the flow.

The Reynolds numbers achieved in the present work lie at the low end of the range of values for turbulent events in the ocean interior (Smyth et al. 2001). They are considerably smaller than those occurring in exchange flows and river outflows and are orders of magnitude below the values at which standard theories of turbulence apply. Future simulations will use a newly developed numerical model in which the scalar field is discretized on a finer grid than the velocity fields, optimizing the use of computer memory for high-Prandtl-number flows. While this gives access to higher Reynolds numbers, values will remain at the low end of the geophysical range for the forseeable future. We do not believe, however, that the main results of the present study are likely to be sensitive to Reynolds number. The evolution of Γi has been checked explicitly for Reynolds number dependence by Smyth et al. (2001) for the case of KH waves, and no such dependence was found. There is no reason to expect that the Holmboe wave will differ in this regard. Similarly, there is no reason to expect that the relative effectiveness of Holmboe and KH waves at mixing the fluid should depend upon Reynolds number. To summarize, future studies must employ larger domains and higher Reynolds numbers as resources permit, but it is unlikely that such studies will alter the conclusions reached here.

The present comparison has covered only a single point in the Holmboe wave parameter space; results may differ at different values of Ri0 and R. Also, we have assumed that the initial flow is symmetric, in the sense that the center of the stratified layer coincides with the center of the shear layer. This is a special case of a more general situation in which the centers of the sheared and stratified layers may be offset from one another (Haigh and Lawrence 1999). As the offset increases from zero, the gradient Richardson number decreases in the half of the shear layer away from the stratified layer, and the component of the disturbance that is centered in that region thus becomes stronger. When the offset becomes large, the stronger disturbance takes the form of a KH instability with small Ri0. The results of Staquet (2000) suggest that mixing efficiency in KH waves increases as Ri0 decreases, and therefore that the limit Ri0 → 0 is singular (because mixing efficiency must vanish in unstratified flow). Dependence of mixing physics on Ri0, R and asymmetry must be explored before the present results can be applied with confidence to more general classes of oceanic flows.

Acknowledgments

The research was funded by the National Science Foundation under Grants OCE0097124, OCE0221057 (WDS) and OCE9617761 (KBW). Time on the Cray T3E at the University of Texas was provided through an allocation from the National Partnership for Advanced Computing Infrastructure.

REFERENCES

REFERENCES
Andreotti
,
B.
,
1997
:
Studying Burgers models to investigate the physical meaning of the alignments statistically observed in turbulence.
Phys. Fluids
,
9A
,
735
742
.
Baines
,
P.
, and
H.
Mitsudera
,
1994
:
On the mechanism of shear flow instabilities.
J. Fluid Mech.
,
276
,
327
342
.
Batchelor
,
G. K.
,
1959
:
Small-scale variation of convected quantities like temperature in turbulent fluid.
J. Fluid Mech.
,
5
,
113
133
.
Batchelor
,
G. K.
,
1967
:
An Introduction to Fluid Dynamics.
Cambridge University Press, 615 pp
.
Browand
,
F.
, and
C.
Winant
,
1973
:
Laboratory observations of shear instability in a stratified fluid.
Bound.-Layer Meteor.
,
5
,
67
77
.
Caulfield
,
C.
,
1994
:
Multiple linear instability of layered stratified shear flow.
J. Fluid Mech.
,
258
,
255
285
.
Caulfield
,
C.
, and
W.
Peltier
,
1994
:
Three-dimensionalization of the stratified mixing layer.
Phys. Fluids
,
6A
,
3803
3805
.
Caulfield
,
C.
, and
W.
Peltier
,
2000
:
Anatomy of the mixing transition in homogeneous and stratified free shear layers.
J. Fluid Mech.
,
413
,
1
47
.
Cortesi
,
A.
,
G.
Yadigaroglu
, and
S.
Banerjee
,
1998
:
Numerical investigation of the formation of three-dimensional structures in stably-stratified mixing layers.
Phys. Fluids
,
10A
,
1449
1473
.
Cortesi
,
A.
,
B.
Smith
,
G.
Yadigaroglu
, and
S.
Banerjee
,
1999
:
Numerical investigation of the entrainment and mixing processes in neutral and stably-stratified mixing layers.
Phys. Fluids
,
11A
,
162
185
.
DeSilva
,
I.
,
H.
Fernando
,
F.
Eaton
, and
D.
Hebert
,
1996
:
Evolution of Kelvin-Helmholtz billows in nature and laboratory.
Earth Planet. Sci. Lett.
,
143
,
217
231
.
Farmer
,
D.
, and
J.
Smith
,
1980
:
Tidal interaction of stratified flow with a sill in Knight Inlet.
Deep-Sea Res.
,
27
,
239
254
.
Farmer
,
D.
, and
L.
Armi
,
1999
:
Stratified flow over topography: The role of small-scale entrainment and mixing in flow establishment.
Proc. Roy. Soc. London
,
455A
,
3221
3258
.
Gregg
,
M. C.
,
1987
:
Diapycnal mixing in the thermocline: A review.
J. Geophys. Res.
,
92
(
C5
)
5249
5286
.
Haigh
,
S.
, and
G.
Lawrence
,
1999
:
Symmetric and nonsymmetric Holmboe instabilities in an inviscid flow.
Phys. Fluids
,
11
,
1459
1468
.
Hazel
,
P.
,
1972
:
Numerical studies of the stability of inviscid parallel shear flows.
J. Fluid Mech.
,
51
,
39
62
.
Head
,
M.
, and
P.
Bandyopadhyay
,
1981
:
New aspects of turbulent boundary layer structure.
J. Fluid Mech.
,
107
,
297
338
.
Hebert
,
D.
,
J.
Moum
,
C.
Paulson
, and
D.
Caldwell
,
1992
:
Turbulence and internal waves at the equator. Part II: Details of a single event.
J. Phys. Oceanogr.
,
22
,
1346
1356
.
Hogg
,
A. M.
, and
G.
Ivey
,
2002
:
The Kelvin–Helmholtz to Holmboe instability transition in stratified exchange flows.
J. Fluid Mech., in press
.
Holmboe
,
J.
,
1962
:
On the behaviour of symmetric waves in stratified shear layers.
Geophys. Publ.
,
24
,
67
113
.
Howard
,
L.
,
1961
:
Note on a paper of John W. Miles.
J. Fluid Mech.
,
10
,
509
512
.
Jimenez
,
J.
,
1992
:
Kinematic alignment effects in turbulent flows.
Phys. Fluids
,
4A
,
652
654
.
Jimenez
,
J.
, and
A. A.
Wray
,
1998
:
On the characteristics of vortex filaments in isotropic turbulence.
J. Fluid Mech.
,
373
,
255
285
.
Klaassen
,
G.
, and
W.
Peltier
,
1985a
:
The evolution of finite-amplitude Kelvin–Helmholtz billows in two spatial dimensions.
J. Atmos. Sci.
,
42
,
1321
1339
.
Klaassen
,
G.
, and
W.
Peltier
,
1985b
:
The onset of turbulence in finite-amplitude Kelvin–Helmholtz billows.
J. Fluid Mech.
,
155
,
1
35
.
Klaassen
,
G.
, and
W.
Peltier
,
1989
:
The role of transverse secondary instabilities in the evolution of free shear layers.
J. Fluid Mech.
,
202
,
367
402
.
Klaassen
,
G.
, and
W.
Peltier
,
1991
:
The influence of stratification on secondary instability in free shear layers.
J. Fluid Mech.
,
227
,
71
106
.
Koop
,
C.
, and
F.
Browand
,
1979
:
Instability and turbulence in a stratified layer with shear.
J. Fluid Mech.
,
93
,
135
159
.
Lawrence
,
G.
,
F.
Browand
, and
L.
Redekopp
,
1991
:
Stability of a sheared density interface.
Phys. Fluids
,
3A
,
2360
2370
.
Lawrence
,
G.
,
S.
Haigh
, and
Z.
Zhu
,
1998
:
In search of Holmboe's instability.
Physical Processes in Lakes and Oceans, J. Imberger, Ed., Amer. Geophys. Union, 295–304
.
Lin
,
S.
, and
G.
Corcos
,
1984
:
The mixing layer: Deterministic models of a turbulent flow. Part 3: The effect of plane strain on the dynamics of streamwise vortices.
J. Fluid Mech.
,
141
,
139
187
.
Miles
,
J.
,
1961
:
On the stability of heterogeneous shear flows.
J. Fluid Mech.
,
10
,
496
508
.
Moin
,
P.
, and
K.
Mahesh
,
1998
:
Direct numerical simulation: A tool in turbulence research.
Annu. Rev. Fluid Mech.
,
30
,
539
578
.
Moum
,
J.
,
1996
:
Energy-containing scales of turbulence in the ocean thermocline.
J. Geophys. Res.
,
101
(
C6
)
14095
14109
.
Patnaik
,
P.
,
F.
Sherman
, and
G.
Corcos
,
1976
:
A numerical simulation of Kelvin–Helmholtz waves of finite amplitude.
J. Fluid Mech.
,
73
,
215
240
.
Pouliquen
,
S.
,
J.
Chomaz
, and
P.
Huerre
,
1994
:
Propagating Holmboe waves at the interface between two immiscible fluids.
J. Fluid Mech.
,
266
,
277
302
.
Scinocca
,
J.
,
1995
:
The mixing of mass and momentum by Kelvin–Helmholtz billows.
J. Atmos. Sci.
,
52
,
2509
2530
.
Seim
,
H.
, and
M.
Gregg
,
1994
:
Detailed observations of a naturally-occurring shear instability.
J. Geophys. Res.
,
99
(
C5
)
10049
10073
.
Smyth
,
W.
,
1999
:
Dissipation range geometry and scalar spectra in sheared, stratified turbulence.
J. Fluid Mech.
,
401
,
209
242
.
Smyth
,
W.
, and
W.
Peltier
,
1989
:
The transition between Kelvin–Helmholtz and Holmboe instability: An investigation of the overreflection hypothesis.
J. Atmos. Sci.
,
46
,
3698
3720
.
Smyth
,
W.
, and
W.
Peltier
,
1990
:
Three-dimensional primary instabilities of a stratified, dissipative, parallel flow.
Geophys. Astrophys. Fluid Dyn.
,
52
,
249
261
.
Smyth
,
W.
, and
W.
Peltier
,
1991
:
Instability and transition in Kelvin–Helmholtz and Holmboe waves.
J. Fluid Mech.
,
228
,
387
415
.
Smyth
,
W.
, and
W.
Peltier
,
1994
:
Three-dimensionalization of barotropic vortices on the f-plane.
J. Fluid Mech.
,
265
,
25
64
.
Smyth
,
W.
, and
J.
Moum
,
2000a
:
Anisotropy of turbulence in stably stratified mixing layers.
Phys. Fluids
,
12
,
1343
1362
.
Smyth
,
W.
, and
J.
Moum
,
2000b
:
Length scales of turbulence in stably stratified mixing layers.
Phys. Fluids
,
12
,
1327
1342
.
Smyth
,
W.
,
G.
Klaassen
, and
W. R.
Peltier
,
1988
:
Finite amplitude Holmboe waves.
Geophys. Astrophys. Fluid Dyn.
,
43
,
181
222
.
Smyth
,
W.
,
J.
Moum
, and
D.
Caldwell
,
2001
:
The efficiency of mixing in turbulent patches: Inferences from direct simulations and microstructure observations.
J. Phys. Oceanogr.
,
31
,
1969
1992
.
Staquet
,
C.
,
2000
:
Mixing in a stably stratified shear layer: Two- and three-dimensional numerical experiments.
Fluid. Dyn. Res.
,
27
,
367
404
.
Strang
,
E.
, and
H.
Fernando
,
2001
:
Entrainment and mixing in stratified shear flows.
J. Fluid Mech.
,
428
,
349
386
.
Thorpe
,
S.
,
1973
:
Turbulence in stably stratified fluids: A review of laboratory experiments.
Bound-Layer Meteor.
,
5
,
95
.
Thorpe
,
S.
,
1985
:
Laboratory observations of secondary structures in Kelvin–Helmholtz billows and consequences for ocean mixing.
Geophys. Astrophys. Fluid Dyn.
,
34
,
175
199
.
Thorpe
,
S.
,
1987
:
Transition phenomena and the development of turbulence in stratified fluids.
J. Geophys. Res.
,
92
,
5231
5245
.
Winters
,
K.
, and
E. A.
D'Asaro
,
1996
:
Diascalar flux and the rate of fluid mixing.
J. Fluid Mech.
,
317
,
179
193
.
Winters
,
K.
,
P.
Lombard
,
J.
Riley
, and
E. A.
D'Asaro
,
1995
:
Available potential energy and mixing in density-stratified fluids.
J. Fluid Mech.
,
289
,
115
128
.
Yonemitsu
,
N.
,
G.
Swaters
,
N.
Rajaratnam
, and
G.
Lawrence
,
1996
:
Shear instabilities in arrested salt-wedge flows.
Dyn. Atmos. Oceans
,
24
,
173
182
.
Yoshida
,
S.
,
M.
Ohtani
,
S.
Nishida
, and
P.
Linden
,
1998
:
Mixing processes in a highly stratified river.
Physical Processes in Lakes and Oceans, J. Imberger, Ed., Amer. Geophys. Union, 389–400
.
Zhu
,
D.
, and
G.
Lawrence
,
2001
:
Holmboe's instability in exchange flows.
J. Fluid Mech.
,
429
,
391
409
.

Footnotes

Corresponding author address: Dr. W. D. Smyth, College of Oceanic and Atmospheric Sciences, Oregon State University, Corvallis, OR 97331. Email: smyth@coas.oregonstate.edu