Abstract

The Lagrangian properties of a high-resolution, three-dimensional, direct numerical simulation of Kelvin– Helmholtz (K–H) instability are examined with the goal of assessing the ability of Lagrangian measurements to determine rates and properties of ocean mixing events. The size and rotation rates of the two-dimensional K–H vortices are easily determined even by individual trajectories. Changes in density along individual trajectories unambiguously show diapycnal mixing. These changes are highly structured during the early phases of the instability but become more random once the flow becomes turbulent. Only 36 particles were tracked, which is not enough to usefully estimate volume-averaged fluxes from the average rates of temperature change. Similarly, time-and volume-averaged vertical advective flux can be estimated to only 20% accuracy. Despite the relatively low Reynolds number of the flow, Rλ ≈ 100, the dissipation rates of energy ɛ and density variance χ are correlated with the spectral levels of transverse velocity and density in an inertial subrange, as expected for high-Reynolds-number turbulence. The Kolmogorov constants are consistent with previous studies. This suggests that these inertial dissipation methods are the most promising techniques for making useful measurements of diapycnal mixing rates from practical Lagrangian floats because they converge rapidly and have a clear theoretical basis.

1. Introduction

Although it has been possible to measure mixing rates in the ocean interior for several decades (Gregg 1991), the methods typically require a ship and are labor intensive and, therefore, have sampled only a minute fraction of the ocean. Accordingly, we have only the crudest measurements of the geography of mixing (Gregg 1998). This situation is unlikely to improve rapidly without the development of robust autonomous techniques to measure ocean mixing rates.

One approach is to use instruments (Lagrangian floats) that are accurately Lagrangian on short time scales. D'Asaro (2003) reviews recent progress in developing such instruments. They measure temperature and their own depth (and sometimes much more), can make measurements for many months, and relay their data via satellite. In energetic quasi-homogeneous turbulence, such as that in the ocean mixed layer, the floats can gather a variety of standard turbulence statistics, including profiles of vertical kinetic energy, temperature variance, vertical heat flux, velocity, and temperature skewness, as well as less commonly measured Lagrangian quantities such as local heating rate and acceleration. Together, these can be used to diagnose terms in the heat and momentum budgets (D'Asaro 2004), relate the turbulence quantities to forcing (Steffen and D'Asaro 2002; D'Asaro 2001), and test the performance of numerical models of mixed-layer turbulence (Harcourt et al. 2002).

