## 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.

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*:

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} m^{2} s^{−1}.

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

and the density evolves in accordance with

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

We assume periodicity in the horizontal dimensions:

in which *f* is any solution field and the periodicity intervals *L*_{x} and *L*_{y} are constants. At the upper and lower boundaries (*z* = ±*L*_{z}/2), we impose an impermeability condition on the vertical velocity:

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

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

### 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 *L*_{K} = (*ν*/*ε*)^{1/4}, where *ε* is the kinetic energy dissipation rate. Grid spacing of 4*L*_{K} is generally sufficient: Moin and Mahesh (1998) list several examples of successful simulations with grid spacing ranging from 3.75*L*_{K} to 14.3*L*_{K}. (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, *L*_{B} = *L*_{K}/Pr^{1/2} (Batchelor 1959). Previous studies of stratified shear layers with high Prandtl number have used grid spacing ≤2.5*L*_{B} (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 *L*_{B}, 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 2*L*_{B}. (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

and the array sizes were set to

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:

The constants *h*_{0}, Δ*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

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

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

where −1 < *r*_{u} < 1 and similarly for *υ*_{rand}.

### d. Parameter values

The constants *h*_{0} 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

The initial macroscale Reynolds number Re_{0} expresses the relative importance of viscous effects. In the present simulations, Re_{0} 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 Ri_{0} 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 *h*_{0}/2, Δ*u*/2, and Δ*ρ*/2, because this choice simplifies some of the mathematical expressions. The choice of *h*_{0}, Δ*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) = Ri_{0}, and increases without bound at large |*z*|. Ri_{0} < 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 Ri_{0}. In fact, instability at Ri_{0} > 1/4 requires *R* > 2.4 (Smyth and Peltier 1989).

To motivate the remaining parameter choices, we now explicitly compute the linear stability characteristics of the stratified shear flow described by (11)–(16) with Re_{0} = 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:

Here, *ϕ* may represent any perturbation field, *α* and *β* are streamwise and spanwise wavenumbers nondimensionalized by 2/*h*_{0}, and *σ* is a (possibly complex) exponential growth rate nondimensionalized by Δ*u*/*h*_{0}. 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 *h*_{0} are fixed so that

The ratio Ri_{0}/*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 Ri_{0}) to strong stratification concentrated in a thin layer (at large *R* and Ri_{0}), 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 (Ri_{0} = 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.

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 Ri_{0} 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 *L*_{x} 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 *h*_{0} was set to 0.1795 m. The (dimensional) wavelength was chosen as 2*π*/*α* × *h*_{0}/2 with *α* = 0.35 (cf. Fig. 3), so that *L*_{x} = 1.6116 m. Here *L*_{z} was chosen as *L*_{x}/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 *L*_{y} was chosen to be *L*_{x}/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 *L*_{y} 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 *h*_{0}/Δ*u* 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 *h*_{0}/Δ*u.*

## 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.)

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 Re_{0} = 1200, Pr = 1, Ri_{0} = 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 Ri_{0} ≤ 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 (Ri_{0} ≤ 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 (Ri_{0} = 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).

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.

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 *x*–*z* 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.

Next, we compare the *x*–*z* 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:

*K*_{3D} 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 Ri_{0} 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 *L*_{y}/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).

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.]

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 *y*–*z* plane reveal the counterrotating structures visible in the vorticity field (Fig. 9).

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 *d**u*/*dz*. We approximate *d**u*/*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} m^{2} 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.

The *x*–*z* structure of the three-dimensional motions is shown in Figs. 11b,d,f via contour diagrams of the kinetic energy *K*_{3D} superimposed on the spanwise-averaged density. The *x*–*z* 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), *K*_{3D} has grown by two orders of magnitude and has become much more symmetric. For the most part *K*_{3D} 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

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

where the integral is over the computational volume. The available potential energy is just the difference: *P*_{a} = *P*_{t} − *P*_{b}.

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

in which *ρ*|^{top}_{bottom} 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 *P*_{b} is monotonic. In contrast, *P*_{t} 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, *P*_{b} grows more rapidly once the flow has become fully turbulent, and *P*_{t} and *P*_{b} 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 *P*_{t} and *P*_{b} become nearly equal, the rate of potential energy increase is only negligibly larger than Φ.

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,

The dissipation rate is given explicitly by *ε* = *ν *Σ^{3}_{i,j=1} (*u*_{i,j} + *u*_{j,i})^{2}/2 and angle brackets denote a volume average over the computational domain. In this computation, Φ is subtracted from *dP*_{b}/*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 *dP*_{b}/*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 *dP*_{b}/*dt* − Φ. As turbulence decays, the two quantities decrease approximately in proportion to one another.

The approximate equality of *dP*_{b}/*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 *P*_{t} and *P*_{b} become equal and *dP*_{b}/*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.

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.

## 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 Ri_{0} (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.

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.

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 *L*_{y} 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* = ±*L*_{z}/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 Ri_{0} 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 Ri_{0}. The results of Staquet (2000) suggest that mixing efficiency in KH waves increases as Ri_{0} decreases, and therefore that the limit Ri_{0} → 0 is singular (because mixing efficiency must vanish in unstratified flow). Dependence of mixing physics on Ri_{0}, *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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## 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