In the stratified ocean interior, the interpretation of float data is more difficult. Here, the turbulence is less energetic so that its dominant scales are often smaller than the size of the float (Lien et al. 1998). This problem will not be explored here. A second problem results from the highly intermittent nature of the turbulence. Mixing occurs as a series of events in which the flow becomes unstable, the instabilities grow, break down into patches of turbulence, and then decay. The flow is thus both spatially and temporally inhomogeneous and only intermittently turbulent. Stratification also allows for large wavelike components in the flow that may have large variances but only small fluxes (D'Asaro and Lien 2000).

The goal here is to investigate the ability of Lagrangian methods to diagnose and quantify the behavior of intermittent turbulence in a stratified flow. This is done by the direct numerical simulation of a three-dimensional Kelvin–Helmholz (K–H) mixing layer in which a modest number of exactly Lagrangian trajectories are computed. The simulations mimic those of Smyth et al. (2001). We briefly describe the model in section 2. Section 3 focuses on diagnosing the large-eddy properties of the flow and mixing estimates based on these. Section 4 focuses on estimating the dissipation rates of energy and density, ɛ and χ, from Lagrangian spectra. The results are summarized in section 5.

2. The numerical model

a. Setup

The analyses presented are based on a direct numerical simulation of a temporally evolving stratified mixing layer. A three-dimensional turbulent flow field is allowed to develop from an initial state comprised of a vertically varying, density-stratified shear flow with velocity profile

 
formula

and density profile

 
formula

where u0 = 0.0133 m s−1, h0 = 0.38 m, and δρ0 = 3.91 10−3 kg m−3. Here z is the vertical coordinate oriented positive upward. The parameters were chosen to match those in Smyth et al.'s (2001) simulation 3 or, equivalently, Smyth and Moum's (2000b) simulation R16P1. The initial state is characterized by the dimensionless parameters Re0 (=u0h0ν−1) ≈ 5000, Ri0(=gδρ h0 ρ−10u−20) = 0.08, and Pr = ν/κ = 1. We have taken ρ0 = 1027 kg m−3 as the reference density. The bulk shear time h0/u0 = 28.57 s. The buoyancy time 1/N = 100.86 s, where N2 = (gδρ)/ρ0h0.

As discussed in Smyth et al. (2001), the initial state is unstable to K–H instability at a streamwise horizontal wavelength equal to one-half the length of the computational domain Lx. Choosing Lx to be twice the instability wavelength allows the development of two billows associated with the primary instability and one pairing event. Smyth et al. (2001) found that this flow was insensitive to the spanwise scale of the computational domain for widths of at least Lx/8. Accordingly, we have used this scale as well. In dimensional values the flow is computed in a domain of size (Lx, Ly, Lz) equal to (5.24, 0.65, 2.62) m using a uniformly spaced grid of (512, 64, 256) points, which gives an approximately isotropic spatial mesh size Δ of about 0.01 m.

The underlying flow simulation consists of solving the Boussinesq Navier–Stokes equations for a density-stratified flow. We assume laterally periodic boundary conditions with free-slip rigid-lid conditions applied to the upper and lower boundaries. Though the description of the flow simulation coincides with that in Smyth et al. (2001) to this point, there are technical differences in the methodologies worth noting. Following Smyth and Moum (2000a), we perturb the ambient states (1) and (2) at horizontal wavelengths corresponding to both the primary and secondary (pairing) modes of instability. We also add random density perturbations with an rms magnitude of about 1% of δρ0 at each grid point within a centered window spanning the vertical extent of the ambient shear. Our perturbations differ slightly from those of Smyth and Moum (2000a) in that we perturb the density field rather than velocity. As expected, these differences are unimportant and both simulated flows progress from primary instability to pairing to three-dimensional collapse. A more significant difference is the accuracy of the discrete differentiation schemes employed in the solution procedure. Smyth et al. (2001) used Fourier representation in the lateral directions coupled with a second-order finite-difference scheme in the vertical. This mixed approach can have subtle but important consequences, with the relatively inaccurate finite-difference scheme significantly degrading the overall resolution properties of the solution. The simulations presented here were computed using a spectral approach in all three directions: Fourier in the lateral directions, and sine and cosine expansions in the vertical. This approach maximizes the effective resolution for a fixed number of equally spaced grid points. An indication of the difference in effective resolution is the appearance of secondary instabilities in the braids (see Fig. 1) that are not evident in the simulations of Smyth et al. (2001). Despite these minor differences, the overall evolution of the flow in the two simulations is not significantly different.

Fig. 1.

Examples of density contours at y = 0.32 m (the center) and positions of Lagrangian particles projected on the xz plane (“o”). Density is scaled between 0.002 (black) and −0.002 kg m−3 (white). Secondary instabilities, which take the form of two small billows on the main braid, are apparent at t = 1510 s

Fig. 1.

Examples of density contours at y = 0.32 m (the center) and positions of Lagrangian particles projected on the xz plane (“o”). Density is scaled between 0.002 (black) and −0.002 kg m−3 (white). Secondary instabilities, which take the form of two small billows on the main braid, are apparent at t = 1510 s

In addition to computing the Eulerian flow fields, we also compute 36 Lagrangian trajectories within the flow. Along these trajectories we sample the underlying flow to generate time series of flow properties. We have taken care to minimize the two sources of numerical error in computing Lagrangian trajectories in time-varying Eulerian fields so that the Lagrangian density equation /Dt = κ2ρ is satisfied. First, the fields are interpolated to points between mesh points by direct Fourier summation. This eliminates interpolation error in the sense that the interpolated values at arbitrary points are as good as the Eulerian values at the mesh points. We also use a higher-order time-stepping (fourth-order Runge–Kutta) method for computing the trajectories than we use for the Eulerian simulation (third-order Adams– Bashforth). The use of this highly accurate scheme was computationally expensive, which is why only 36 trajectories were computed.

b. Flow evolution

Figure 1 shows 12 snapshots of the simulated density field on the xz plane at y = 0.325 m (midplane in the spanwise direction). The flow can be characterized by an instability phase (approximately 0–2000 s) and a turbulent phase (2000–4128 s). The Smyth et al. (2001) simulations, which continue to approximately 10 000 s, also show a decaying turbulence phase. During instability two K–H billows form and grow until about 1141 s. Secondary instabilities of the braid are apparent between the two billows by 1510 s. The flow remains nearly two-dimensional in xz. By 1880 s the billows have begun to collapse and merge, and the flow is in transition into three-dimensional turbulence. Parameters of the turbulent phase can be drawn from Smyth and Moum (2000a,b). During this phase the bulk Richardson number gδρ0 h0/u20 increases from about 0.15 to about 0.22. The buoyancy Reynolds number ɛ/νN2 and the microscale Reynolds number Rλ = /ν (q2 is the kinetic energy, λ2 = 5νq2/ɛ) are about 100, just barely large enough for the dissipation scales to be approximately isotropic. The rms Thorpe scale, drawn from either 1D or 3D sorting, is about 40 cm. The Osmidov scale is about 15 cm. The Cox number is about 500. The mixing efficiency, computed as in Winters et al. (1994), drops from about 0.8 at 2000 s to less than 0.2 by the end of the simulation.

3. Lagrangian diagnosis of large-eddy properties

a. Trajectories and structure

Lagrangian data were computed for 36 infinitesimally small particles seeded uniformly along the diagonal from (x, y, z) = (0, 0, 1.0 m) to (Lx, Ly, 1.62 m), that is, within the initial region of high shear and stratification (Fig. 2a). Density and all three velocity components at both the instantaneous positions of particles (xf, yf, zf) as well as at the nearby positions below (xf, yf, zf − Δz) are sampled at an approximately 1-s interval.

Fig. 2.

Vertical profiles of (a) velocity and (b) density at the initial time (black solid curves) and the final time (dashed curves) of the model simulation, and (c) histogram of vertical positions of Lagrangian particles. Solid dots in (a) and (b) indicate the initial velocity and density of the 36 seeded particles. Open circles show the final velocity and density of the particles. Two horizontal lines on the right axes of (c) encompass 95% of vertical positions of Lagrangian particles

Fig. 2.

Vertical profiles of (a) velocity and (b) density at the initial time (black solid curves) and the final time (dashed curves) of the model simulation, and (c) histogram of vertical positions of Lagrangian particles. Solid dots in (a) and (b) indicate the initial velocity and density of the 36 seeded particles. Open circles show the final velocity and density of the particles. Two horizontal lines on the right axes of (c) encompass 95% of vertical positions of Lagrangian particles

The Lagrangian trajectory data (Figs. 3 and 4) clearly show the evolution of the flow and its dominant large eddies. The growth of the initial instability is apparent in the z-time trajectories of Fig. 3a; the transition to three-dimensional flow at about 1500 s is apparent in the y-time trajectories of Fig. 3b. The initial two-vortex structure and their pairing into a single vortex is clear in the upper two panels of Fig. 4. The persistence of this vortex to the end of the simulation is apparent in the bottom panels of Fig. 4. The large-eddy properties are also apparent in the individual trajectories (Fig. 5).

Fig. 3.

Time series of Lagrangian trajectories: (a) z time, (b) y time, (c) ρ time

Fig. 3.

Time series of Lagrangian trajectories: (a) z time, (b) y time, (c) ρ time

Fig. 4.

Lagrangian trajectories (xz) within selected time intervals. The dashed box indicates the width of the domain, 5.24 m, in the x direction

Fig. 4.

Lagrangian trajectories (xz) within selected time intervals. The dashed box indicates the width of the domain, 5.24 m, in the x direction

Fig. 5.

Selected individual xz Lagrangian trajectories. The horizontal line indicates the width of the domain in the x direction

Fig. 5.

Selected individual xz Lagrangian trajectories. The horizontal line indicates the width of the domain in the x direction

Notice that the Lagrangian trajectories do not uniformly sample the model volume because they are initially seeded within the central sheared region. Accordingly, during the initial vortex roll-up they are concentrated within the billows and do not sample the “braid” region between the vortices. The particles only approach a horizontally uniform distribution by the end of the simulation (Fig. 1).

This inhomogeneity in the Lagrangian sampling makes some quantitative comparisons between the Eulerian and Lagrangian measurements difficult. For example, the average vertical kinetic energy measured by the Lagrangian particles is 2.7 × 10−6 m2 s−2; the Eulerian average over the entire time and domain is 1.1 × 10−6 m2 s−2. A more uniform distribution of Lagrangian measurements would undoubtedly improve this comparison.

Comparisons with overturning scales are more useful because most of the overturning occurs in the well-sampled vortex. For example, the vortex has a diameter of about 1.1 m; assuming that the vortex is circular and rotates the preexisting stratification as a solid body with no mixing, the maximum overturning scale is 1.1 m and the rms displacement is 0.55 m. The rms z displacement of a trajectory during the turbulent phase is 0.34 ± 0.09 m. The rms Thorpe scale, 0.4 m, is between these two values.

b. Mixing rates

1) Indicators

A surprising result of these simulations, and those of Smyth et al. (2001), is that a large amount of mixing occurs during the two-dimensional instability period before 2000 s. This is apparent from the large volume of “gray” fluid that is created in the first few frames of Fig. 1. Similarly, the volume average of χ increases rapidly in the first 700 s of the simulation, is strongest between 800 and 1800 s, and decreases gradually afterward (Fig. 6). The clearest Lagrangian signal is the rapid change in ρ(t) from 500 to 1000 s (Fig. 3c). Note that ρ can only change because of mixing along a Lagrangian trajectory, so changes in ρ are a direct measure of mixing rates. In Fig. 3c the lightest particles get heavier and the heaviest ones get lighter, so that by 1500 s the range of ρ is decreased dramatically. This corresponds exactly to the creation of “gray” fluid from “white” and “black” fluids in Fig. 1.

Fig. 6.

Time series of ɛ (black) and χ (gray) averaged over the whole model volume (thin curves) compared with those averaged over the 36 Lagrangian trajectories (thick curves)

Fig. 6.

Time series of ɛ (black) and χ (gray) averaged over the whole model volume (thin curves) compared with those averaged over the 36 Lagrangian trajectories (thick curves)

2) Mass balance

A more formal measure of the creation of mixed fluid is the density-sorted depth coordinate z∗(ρ) (Winters et al. 1994), which measures the volume of fluid greater than a given density. The changes in z∗ can be used to infer the rate of mixing of fluid between density classes and the potential energy changes associated with this mixing (Smyth et al. 2001, their Fig. 9). Although the density trajectories in Fig. 3c clearly show the transfer of fluid between density classes, they do not sample the density field with sufficient density and uniformity to usefully estimate z∗(ρ) or its rate of change.

3) Density flux

A common measure of mixing is the advective flux ϕa = , where the overbar indicates a time and/or space average. The average advective flux ϕa is computed by averaging ′ over the whole model domain (Fig. 7), where ρ′ = ρρ0. The Lagrangian equivalent ϕLa is computed by averaging measurements along all Lagrangian trajectories. Both ϕa and ϕLa fluctuate in time with positive and negative signs and exhibit large fluctuations in the preturbulent period, although their temporal variations are rather different. The magnitude of ϕLa is generally greater than ϕa because particles were initially seeded in a region of large advective flux. Averaged over time, however, 〈ϕLa〉 is 5.61 × 10−8 m s−1 kg m−3 and 〈ϕa〉 is 7.32 × 10−8 m s−1 kg m−3. If we compute ϕa by averaging ϕa horizontally over the whole domain but vertically bounded by the particles' vertical positions, then the 〈ϕa〉 is 6.66 × 10−8 m s−1 kg m−3, closer to 〈ϕLa〉.

Fig. 7.

Time series of advective density flux computed from Eulerian data over the entire domain (thick gray) and from all Lagrangian data (thin solid)

Fig. 7.

Time series of advective density flux computed from Eulerian data over the entire domain (thick gray) and from all Lagrangian data (thin solid)

4) Density equation

Terms from the horizontally averaged advective–diffusive equation of density,

 
formula

where 〈 〉 represents horizontal averaging over the entire time period and D/Dt is the Lagrangian derivative, are shown in Fig. 8.

Fig. 8.

Profiles of (a) advective density flux computed from Eulerian (thick gray, E) and Lagrangian (thin solid, L) data. Vertical lines on top axis show depth averages. (b) ∂ϕa/∂z and ∂〈ρ〉/∂t; Eulerian estimates are gray and thick dashed; Lagrangian estimates are thin solid and circles. These cannot be computed at all depths. (c) Eulerian (thick gray) and Lagrangian (circles) estimates of 〈/Dt〉 = ν2ρ. Small circles show Lagrangian estimates that are off scale

Fig. 8.

Profiles of (a) advective density flux computed from Eulerian (thick gray, E) and Lagrangian (thin solid, L) data. Vertical lines on top axis show depth averages. (b) ∂ϕa/∂z and ∂〈ρ〉/∂t; Eulerian estimates are gray and thick dashed; Lagrangian estimates are thin solid and circles. These cannot be computed at all depths. (c) Eulerian (thick gray) and Lagrangian (circles) estimates of 〈/Dt〉 = ν2ρ. Small circles show Lagrangian estimates that are off scale

The vertical profile of 〈wp′〉 (Fig. 8a) computed from the Lagrangian trajectories shows the correct (Eulerian) shape but only one-half the true magnitude near the center of the model domain and a completely different shape elsewhere. The second and third terms in (3) balance almost exactly (Fig. 8b, gray and dashed lines). The Lagrangian estimate of ∂〈ρ〉/∂t (circles) is quite accurate in the center of the domain; the Lagrangian estimate of ∂〈wp′〉/∂z is accurate to about 40% in this region; both are poorly estimated outside of 0.2 of the center. The last term in (3) is small and poorly estimated by the floats.

4. Microscale flux estimates

a. Trends

The most commonly used methods to infer mixing rates within the stratified ocean use measurements of the dissipation rate of turbulent kinetic energy ɛ and the scalar (density) dissipation rate χ. These can be used to infer turbulent diffusivities via the Osborn (1980) method,

 
KργN−2
(4)

where γ is a “mixing efficiency” (see Smyth et al. 2001) and N is the buoyancy frequency, and the Osborn and Cox (1972) method,

 
formula

where g is the acceleration of gravity.

As noted above, the dissipation of density variance χ is largest during the two-dimensional phase of the flow (Fig. 6). In contrast, the volume average of ɛ remains small in the first 1200 s, increases rapidly, and reaches its maximum at 1800 s, then decreases very slowly, by less than a factor of 2, until the end of the simulation. Over the course of the simulation, most of χ results from the preturbulent billows, whereas most of the ɛ occurs in the turbulent phase.

The values of χ and ɛ along the trajectories were computed at the grid point closest to each particle position. These values are greater, on average, than the volume averages of ɛ and χ (Fig. 6). This most likely results from the initial seeding of the particles within the shear-stratified layer of the flow, where the turbulence levels are higher. The highly patchy distribution of χ and ɛ also affects these averages. For example, χ on the floats is lower than the Eulerian average during the second half of the simulation, when strong χ exists only in sparse regions with large strain, which are undersampled by the particles. During the same period, strong ɛ exists mostly within the interior region, which is well sampled by the particles.

b. Energy dissipation rate from Lagrangian acceleration

For scales smaller than those of the large eddies, Kolmogorov scaling predicts that velocity spectra will scale with ɛ, ν, the wavenumber. This leads to the well-known “−5/3” wavenumber spectrum for velocity. Corrsin (1963) first applied these same ideas to the Lagrangian spectrum of acceleration Φa(ω) as a function of frequency ω. For ω0ωων, where ω0 is a “large-eddy” frequency and ων = ɛ1/2ν−1/2 is the Kolmogorov frequency, an inertial subrange should exist in which Φa(ω) = βɛ. Lien and D'Asaro (2002) review existing measurements and find strong evidence for a Lagrangian inertial subrange. The value of β is between 1.8 and 2.2. Thus, in homogeneous turbulence with sufficiently high Reynolds number, the Lagrangian acceleration spectrum can be used to estimate ɛ with an accuracy of better than 10%. Here, we test how well ɛ can be estimated from acceleration spectra in a moderate-Reynolds-number stratified flow with large coherent structures.

Lagrangian velocity spectra were computed for all particles during the turbulent part of the simulation (2000–4128 s). During this period ɛ is nearly steady both averaged over the domain and at the particles (Fig. 6). The Kolmogorov frequency ων is about 0.03 s−1. The large-eddy frequency ω0 is obviously the rotation rate of the main vortex, about 0.005 s−1. The frequency inertial subrange is therefore less than 1 decade, that is, 0.005–0.03 s−1 wide, as would be expected for a flow with a low Reynolds number.

Figure 9 shows two sample acceleration spectra with different energy levels. At low frequencies, the vertical (W) spectrum is higher than the transverse (V) spectrum, particularly for the more energetic of the two spectra. This reflects two-dimensional (u, w) polarization of the dominant vortex; the W spectrum shows a contribution from the vortex, while the V spectrum does not. At high frequencies the two spectra are nearly identical, reflecting the isotropy of the small scales (Smyth and Moum 2000b). In these plots an inertial subrange would appear as a plateau with a spectral level of βɛ (see Lien et al. 1998 for examples). Rather than a plateau, only a broad maximum is seen. The inertial subrange occurs for frequencies between ω0 and ων. For the relatively low Reynolds number of this simulation, ων/ω0 is not very large and the inertial subrange covers only a short frequency range. Any inertial subrange should be more apparent in the V spectrum than in the W spectrum since the V spectrum is insensitive to the large vortex and thus should reflect the “universal” properties of turbulence more strongly.

Fig. 9.

Lagrangian acceleration spectra of two typical floats, one encountering strong turbulence (solid curve) and the other weaker turbulence (dashed curve). Both vertical acceleration spectra (light, W) and transverse acceleration spectra (heavier, V) are shown. Shading indicates 95% confidence intervals for fitted ɛ computed using Monte Carlo simulations. The heavy solid and dashed curves show the model Lagrangian spectra fitted to these spectra

Fig. 9.

Lagrangian acceleration spectra of two typical floats, one encountering strong turbulence (solid curve) and the other weaker turbulence (dashed curve). Both vertical acceleration spectra (light, W) and transverse acceleration spectra (heavier, V) are shown. Shading indicates 95% confidence intervals for fitted ɛ computed using Monte Carlo simulations. The heavy solid and dashed curves show the model Lagrangian spectra fitted to these spectra

We tested whether the maximum of the V acceleration spectrum scales like an inertial subrange and can therefore be used to estimate ɛ. First, the spectral shape in the viscous region must be determined. Figure 10a shows all spectra normalized for large-eddy scaling Φa = ɛF0(ω/ω0); Fig. 10b shows the spectra normalized for viscous scaling Φa = ɛF(ω/ων). The solid lines show an empirical spectrum,

 
formula

Here β = 1.8 is the Kolmogorov constant. The second factor is the form of the spectrum derived by Lien et al. (1998) for homogeneous turbulence excluding the viscous range. The third factor, empirically formed to fit the data in Fig. 10b, accounts for the falloff in the viscous range. For ων/ω0 = 10 the plots of ΦMa/ɛ show a clear inertial subrange. For ων/ω0 ≤ 5 only a maximum appears; for ων/ω0 = 2 the maximum becomes sharper and falls below the intertial subrange level.

Fig. 10.

(a) Large-eddy and (b) viscous scaling of normalized transverse (y) Lagrangian acceleration spectra. The thick solid curves are averages of computed Lagrangian spectra (dots); shading represents 95% confidence intervals of spectral estimates. The thin curves show model spectra for values of ων/ωo = 2, 5, and 10

Fig. 10.

(a) Large-eddy and (b) viscous scaling of normalized transverse (y) Lagrangian acceleration spectra. The thick solid curves are averages of computed Lagrangian spectra (dots); shading represents 95% confidence intervals of spectral estimates. The thin curves show model spectra for values of ων/ωo = 2, 5, and 10

The value of ε̂ is estimated for each particle by fitting (6) to its V spectrum. The V variance, and thus these fits, are dominated by the spectral values in the broad spectral peak, that is, in the region of the inertial subrange. Comparison of ε̂ with the true ɛ (Fig. 11) shows agreement within a factor of 2 over a range in ɛ of about 8. The average ratio of ε̂ and ɛ is 1.05 with 95% bootstrap confidence limits of 0.75 and 1.2.

Fig. 11.

Comparison of ɛ estimated by fitting Lagrangian spectra with the true ɛ averaged along the trajectory. Vertical lines indicate 95% confidence intervals for fitted ɛ computed using Monte Carlo simulations. Horizontal lines represent 95% confidence intervals for true ɛ computed by the bootstrapping method based on the values of ɛ at the points surrounding each point of the trajectory

Fig. 11.

Comparison of ɛ estimated by fitting Lagrangian spectra with the true ɛ averaged along the trajectory. Vertical lines indicate 95% confidence intervals for fitted ɛ computed using Monte Carlo simulations. Horizontal lines represent 95% confidence intervals for true ɛ computed by the bootstrapping method based on the values of ɛ at the points surrounding each point of the trajectory

c. Density variance dissipation rate from Lagrangian density spectra

The properties of Lagrangian density spectra are not well understood. Dimensional analysis and analogy with velocity spectra suggests that the Lagrangian spectrum of ρt = /Dt should be white in the inertial subrange, with a level proportional to χ (Yeung 2001), that is, Φρt = βθχ. Lien and D'Asaro (2004) show that this shape is consistent with (5) if βθ = 1/π. Lien et al. (2002), using data from the equatorial undercurrent, compute χ from Φρt, ɛ from Φa, and find (4) consistent with (5) for βθ = O(1). Here, we examine the shape of Φρt and test whether its spectral level scales with χ.

The analysis is similar to that for Φa. Spectra are computed for all trajectories between 2000 and 4128 s. Unlike ɛ, however, χ decreases by more than a decade during this time (Fig. 6), and so, compared to ɛ, χ is inhomogeneous and unstationary. The spectra (Fig. 12a) vary by more than 3 decades, 10−13 to 10−10 kg2 m−6 s−1 at low frequency. Normalizing by the known average χ along each trajectory (Fig. 12b) collapses the spreading of spectral levels in the inertial subrange by one-half. We estimate χ̂ from the spectral level of Φρt, using either the second spectral estimate (dark shading in Fig. 12a) or the average of the second and third estimates (lighter shading in Fig. 12a). The value of χ̂ is compared with the true χ in Fig. 13. The correlation coefficient is 0.86 with χ̂ ≈ 1/4χ, very close to the value of 1/π estimated theoretically by Lien and D'Asaro (2004).

Fig. 12.

(a) Spectra of /Dt from Lagrangian trajectories. Shaded regions show spectral values used in Fig. 13. (b) Same as (a) but normalized by χ along these trajectories

Fig. 12.

(a) Spectra of /Dt from Lagrangian trajectories. Shaded regions show spectral values used in Fig. 13. (b) Same as (a) but normalized by χ along these trajectories

Fig. 13.

Comparison of χ estimated from the second spectral value in Fig. 12 (dots) or the average of second and third spectral values (crosses) with true ɛ along the trajectory. Diagonal lines give values of 1/βρ

Fig. 13.

Comparison of χ estimated from the second spectral value in Fig. 12 (dots) or the average of second and third spectral values (crosses) with true ɛ along the trajectory. Diagonal lines give values of 1/βρ

5. Summary

The Lagrangian characteristics of a stratified shear evolving through Kelvin–Helmholtz instability were studied using direct numerical simulation. The simulation is nearly identical to that of Smyth and Moum (2000b), which was therefore used to provide background information on the flow. The trajectories of 36 particles were computed, and the ability of these Lagrangian data to describe and quantify the flow, especially during its turbulent phase, was studied. The major results are as follows.

  • The particle trajectories clearly show the initial instability, the evolution of the two initial vortices into a single spanwise vortex, the development of three-dimensionality and turbulence, and the persistence and structure of the vortex during the turbulent phase of the flow.

  • The Lagrangian trajectories in this simulation cannot accurately compute the mean turbulent kinetic energy of the flow because they are not uniformly distributed throughout the volume, but rather concentrated in the vortex. We expect that a sufficient number of uniformly or randomly distributed trajectories would be able to accurately estimate the kinetic energy.

  • The rms vertical displacement of the Lagrangian trajectories during the turbulent phase (0.34 m) is close to the rms Thorpe scale (0.4 m).

  • Individual trajectories cleanly show the size, rotation rate, and evolution of the large vortex.

  • The volume mean advective flux estimated from the Lagrangian trajectories is within 20% of the (correct) Eulerian value. Its vertical divergence ∂z〉 is accurately reproduced within the depth range where the particles are concentrated. The time variation and the vertical profile of the advective flux, however, are not well reproduced by the particles.

  • The evolution of density along the trajectories clearly shows the mixing of light and heavy water within the shear layer and its change over time. Despite this, there are not enough trajectories to usefully estimate fluxes from the mass balance of density.

  • Lagrangian spectra of acceleration are anisotropic at low frequencies, showing the dominance of the two-dimensional vortex, but isotropic at and above the Kolmogorov frequency. The spectrum of transverse acceleration, which is insensitive to the vortex, has a form close to that expected for homogeneous turbulence, but the spectrum of vertical acceleration does not. In the dissipative frequency range, ω > ων, normalized acceleration spectra show a nearly universal and isotropic form with a spectral slope of −5.

  • Although a clear inertial subrange for acceleration does not exist because the Reynolds number is too low (Rλ ≈ 100), the spectral level near the peak of the transverse acceleration spectrum scales with the turbulent kinetic energy dissipation rate ɛ consistent with the expected value of the Kolmogorov constant.

  • Lagrangian spectra of /Dt show a short region with a flat slope, consistent with the theoretically expected white slope in an inertial subrange. The spectral levels here are correlated with the χ, the dissipation rate of density variance, with a Kolmogorov constant (1/4) close to that (1/π) derived theoretically by Lien and D'Asaro (2004).

6. Discussion and conclusions

The major practical goal of this study was to assess the ability of Lagrangian floats to measure mixing rates within the stratified ocean. The shear instability simulation was thus used to approximate a typical ocean mixing event. Although it displays many features in common with such events, it is somewhat unsatisfactory in several ways. First, it is a very weak event. Although the shear and stratification are typical of the upper ocean, the dissipation and mixing rate is typical of weak open-ocean “Garrett–Munk” mixing. Practical Lagrangian floats used in the ocean are only useful at much higher mixing rates, where the overturning scales become comparable to the size of the floats. Second, a significant fraction of the mixing in the simulation occurs during the two-dimensional instability phase, during which there is no turbulence. This would appear to be a classic symptom of a low-Reynolds-number flow. Third, the simulation, unlike the ocean, has Pr = 1. These limitations are partially overcome by confining the bulk of the analysis to the second half of the simulation, which is more turbulent, and by understanding that the main questions posed are semiqualitative: which techniques for estimating mixing rates from Lagrangian data are most promising?

Lagrangian floats offer the possibility of measuring ocean mixing remotely and for long periods of time. For this application, methods that require several floats within the same mixing event are highly impractical. This eliminates mass balance methods [section 3b(2)], which also do not appear to work well even with the tens of floats used in this simulation. Accurate three-dimensional tracking is also difficult over long distances, while accurate measurements of depth and density are easily accomplished. Thus, although the vertical motion of the floats appears to correlate well with overturning scales, accurate discrimination between vertical advection by internal waves, which results in no mixing, and overturning, which results in mixing, would appear to require three-dimensional tracking, so that pictures similar to those in Fig. 4 could be constructed. This does not appear to be practical.

An array of Lagrangian floats could, in principle, be used as a flow visualization tool in the ocean. Figure 5 clearly shows that even individual trajectories can reveal the structure of the K–H instability. Figure 4 shows that an array of floats within the same feature does an excellent job of visualizing the flow. This should be no surprise because such techniques have been used in laboratory flows for many decades. These methods, however, present many practical difficulties, including the tendency for floats to rapidly disperse while waiting for a feature to appear, the relatively high cost of individual floats, and, as discussed above, the difficulties of long-term three-dimensional tracking.

The two microscale methods for estimating mixing appear highly promising for long-term measurements. Each is widely used (Moum and Caldwell 1995) by the microstructure community and has a strong theoretical basis. We would expect that floats deployed for long periods in a stratified fluid would randomly sample mixing events and, over time, obtain accurate average values of the spectral levels of acceleration and temperature at inertial subrange frequencies. The main issue, then, is whether the floats can estimate ɛ and/or χ from these levels.

For high-Reynolds-number unstratified turbulence it is clear that an acceleration inertial subrange exists and that its level is proportional to ɛ (Lien and D'Asaro 2002). The flow simulated here is stratified and has moderate Reynolds number with only a hint of an inertial subrange. The level of the inertial subrange is correlated with ɛ as expected for the transverse component, but not for the vertical component, which is influenced by the energetic vortex. This shows that at moderate Reynolds number the spectral width of the inertial subrange is short and its level can be contaminated by coherent structures or, equivalently, by the form of the initial instability. For measurements by practical floats for which the Reynolds number will be much higher, this may not be an important issue. Sreenivasan (1995), for example, reviews a variety of turbulent flows and finds that for Rλ > 1000 all velocity components show a robust −5/3 inertial subrange.

Similarly, although the evidence for a Lagrangian inertial subrange for temperature is limited, the spectral level of the temperature derivative in this simulation appears to scale with χ, as predicted by Lien and D'Asaro (2004). These measurements are easy to make and have a low noise level. This, plus the measurements of Lien et al. (2002) and the theoretical discussion of Lien and D'Asaro (2004), suggests that a better understanding of Lagrangian temperature spectra could yield a robust mixing measurement.

Direct estimation of advective flux 〈〉 may also be practical. However, because this quantity has fluctuations much larger than the mean, a large amount of averaging is required to make it converge. A single float would probably have to be averaged over many dozens of mixing events to obtain a stable value of 〈〉. This may be possible in long records.

These results for Lagrangian measurements are qualitatively similar to those well known for undersampled Eulerian measurements of stratified turbulence, such as those made from oceanic profilers. Microstructure measurements yield reliable and rapidly converging measurements of fluxes. Covariance fluxes and fluxes from scalar budgets converge slowly so that it is often difficult to make sufficient measurements to obtain useful results.

Acknowledgments

This work was supported by NSF Grant OCE9617671. Simulations were done on the CRAY T3E at the Texax Advanced Computing Center under a National Partnership for Advanced Computational Infrastructure (NPACI) resource allocation.

REFERENCES

REFERENCES
Corrsin
,
S.
,
1963
:
Estimates of the relations between Eulerian and Lagrangian scales in large Reynold's number turbulence.
J. Atmos. Sci
,
20
,
115
199
.
D'Asaro
,
E. A.
,
2001
:
Turbulence intensity in the ocean mixed layer.
J. Phys. Oceanogr
,
31
,
3530
3537
.
D'Asaro
,
E. A.
,
2003
:
Performance of autonomous Lagrangian floats.
J. Atmos. Oceanic Technol
,
20
,
896
911
.
D'Asaro
,
E. A.
,
2004
:
Air–sea heat flux measurements from nearly neutrally buoyant floats.
J. Atmos. Oceanic Technol., in press
.
D'Asaro
,
E. A.
, and
R. C.
Lien
,
2000
:
Lagrangian measurements of waves and turbulence in stratified flows.
J. Phys. Oceanogr
,
30
,
641
655
.
Gregg
,
M.
,
1991
:
The study of mixing in the ocean: A brief history.
Oceanography
,
4
,
39
45
.
Gregg
,
M.
,
1998
:
Estimation and geography of diapycnal mixing in the stratified ocean.
Physical Processes in Lakes and Oceans, J. Imberger, Ed., Coastal and Estuarine Studies, Vol. 54, Amer. Geophys. Union, 305–338
.
Harcourt
,
R. R.
,
E. L.
Steffen
,
R. W.
Garwood
, and
E. A.
D'Asaro
,
2002
:
Fully Lagrangian floats in Labrador Sea deep convection: Comparison of numerical and experimental results.
J. Phys. Oceanogr
,
32
,
493
510
.
Lien
,
R. C.
, and
E. A.
D'Asaro
,
2002
:
The Kolmogorov constant for the Lagrangian velocity spectrum and structure function.
Phys. Fluids
,
14
,
4456
4459
.
Lien
,
R. C.
, and
E. A.
D'Asaro
,
2004
:
Lagrangian spectra and diapycnal mixing in stratified flow.
J. Phys. Oceanogr
,
34
,
978
984
.
Lien
,
R. C.
,
E. A.
D'Asaro
, and
G. T.
Dairiki
,
1998
:
Lagrangian frequency spectra of vertical velocity and vorticity in high-Reynolds number oceanic turbulence.
J. Fluid Mech
,
362
,
177
198
.
Lien
,
R. C.
,
E. A.
D'Asaro
, and
M. J.
McPhaden
,
2002
:
Internal waves and turbulence in the upper central equatorial Pacific: Lagrangian and Eulerian observations.
J. Phys. Oceanogr
,
32
,
2619
2639
.
Moum
,
J. N.
, and
D. R.
Caldwell
,
1995
:
Turbulence and mixing in the ocean.
Rev. Geophys
,
33
,
(Suppl. S),
.
1385
1394
.
Osborn
,
T. R.
,
1980
:
Estimates of the local rate of vertical diffusion from dissipation measurements.
J. Phys. Oceanogr
,
10
,
83
89
.
Osborn
,
T. R.
, and
C. S.
Cox
,
1972
:
Oceanic fine structure.
Geophys. Fluid Dyn
,
3
,
321
345
.
Smyth
,
W. D.
, and
J. N.
Moum
,
2000a
:
Anisotropy of turbulence in stably stratified mixing layers.
Phys. Fluids
,
12
,
1343
1362
.
Smyth
,
W. D.
, and
J. N.
Moum
,
2000b
:
Length scales of turbulence in stably stratified mixing layers.
Phys. Fluids
,
12
,
1327
1342
.
Smyth
,
W. D.
,
J. N.
Moum
, and
D. R.
Caldwell
,
2001
:
The efficiency of mixing in turbulent patches: Inferences from direct simulations and microstructure observations.
J. Phys. Oceanogr
,
31
,
1969
1992
.
Sreenivasan
,
K.
,
1995
:
On the universality of the Kolmogorov constant.
Phys. Fluids
,
7
,
2788
2784
.
Steffen
,
E.
, and
E.
D'Asaro
,
2002
:
Deep convection in the Labrador Sea as observed by Lagrangian floats.
J. Phys. Oceanogr
,
32
,
475
492
.
Winters
,
K. B.
,
P. N.
Lombard
,
J. J.
Riley
, and
E. A.
D'Asaro
,
1994
:
Available potential energy and mixing in stratified fluids.
J. Fluid Mech
,
289
,
115
128
.
Yeung
,
P.
,
2001
:
Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations.
J. Fluid Mech
,
427
,
241
274
.

Footnotes

Corresponding author address: Eric A. D'Asaro, Applied Physics Laboratory, University of Washington, 1013 NE 40th St., Seattle, WA 98105. Email: dasaro@apl.washington.edu