Документ взят из кэша поисковой машины. Адрес оригинального документа : http://ocean.phys.msu.ru/courses/stathyd/2005%20Smyth%20et%20al.%2C%20Differential%20diffusion%20in%20breaking%20Kelvin-Helmholtz%20billows.pdf Дата изменения: Wed Nov 11 15:01:42 2009 Дата индексирования: Mon Oct 1 19:55:43 2012 Кодировка:

S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

Differential diffusion in breaking Kelvin-Helmholtz billows
W.D. Smyth , J.D. Nash and J.N. Moum College of Oceanic and Atmospheric Sciences Oregon State University, Corvallis OR Submitted to J. Phys. Oceanogr. July 20, 2004 Revised: November 22, 2004
Corresponding author. Telephone: (541) 737-3029; email: smyth@coas.oregonstate.edu


1. Introduction
The density of seawater is controlled by a pair of scalar quantities, temperature and salinity, whose molecular diffusivities differ by two orders of magnitude. Despite this difference, we customarily assume that the turbulent diffusivities of temperature and salinity are the same. This assumption is grounded in the classical theory of stationary turbulence in the limit of infinite Reynolds number (e.g. Corrsin, 1951). However, much of the ocean interior is mixed by turbulent events for which the Reynolds number is decidedly finite (e.g. Moum, 1996b) and the turbulence is nonstationary. In a mixing event of finite duration, vertically displaced fluid parcels may return to an equilibrium configuration after mixing only partially with the surrounding fluid. The lower the molecular diffusivity, the greater the tendency for incomplete mixing. The large difference between the molecular diffusivities of heat and salt therefore suggests that heat and salt could mix differently in turbulent events of finite duration, i.e. that turbulent seawater may exhibit differential diffusivity. In the present study, we assess the potential for differential diffusivity via direct numerical simulations (DNS) of turbulent Kelvin-Helmholtz (KH) billows.

Abstract

Direct numerical simulations are used to compare turbulent diffusivities of heat and salt during the growth and collapse of Kelvin-Helmholtz billows. The ratio of diffusivities is obtained as a function of buoyancy Reynolds number Reb and of the density ratio R (the ratio of the contributions of heat and salt to the density stratification). The diffusivity ratio is generally less than unity (heat is mixed more effectively than salt), but it approaches unity with increasing Reb and also with increasing R . Instantaneous diffusivity ratios near unity are achieved during the most turbulent phase of the event even when Reb is small; much of the Reb dependence results from the fact that, at higher Reb , the diffusivity ratio remains close to unity for a longer time after the turbulence decays. An explanation for this is proposed in terms of the Batchelor scaling for scalar fields. Results are interpreted in terms of the dynamics of turbulent Kelvin-Helmholtz billows, and are compared in detail with previous studies of differential diffusion in numerical, laboratory and observational contexts. The overall picture suggests that the diffusivities Several large-scale modeling studies (e.g. Gargett and become approximately equal when Reb exceeds O(102 ). Holloway, 1992; Merryfield et al., 1999) have revealed The effect of R is significant only when Reb is less than that a difference in the assumed diffusivities of heat and this value. salt can lead to significant differences in computed large scale circulation, so the issue is potentially important 1


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

for the development of accurate turbulence parameteriza- of the density ratio and the important limit of zero net tions. stratification. A review of the subject has been provided Differential diffusion was first demonstrated in the lab- by Gargett (2003). oratory experiments of Turner (1968), who measured entrainment fluxes in a fluid where turbulence was generated by an oscillating grid. The working fluid was stratified by either temperature or salinity, but not by both. A significant difference in turbulent diffusivities was evident. Altman and Gargett (1990) repeated Turner's experiments, this time using thermal and saline stratification simultaneously. Like Turner, they found diffusivity ratios significantly different from unity. Individual entrainment rates were independent of the presence of the other density component, i.e. no dependence on the density ratio was detected. In the laboratory experiments of Jackson and Rehmann (2003), a fluid stratified by both salinity and temperature was stirred by oscillating rods, with special care taken to insulate the boundaries against heat loss. A distinct dependence on the buoyancy Reynolds number (defined below) was identified. Hebert and Ruddick (2003) measured differential diffusion of dynamically passive chemical dyes in breaking internal gravity waves, and again found a dependence on the buoyancy Reynolds number. Nash and Moum (2002) have made a similar assessment using in situ measurements of ocean microstructure. Statistical analysis of many turbulent events indicated a tendency for heat to diffuse more rapidly than salt, but the ratio of diffusivities was within experimental error of unity. No dependence upon the buoyancy Reynolds number was evident. The resolution of weakly diffusive scalars in a turbulent flow presents an extreme challenge for DNS. The first study to attempt this was Merryfield et al. (1998), in which flow was restricted to two dimensions to save memory. Those simulations were successful in detecting differential diffusion and they served as an important prelude to the first fully three-dimensional numerical realizations of the phenomenon, those of Gargett et al. (2003; hereafter GMH). To facilitate simulation in three dimensions, the diffusivity of salt was artificially increased (as it has been in all subsequent DNS studies including the present work). The results of GMH have recently been extended by Merryfield (2004, hereafter M04) to include variation 2 Here, we assess the potential for differential diffusion in turbulent KH billows. We do so using DNS of shear flows stratified by both heat and salt. KH billows have proven to be a useful model for shear-driven overturns observed in the ocean. Direct observations by Woods (1968) showed billows forming on the crests of larger-scale internal waves. Vivid images of KH-like billows have been obtained via echosounder in flow over topography (e.g. Seim and Gregg, 1994; Farmer and Armi, 1999) and in large amplitude internal waves (e.g. Moum et al., 2003). Smyth et al. (2001) have compared turbulence statistics from DNS of KH billows with measurements of turbulent events in the thermocline, and found that the two are statistically indistinguishable (except for the generally lower Reynolds numbers of the simulated flows, which reflects the limitations of existing computer technology, not of the KH model). Given this evidence for the importance of KH-like dynamics in ocean mixing events, we are motivated to learn whether, and if so under what conditions, turbulent KH billows exhibit differential diffusion. Section 2 describes the numerical model used for the simulations. A general overview of the KH life cycle as realized in these experiments is given in section 3. In section 4, we describe the scalar fields in terms of gradient spectra, and compare the results with both the ocean observations of Nash and Moum (2002) and the theoretical spectrum due to Kraichnan (1968). The main results are in section 5, where potential energy components, scalar variances and turbulent diffusivities for the two scalars are examined. In section 6, results are described in the context of previous work. A summary is given in section 7.

2. Methodology
a. The mathematical model
The mathematical model is based on the field equations for nonrotating, incompressible flow in the Boussinesq limit, together with advection-diffusion equations for the


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

two scalars, v iz .: ui t uj xj T t S t = = -u =0 T + T xj S = -uj + S xj 0 + T + S . = -u
j 2 2 j

ui 1 p - 0 - -g i3 + xj 0 xi 0

2

u

i



T



S

(1) The vector ui contains the components of the velocity field and p and represent pressure and density, respectively. The constant 0 is a reference density from which deviations are assumed to be small (so that the Boussinesq approximation applies and the equation of state is linear). Accordingly, the thermal and saline contributions to the density anomaly - 0 are represented by T = -0 (T - T0 ) and S = 0 (S - S0 ), where T0 and S0 refer to the reference state and and are the (constant) expansion and contraction coefficients for heat and salt in water. The molecular diffusivities of heat and salt in water are represented by the constants T and S . The constants and g represent kinematic viscosity and gravitational acceleration. The field equations (1) are solved in the computational box 0 x Lx , 0 y Ly , 0 z Lz . Boundary conditions are periodic in the horizontal directions, i.e. f (x + Lx , y , z , t) = f (x, y , z , t) = f (x, y + Ly , z , t), (2) where f represents any field variable. At the upper and lower boundaries z = 0, Lz , vertical velocity and vertical fluxes of heat, salt and horizontal momentum are required to vanish.

The Winters model has been extended for use in ocean DNS via the addition of a second active scalar, here representing salinity. The second scalar is resolved on a fine grid with spacing equal to one half the spacing used to resolve the other fields (as was done by GMH). Interpolations and decimations between grids are accomplished using Fourier transforms. Aliasing errors are reduced by applying to both grids at every time step an isotropic filter having a cosine-bell shape that decreases gradually from amplitude 1 to 0.6 over the range 0.8 to 1 times the Nyquist wavenumber. This gradual decrease minimizes the effect of dealiasing on the resolved fields. The multiple grid approach described above allows the efficient resolution of weakly diffusive scalars such as temperature and salinity in seawater. The memory requirement is about 1/3 of that required if all fields are resolved on the same grid. It is possible to increase the difference in resolution between the coarse and fine grids, but further increases yield only small improvements in efficiency.

c. Initial conditions and parameter values
For the present experiments, the initial conditions describe a pair of water masses separated by a horizontal transition layer: T S u =- =- = tanh u T S z - Lz /2 h0 . (3)

b. Numerical methods
The numerical code is an extension of that described by Winters et al. (2003). It uses Fourier pseudospectral discretization in all three dimensions. Time stepping is via the third-order Adams-Bashforth operator, with timestep determined by a Courant-Friedrichs-Lewy stability condition. Viscous and diffusive terms are integrated exactly. MPI routines are used for parallelization. 3

Here, h0 is the initial half-depth of the transition layer, and u is the half-velocity difference. T and S are the absolute values of the half-differences of the density components T and S , respectively, so that the absolute halfdifference of density across the layer is = T + S . With these choices, the initial stratification is both statically and diffusively stable. Dynamic (shear) instability depends on the relative values of h0 , u, T and S as discussed below. The horizontal periodicity intervals were determined according to the fastest-growing modes of linear theory. The domain length Lx was generally twice the wavelength of the fastest-growing KH mode, though a single wavelength was used for some experiments. For the profiles (3) with the parameter values used here, the fastest growing wavelength is closely approximated by F GM =


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

h0 в 2 /0.44. The domain width Ly was F GM /2, which choices of four nondimensional parameters: is approximately three times the spanwise wavelength of g h0 ; S c = /S ; Ri0 = the fastest-growing three-dimensional instability of KH 0 u 2 billows in air as described by Klaassen and Peltier (1991). uh0 T (Note that this wavelength is partly controlled by diffuRe0 = ; R = . (4) S sion, so we expect it to be smaller in seawater.) The Schmidt number, S c, for salt in seawater ranges beIn addition to the profiles described above, the initween 700 and 1000. In order to attain a significant level tial conditions included a two-part perturbation designed of turbulence in the computed flows, we have reduced this to efficiently stimulate both the KH mode and its secvalue to 50. Equivalently, one may express saline diffuondary instabilities. First, disturbances proportional to the sivity in terms of the inverse Lewis number, = S /T , fastest-growing KH mode and the KH mode with twice which is of order 10-2 in seawater but is 0.14 in these that wavelength were added. The amplitude of the fastestsimulations. (GMH used a similar value: S c = 70, or growing mode was chosen so that its maximum vertical =0.1.) Even with this compromise, "salinity" still difdisplacement was 0.2h0 . The maximum vertical displacefuses an order of magnitude more slowly than does heat, ment associated with the subharmonic mode was 0.1h0 . so the effects of the different molecular diffusivities ought These amplitudes are large enough to efficiently stimulate to be evident, although those effects are likely to be underprimary and subharmonic modes, yet small enough to be estimated. For simplicity we will refer to the scalar corwell described by linear perturbation theory. The phases responding to the density S as "salinity", even though of the primary and subharmonic modes were chosen to it actually represents a fictitious solute that diffuses more induce pairing at the streamwise boundary of the (perirapidly than sea salt. odic) computational domain, so that the inner core would The intensity of turbulence attained in the stratified be easily visible in volume renderings (e.g. figure 1). Secshear layer (3) is governed mainly by the initial Richardond, a random velocity field was added in order to excite son and Reynolds numbers, Ri0 and Re0 . The primary three-dimensional motions. At each point in space, the KH mode is inviscidly unstable provided that Ri0 < 1/4 three components of the velocity increment were chosen (Miles, 1961; Howard, 1961); for the present simulations, from a list of random numbers whose probability distriRi0 was in the range 0.10-0.12. The initial Reynolds numbution was uniform between the limits ±0.1u. During ber controls the range of scales in the resulting flow. A the first time step, the random motions were automatically standard compromise in DNS of geophysical flows, occamade solenoidal by the pressure gradient force. sioned by limitations of computer technology, is that the The computations were done using MKS units. To Reynolds number cannot normally be made as large as represent flow in terrestrial oceans, the gravitation ac- one would like. In this case, the slow diffusion of salinity celeration, characteristic density, molecular viscosity and requires that Re be set to 300 or smaller. Initial Richard0 thermal diffusivity were set to g = 9.81m/s2 , 0 = son and Reynolds numbers in this range lead to turbu1027k g /m3 , = 1.0 в 10-6 m2 /s, and T = 1.43 в lent patches whose intensity (as measured by the buoy10-7 m2 /s, respectively. Note that the choices of and ancy Reynolds number to be defined below) is within, but T correspond to a Prandtl number P r = /T of 7, near the weak end of, the range observed in thermocline a typical value for seawater. Approximate correspon- patches (Smyth et al., 2001). dence to a typical turbulent patch in the thermocline was The relative importance of heat and salt in determining achieved by setting the initial turnover time for the shear the initial density stratification is expressed by the denlayer TS = h0 /u0 to the value 28.28s. With this choice, sity ratio R . There are several conventions in current use KH billows were found to grow and decay over a time for defining R . With the definition given in (4), R is span of 1-3 hours. (Results can be converted to any other positive when both thermal and saline components of the time scale as necessary.) stratification are stable. Turbulent patches in the thermoThe remaining parameter values were determined via cline typically exhibit values of R between 0.2 and 5. 4


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

Parameter Re0 R Ri0 h0 u Lx Ly Lz Nx Ny Nz

Unit

10-3 m 10-3 ms-1 10-3 k g m-3 m m m

1 300 1.0 0.10 92.0 3.32 1.21 2.62 0.65 0.88 512 128 192

2 240 0.2 0.10 82.4 2.96 1.08 1.17 0.59 0.78 192 96 128

3 240 1.0 0.10 82.4 2.96 1.08 1.17 0.59 0.78 192 96 128

4 240 5.0 0.10 82.4 2.96 1.08 1.17 0.59 0.78 192 96 128

5 200 0.2 0.10 75.0 2.66 0.984 2.15 0.54 0.72 384 96 128

6 240 1.0 0.12 82.4 2.91 1.29 2.34 0.59 0.78 384 96 128

7 180 0.2 0.12 71.3 2.52 1.12 1.02 0.51 0.68 192 96 128

8 180 1.0 0.12 71.3 2.52 1.12 1.02 0.51 0.68 192 96 128

9 180 5.0 0.12 71.3 2.52 1.12 1.02 0.51 0.68 192 96 128

10 100 1.0 0.10 53.2 1.88 0.696 1.52 0.38 0.51 256 64 96

Table 1: Parameters for numerical simulations. Re0 and Ri0 represent the initial Reynolds and Richardson numbers, and R is the density ratio. The variable h0 is the half-thickness of the initial shear layer. The half-changes in horizontal velocity and net density are u and . Lx , Ly and Lz are the domain dimensions in the streamwise, cross-stream and vertical directions, respectively, and Nx , Ny and Nz are the corresponding array sizes. The dimensions of the fine array are 2Nx , 2Ny and 2Nz . For all simulations, P r = 7 and S c = 50.

Choices for the parameter ble 1. Most of the analysis remaining runs are included hensive view of the factors sion.

values are summarized in ta- dimensional components: will focus on runs 1-4; the 1 K2d (t) = to provided a more compre0 u governing differential diffu1 K3d (t) = 0 u

2

u u

2d

·u ·u

2d V

; . (6)

2

3d

3d V

3. Overview of flow evolution

The velocity fields associated with two- and threedimensional motions are

u2d (x, z , t) = u y - u xy ; The growth, breaking and decay of the KH billow in run 1 is illustrated in figure 1 via volume renderings of the S u3d (x, y , z , t) = u - u y , (7) field at selected times, and in figure 2 via the evolution of where subscripts on the angle brackets indicate spatial three energy reservoirs that we now define. The potential energy is given in nondimensional form averages over the specified dimensions. The velocity field u2d describes the primary KH billows, and other by large-scale, wavelike motions, while u3d is associated g P (t) = (z - Lz /2) (T + S ) V - P0 , (5) with longitudinal secondary instabilities (e.g. Klaassen 0 u 2 and Peltier, 1991) and turbulence. in which angle brackets indicate a volume average over Figure 1a shows the salinity field from run 1 at t = 0. the computational domain V , and P0 is the potential en- The transition layer was horizontal except for the smallergy of the initial profiles (3). P (t) evolves in response amplitude linear eigenfunction and the random noise field to both reversible and irreversible processes. Irreversible (figure 1a). Subsequently, both the potential and twopotential energy changes will be examined in section 5. dimensional kinetic energy fields showed rapid growth The kinetic energy is partitioned into two- and three- (figure 2, solid and dashed curves). By t = 610s, the 5


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

0.5

y

0 0.5

y

0

z
0.5

0

(a) t=0

) t= 2150s ((ff) t=2150s. .

z
0.5

0

(b) t=610s.

(g) t=2532s.

z
0.5

0

(c) t=1018s.

(h) t=2940s.

z
0.5

0

(d) t=1425s.

(i) t=3834s.

z
0.5

0

0

1

x

2

(e) t=1675s.

0

1

x

2

(j) t=5459s.

Figure 1: Evolution of the salinity field S for run 1. Values colored range from -0.4S (red) to 0.4S (dark blue). Values outside
this range are transparent. Times are as marked; note that the interval between frames is longer in the later part of the life cycle.

6


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

primary KH billows had rolled up and were approaching their maximum amplitude (figure 1b). In figure 1c, the 0.04 braid separating the billows at the center of the computational domain is visibly longer than that crossing the peP riodic streamwise boundary. This corresponds to a merging of the primary billows across the periodic boundary K3d due to the subharmonic pairing instability (e.g. Collins 0.02 and Maslowe, 1988). The merging process was nearly complete at t = 1425s (figure 1d). Also visible at this time was the emergence of three-dimensionality in the K 2d cores associated with the secondary instability described by Klaassen and Peltier (1985a,b, 1991). Four spanwise 0 wavelengths of the Klaassen-Peltier (hereafter KP) mode 0 1000 2000 3000 4000 5000 6000 t [s] are visible near the right-hand side of figure 1d. This instability was also manifested in rapid growth of the threeFigure 2: Selected energy reservoirs for run 1: potential energy dimensional kinetic energy (figure 2, dotted curve). Beyond this time, K2d decreased sharply. Most of this decrease was transferred to the mean flow as the quasielliptical billow core rotated to a more nearly horizontal orientation (figure 1e). The potential energy continued to grow for a short time after this due to the rollup of streamwise vortices associated with the KP instability (figures 1d,e); however, it too exhibited a rapid decrease around t = 2000s that coincided with rapid growth of three-dimensional structure (figure 1f, dotted curve in figure 2). This phase is referred to as the "breaking" of the KH billow.
(solid), kinetic energy of two-dimensional flow (dashed) and kinetic energy of three-dimensional flow (dotted). All energies are nondimensionalized by 0 u2 as described in the text. Potential energy is shown minus its initial value. Letters above the figure indicate times shown in figure 1. The numeral 6 indicates the time shown in figure 6.

sult of this thickening, the minimum Richardson number had increased to a value greater than 1/4, and the flow was therefore dynamically stable. This irreversible thickening of the transition layer is evident in figure 2 as a permanent The breaking billow cores ejected jets of turbulent fluid increase in potential energy after the disturbance kinetic horizontally toward the center of the domain (figure 1f), energies have decayed. where they engulfed the intervening braid. Figure 1g The flow evolution in runs 2-4 (figure 3) was simpler shows a second pair of jets being ejected from the tur- due to reduced Reynolds number and the suppression of bulent core. This ejection coincided with a second rapid pairing. The growth of the KP mode was in general more decrease in potential energy as the billow rotated again rapid because it did not compete for energy with the pairinto the horizontal orientation. The meeting of the second ing mode (also see Metcalfe et al., 1987). The growth rate pair of jets at the domain center (figure 1h) induced an of the primary KH instability was independent of R , as intense burst of turbulence. Shortly after this, turbulence is evident from the initial evolution of P and K (figure 3d began to decay under the influence of viscosity, as shown 3a,b). In contrast, the initial growth rate of the KP mode by the rapid decrease in both components of the kinetic was a strong function of R , as shown by the divergence energy (figure 2). Due to its low diffusivity, the "salinity" of the curves in figure 3c near t = 1500s. This variation field retained significant small scale structure even in the with R appears to be related to the P r dependence of the late stages of turbulence decay (figure 1i,j). growth rate of the KP mode described by Klaassen and Ultimately, the decay process left behind a sheared, Peltier (1985a). When R < 1 (dashed curve), the dentwo-layer flow similar to the initial condition, except that sity is dominated by the slowly-diffusing salinity compothe transition region had thickened due to mixing. As a re- nents, and the density gradients that drive convection are 7


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

(a)
0.01

2.5 1 2 3

P

0

(b) K2d
0.01

2
2: R=0.2 3: R=1 4: R=5

4

h/h

0

1.5
0.01

(c)

K3d

1 0
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

1000

2000

3000

4000

5000

time [s]

t [s]

Figure 3: Selected energy reservoirs for runs 3 (dashed),
4 (solid) and 5 (dotted): potential energy (a), kinetic energy of two-dimensional flow (b) and kinetic energy of threedimensional flow (c). All energies are nondimensionalized by 0 u2 as described in the text.

Figure 4: Evolution of the transition layer depth for simulations
1-4.

three-dimensional motions ultimately became strongest in the temperature-dominated case. We conclude this overview of KH breaking and turbulence with an examination of energy dissipation via viscous friction. Since a substantial fraction of our computational domain was occupied by laminar flow above and below the mixing layer, higher-order statistics such as the kinetic energy dissipation rate, when computed using simple volume averages over the domain, are not representative of the turbulent region. Instead, we take advantage of the fact that the turbulent layer coincides roughly with the transition layer identified previously, and is therefore delineated effectively by two isosurfaces of the total density field T + S . We choose isosurfaces upon which the density had the values ± tanh (1)(). The subvolume enclosed by these surfaces is denoted VT . At t = 0, the mean half-thickness of VT (denoted h(t)) was equal to h0 , the initial half-thickness of the transition layer. Averages over VT contain very little contribution from the laminar regions. As each simulation progressed, h(t) increased monotonically as a result of the irreversible mixing of density (figure 4). The degree of thickening was greatest in cases where mixing was most vigorous. It was this thickening that caused the increase of the bulk Richardson number to a stable value and hence the ultimate decay of the turbulence. 8

therefore sharper. In the R > 1 case (dotted curve), the converse was true: density was dominated by the rapidly-diffusing temperature field. Temperature dominance also caused the damping action of buoyancy on the primary KH billow to be reduced slightly, as shown by the increased amplitude and duration of the peaks in potential and two-dimensional kinetic energy (figures 3a,b, t 800 - 1400s). As in run 1, the breaking billows transferred much of their energy to the growing three-dimensional mode. This transfer occurred in two stages. In the cases with R 1 (solid and dotted curves), the second stage was considerably longer than the first and resulted in considerably greater growth in K3d . This is because (a) the primary billows lost less energy to three-dimensional motions during the first stage of collapse, and (b) the growth and subsequent rolling motion of the primary billows was less constrained by gravity when density was dominated by the rapidly-diffusing temperature field. The latter effect is illustrated by the large increase of potential energy between t = 1500 and 2000s (dotted curve on figure3a). That potential energy was released as three-dimensional kinetic energy between t = 2000 and 2300s. Therefore, despite the relatively low initial growth rate of the KP mode,


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

10

3

1 2 3
10
2

DNS

OBS

4

y
0.5

0

Reb

z
0.5

10

1

0

10

0

0

1000

2000

3000

4000

5000 0

0.5

1

(a) rS z
z
0.5

time [s]

PDF

Figure 5: Left frame: evolution of the buoyancy Reynolds number for simulations 1-4. Right frame: probability distribution function (PDF) for buoyancy Reynolds number from observations in the main thermocline.
0 1

0

x

2

(b) rT

The buoyancy Reynolds number provides a useful description of the range of scales in stratified turbulence. It is defined as the ratio of the squared Kolmogorov eddy turnover rate, / , to the squared buoyancy frequency N 2 : Reb =
VT 2 VT

N

,

(8)

Figure 6: Partial densities S (a) T (b) for run 1 at t = 2286s. where the subscripts indicate volume averages over the This flow state is intermediate between those shown in figures turbulent subvolume VT . The turbulent kinetic energy dis- 1e and f. Values colored range from -0.4 (red) to 0.4 (dark sipation rate is defined locally as blue). Values outside this range are transparent. = 2 sij sij , in which sij = 1 2 uj ui + xj xi (9)

(10)

is the strain rate. Primes indicate fluctuations about the horizontally-averaged velocity u xy . When Reb is large, turbulent eddies are too energetic to be affected by buoyancy. Run 1 reached a buoyancy Reynolds number slightly in excess of 40 (figure 5, thick curve). The corresponding flow state is illustrated in figure 6. The remaining three runs shown in figure 5 were restricted to lower Reb , in part because the pairing instability was suppressed. Nevertheless, these runs are expected to give an accurate indication of the influence of the density ratio on turbulent diffusion. Note the slight difference in the evolution of Reb between the low and high density ratios (dashed and dotted curves on figure 5). 9

Also shown in figure 5 is a histogram of Reb taken from observations in the thermocline off Northern California (Moum, 1996b). A set of 994 profiles extending from 200m to a maximum of 600m depth was binned to yield 144246 1m segments, from which the statistics of Reb were computed. No attempt was made to isolate overturns or other regions of elevated turbulence. Values generally ranged between 1 and 103 ; the median was 29. Thus, the buoyancy Reynolds numbers attained in the DNS runs reported here appear to be representative of weakly turbulent regions of the ocean thermocline.


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

4. The scalar gradient spectrum
Power spectra of the scalar gradient fields provide a sensitive test of numerical resolution, as well as insights into the physics of turbulent mixing. Here, spectra are computed for the high Reynolds number case (run 1) at a time when turbulence was at its most intense (i.e. Reb was a maximum), so that the demand placed on spatial resolution was high. Figure 6 shows the two scalar fields from run 1 at this time, which is intermediate between the times shown in figure 1f and g. The right-going and left-going ends of the collapsing core are just beginning to interact at the domain center. The T field (figure 6b) displays a structure similar to the salinity, but with markedly less small scale variability (cf. GMH's figure 2). Figure 7 shows spectra of the scalar gradients / z versus the vertical wavenumber k . Spectra were computed in the vertical direction to facilitate comparison with profiler measurements. The symbol is used to denote either temperature or salinity. Spectra were computed from 500 vertical profiles sampled randomly within the domain. For each profile, the turbulent region was selected, was first-differenced, Hanning-windowed, and Fourier transformed to obtain the power spectral density z . A correction was applied to recover variance lost by firstdifferencing. Each spectrum was normalized prior to averaging, using the isotropic variance dissipation rate,


-1 k b ( q/2) -1/2 z

10

-1

10

-2

Sz Tz Sz Tz
-3

OBS OBS DNS DNS

10

-3

10

10

-2

10

-1

10

0

10

1

=( 2 q) 1/2k /k b

Figure 7: Normalized vertical gradient spectra of temperature
(thick dashed curve) and salinity (thick solid curve) for run 1 at t=2286s. Shown for comparison are gradient spectra of salinity (triangles) and temperature (circles) from 350 ocean turbulence patches (Nash and Moum, 2002). The thin solid curve is the Kraichnan (1968) universal form for the viscous convective and viscous diffusive subranges. The value 7.3 was used for the constant q (Smyth, 1999).

=

0

z dk ,

(11)

and the Batchelor scale kb = ( / 2 )1/4 (Batchelor, 1959). Also shown on figure 7 are spectra computed from observational data (Nash and Moum, 2002) and the theoretical spectral form of Kraichnan (1968). The T spectrum extends further into the small scales than does the S spectrum because the former field is somewhat better resolved with respect to its Batchelor scale (the ratio of Batchelor scales for the two scalars is 7 = 2.65; the ratio of grid spacings is 2.0). The spectra of small scale gradients determined from these simulations agree very well with both the observations and the theory. This indicates that the model is reproducing the small scale physics accurately, and in particular that the spatial grid resolution is adequate. At larger scales, correspondence is not as close. Large-

scale gradients are strongly affected by the evolving fields associated with gravity waves (GMH) and with the KH instability. In contrast, the theory assumes that the flow is in statistical equilibrium, and the observations have considerably larger Reynolds numbers, and thus less influence of the forcing scales in the viscous-convective and viscousdiffusive subranges, relative to the DNS. The DNS salinity spectrum peaks at a higher value than either the temperature spectrum, the theory or the observations. Both the salinity spectrum and the observations are systematically higher than the Kraichnan spectrum in the viscousconvective (k +1 ) range.

5. Potential energy, scalar variances and turbulent diffusivity
Our objective is to compare the turbulent diffusivities of the thermal and saline density components T and S in various parameter regimes. Here, we describe two approaches to this comparison, focusing first on the evolution of the horizontally-averaged scalar profiles and later

10


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

on an alternative approach that isolates irreversible mixing processes. Additional insight into the physics of differential diffusion is gained through examination of scalar variances, whose dissipation rates are used to estimate turbulent diffusivities in observational studies.

0.08 0.07 0.06 0.05

a. Component potential energies
In the context of vertical mixing of a scalar (which may represent either T or S ), computation of turbulent diffusivity requires fitting the evolution of the scalar field to a one-dimensional diffusion model, e.g.: Ї Ї = K . t z z (12)

0.04 0.03 0.02 0.01 0 0 PT PS PbT PbS 1000 2000 3000 4000 5000 6000

t [s]

In (12), the diffusion model is expressed in terms of the horizontally-averaged profile xy . Ї There are a number of ways to invert (12) in order to obtain a single, characteristic value for diffusivity at any given time. Here, we begin by considering changes in the specific potential energy associated with the evolution of each density component: g (z - Lz /2) V . (13) P (t) = 0 u 2 In contrast to and N 2 (cf. section 3), P is a global property of the flow. Accordingly, we make no attempt to isolate the turbulent region, but instead compute the average over the entire computational domain. Besides providing a route to the computation of K , P is in itself a useful descriptor of the flow physics. Potential energy components associated with temperature and salinity for run 1 are shown by the thick curves on figure 8. (Thin curves on figure 8 represent background potential energies, to be defined below.) For each scalar, the potential energy rose to a maximum, then decreased rapidly as the primary KH billows paired and subsequently collapsed. The potential energy then oscillated a few times before settling down to an approximately steady state. The oscillations indicate reversible transfers between the potential and kinetic energy reservoirs, associated with interference between leftgoing and rightgoing internal waves generated by the collapsing KH billow (e.g. figure 1f,g, figure 6). During the initial growth and pairing of the KH billows, the total potential energy stored in the temperature and

Figure 8: Evolution of the scaled potential energy components
for run 1. Thick curves: component potential energy. Thin curves: background potential energy. Solid curves: temperature component. Dashed curves: salinity component.

salinity fields increased at nearly equal rates, indicating that the two scalars were advected together. As the billows collapsed (the phase of rapid decrease in total potential energy), the component potential energies diverged. After turbulence had decayed, the temperature field contained more potential energy than did the salinity field. This indicates that salinity restratified more completely than temperature or, equivalently, that temperature mixed more thoroughly. It can be shown that, if the mean density evolved according to (12) with K independent of z , then K would be proportional to the time derivative of the component potential energy, v iz .: K (t) = Lz 0 u2 d P . 2g dt (14)

As a definition of K , (14) has a serious shortcoming: the resulting diffusivity is negative during times when P is decreasing. Negative diffusivity implies "unmixing" of a mixed fluid, an apparent violation of the second law of thermodynamics. The real problem, of course, is that the diffusion equation (12) is a poor model for the evolution of the mean profiles, because that evolution reflects not only diffusion but also the effects of gravity waves and other reversible processes. The roll-up and subsequent

11


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

breaking of the KH billows is an example: breaking does zontal mean but to the background state. The significance not represent a reversal of the diffusion process; in fact, it of these distinctions will be assessed below. is a time of extraordinarily rapid diffusion, as we show in To invert (15) and thereby obtain a characteristic value the next subsection. for K describing only irreversible processes, we first define the contribution to the minimum, or "background" potential energy, Pbn (Winters et al., 1995), associated b. Background potential energies and turbulent with the density component , v iz :

diffusivities

We now describe an alternative definition for the turbulent diffusivity that filters out reversible effects. We begin by defining the reordered height coordinate z (x, y , z , t), which is the height a fluid parcel would end up at if the partial density distribution was allowed to relax adiabatically to a state where the corresponding component potential energy was a minimum. (Note that this reordering is done in three spatial dimensions, not in one dimension as in the calculation of the Thorpe scale, (e.g. Thorpe, 1977).) Changes in this state reflect the effects of irreversible mixing alone (Winters et al., 1995; Scinocca, 1995; Winters and D'Asaro, 1996). A diffusion model that describes the evolution of the minimum potential en ergy state, (z , t), contains only irreversible effects: = t z


Pb (t) =

g 0 u

2

(z - Lz /2)

V



,

(17)

K



. z

(15)

The definition of the turbulent diffusivity implicit in (15) has a number of appealing properties. K is positive definite; in fact its lower bound is the molecular viscosity, achieved when partial density distribution is statically stable and the fluid is motionless. The ratio of turbulent to molecular diffusivity on any isosurface of is equal to the square of the ratio of the area of that isosurface to its area in the stable, motionless state (Winters and D'Asaro, 1996). That ratio is also equal to a ratio of gradients very similar in form to the Cox number which appears via the standard Osborn-Cox formulation for stratified turbulence: | |2 z K = (16) )2 . ( / z Note, however, that the right-hand side of (16) differs from the usual Cox number in that the squared gradient is averaged not over coordinate planes but over isoscalar surfaces (since z is a function of only). Also, the vertical gradient in the denominator pertains not to the hori-

where the subscript on the angle brackets indicates a volume average taken over the background state (or, equivalently, over isoscalar surfaces instead of coordinate planes). To distinguish it from Pb , the potential energy component P defined earlier is referred to as the "total" potential energy due to the density component . The difference P - Pb is called the available potential energy, as it is available for conversion to kinetic energy. Note that P is the "total" potential energy only in the sense that it includes both the background and the available potential energies; it nevertheless refers only to the contribution of the density component . The potential energy contained in the complete density field (discussed in section 3) is given in terms of the component potential energies by P = PT + PS . The background potential energy components PbT and PbS (thin curves on figure 8) respond only to irreversible processes, and they therefore increased monotonically throughout run 1 (and all other runs). The thermal component increased more rapidly than the saline component right from the beginning of the primary growth phase. About one half of the eventual divergence of PbS and PbT occurred before the transition to turbulence was complete. For both scalars, the increase in background potential energy was steepest (i.e. irreversible mixing was most rapid) during the collapse of the billow between t = 1700s and t = 2100s, as indicated by the rapid loss of total potential energy. Throughout this early period of differential diffusion, the total potential energies stored in the temperature and salinity fields (thick curves on figure 8) remained nearly equal. This shows that the difference in the background potential energy increases was compensated in the available potential energies. Vertically displaced fluid parcels create background potential energy by mix-

12


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

10

-5

(a)

1

KT [m2/s]

1 2 3 4

10

-6

0.8

10 10

-7

(b) KS [m2/s]
10
-6

1 2 3 4

KS/KT
6000

-5

0.6

0.4

10

-7

0.2

10

-8

0

1000

2000

3000

4000

5000

0 0

1000

2000

3000

4000

5000

6000

t [s]

t [s]

Figure 9: Evolution of the instantaneous turbulent diffusivities of temperature (a) and salinity (b) for runs 1 (thick solid) 2 (dashed), 3 (thin solid) and 4 (dotted). Horizontal lines indicate the molecular diffusivities.

Figure 10: Ratio of instantaneous turbulent diffusivities for
runs 1 (thick solid), 2 (dashed), 3 (thin solid) and 4 (dotted). Horizontal lines indicate unity and the molecular diffusivity ratio 0.14.

ing with their surroundings, but at the same time give up available potential energy. Only after parcels lose their available potential energy via restratification does the difference in the irreversible mixing of heat and salt show up as a difference between the total potential energies. As turbulence decayed, the available potential energy stored in each scalar field due to waves and turbulent eddies dropped to zero, and hence the total and background potential energies for each scalar became equal. The temperature component of the background potential energy increased more rapidly than the salinity component throughout the run. The net amount of temperature mixing, as indicated by the net change in the associated background potential energy, was greater than that due to salinity, signalling differential diffusion. If the background density evolved according to (15) with K independent of z , then K would be given by K (t) = Lz 0 u2 d Pbn . 2g dt (18)

We adopt (18) as our definition of the instantaneous turbulent diffusivity. Figure 9 shows the instantaneous turbulent diffusivities for runs 1-4. Initially, the diffusivities for the different runs increased together, reflecting the very similar values of KH growth rate in the four cases. As the

billows reached large amplitude, however, the results diverged. In the three cases where pairing was suppressed the thermal diffusivity rose to about 25 times its molecular value near t = 1400s, then decreased. In run 1, thermal diffusivity continued to rise due mainly to additional mixing resulting from the pairing instability, eventually peaking at 40 times its molecular value near t = 2200s. Note that this time coincides with the time of maximum buoyancy Reynolds number (cf. figures 5 and 6). The saline diffusivity was generally smaller, though larger in proportion to its molecular value. The ratio of instantaneous diffusivities (figure 10) increased initially from its molecular value of 0.14 toward values near unity. This increase occurred mostly during a dramatic jump that coincided roughly with the appearance of the KP instability (figures 2, 3). In run 1, the increase was spread out, presumably owing to the influence of the pairing mode. The maximum value of KS /KT depended heavily on R , exceeding unity for the case R = 5. In all cases, KS /KT eventually decreased. (Simulations continued to very long times have confirmed that KS /KT eventually returns to its molecular value.) The high Reynolds number case (run 1) did not achieve the highest maximum KS /KT ; however, the ratio remained close to unity long after it had begun to decrease

13


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

in the other cases. We will see below that this difference tended to reduce the difference in the cumulative diffusion rates of heat and salt in the higher Reb case. To understand the differences in the decay of KS /KT in the four cases shown in figure 10, we must first explore the physics of scalar mixing in terms of temperature and salinity variances.

d/dt
2

-

2

[10-4 s-1]

0

0

-2

-2

(a) T
0 1000 2000 3000 4000 5000 6000 0

(b) S
1000 2000 3000 4000 5000 6000

t [s]

t [s]

c. Scalar variances
In these simulations, volume-averaged scalar variances evolve according to d = - , dt where =
2 V

Figure 11: Terms in the variance equation (19) for temperature
(a) and salinity (b).

(19) of the diffusivity ratio for irreversible mixing processes regardless of the validity of the Osborn-Cox theory. /2 , The temperature and salinity variances in run 1 evolved in very similar fashion (figure 11). Early in the run, d/dt increased rapidly due to strong production. The dissipation term became important gradually as gradients sharpened. A second peak in the production rate corresponded to the pairing instability. The two breaking events (rapid decreases in potential energy in figure 2) were characterized by strong negative production as the rolling of the KH vortex reduced scalar fluctuations about the horizontal mean. The dissipation rate reached a maximum during this time. The late evolution was dominated by dissipative decay, with only weak and fluctuating production rates. In figure 12 we show the scaled ratio of dissipation rates, d , along with the instantaneous diffusivity ratio and the buoyancy Reynolds number. Here, the dissipation rates S and T are computed using the full scalar fields (including mean profiles) to ensure that their ratios remain well-defined when turbulence is weak. Note first that dissipation and diffusivity ratios remained very nearly equal over most of each run, diverging only by 10% as the flow reached its most turbulent state. This is somewhat surprising, since figure 11 shows that the productiondissipation balance assumed in the Osborn-Cox formulas is satisfied only when averaged over the whole event; the instantaneous production and dissipation rates show no relationship whatsoever. Recall, however, that the relationship between dissipation rates and diffusivities (16) does not require the fields to be in equilibrium if the latter represents irreversible processes only, as it does here. The relationship remains imperfect because is averaged over

/ 2 ,

= -2 w
V

V

= 2 | |2

/

2

represent the variance, production rate and dissipation rate of the scalar . Primes indicate fluctuations about the horizontal mean. All quantities are normalized by 2 to facilitate comparison between temperature and salinity variance budgets. In the Osborn-Cox formulation for stratified turbulence (Osborn and Cox, 1972), turbulent diffusivities are proportional to scalar variance dissipation rates. This relationship requires that the scalar field be in statistical equilibrium, so that = , and that the scalar flux represented by obey a flux-gradient relationship. This formulation is used to estimate turbulent diffusivities from ocean microstructure measurements, and the ratio d = S R T
2

(20)

was used as a surrogate for KS /KT in the observational analyses of Nash and Moum (2002). A more general formulation by Winters and D'Asaro (1996) resulted in (16), which relates the irreversible scalar flux (and hence K ) to the dissipation rate averaged on isoscalar surfaces without the need for an equilibrium assumption. The WintersD'Asaro formulation is applicable to three-dimensional solutions but cannot be realized directly from field data, as the latter is generally one-dimensional. However, results given below suggest that d is actually a useful estimate

14


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

1

1

(a) d , KS/KT

(b)

the accompanying discussion): k , kb (21) where f is an unspecified nondimensional function. The superscript "0" on 0 , 0z (k ) and 0 indicates evalua tion at t = 0. Now suppose that, at t = 0, all motion is brought instantaneously to a halt, leaving the scalar field to diffuse with no turbulent straining. Scalar evolution is now governed by a simple, linear diffusion equation. Each Fourier mode decays exponentially, and the gradient spectrum therefore evolves according to 0z (k ) = f ( ); = (2q )1
/2

0

0

0 1

0

1

(c)

(d)

40

d , KS/KT

20

0 0

1000

2000

3000

4000

00 5000 0

1000

2000

3000

4000

0 5000

t [s]

t [s]

Reb

10

Reb

10

10

0 kb

q 2

1 /2

Figure 12: Comparison of the scaled dissipation ratio d (solid
curve) with the buoyancy Reynolds number (dashed curve) and the ratio of turbulent diffusivities (dotted curve) for runs 2 (a), 3 (b), 4 (c) and 1 (d).

z (k , t) = 0z (k )e-

k2 t

.

(22)

We may then calculate the evolution of the dissipation rate (t) using (11), which under the Batchelor scaling becomes (t) = where e = 1 q
0

0 2



f ( )e-e
0

2 t

d ,

(23)

coordinate planes rather than isoscalar surfaces as in (16), but that discrepancy is evidently important only during a brief phase when turbulence is strongest. Note also that the ratios remain close to their maximum values for a significant time after turbulence intensity, as measured by Reb , has begun to decrease. As was seen for the diffusivity ratio in figure 10, there is a marked difference among runs in the time taken for the ratios to drift away from their maximum values. In particular, the scalar fields in run 1 retained the characteristic that d KS /KT 1 long after turbulence had decayed. The fact that the scalar field retains the characteristics of turbulence for a time after turbulence has decayed is not surprising since both scalars diffuse less rapidly than does momentum. However, the origin of the differences between runs is less obvious. To understand the fact that the diffusivity ratios remain high for so long in run 1, consider the following thought experiment. Suppose that, at some time that we arbitrarily designate as t = 0, a turbulent flow with energy dissipation rate 0 carries a passive scalar whose gradient spectrum obeys the Batchelor(1959) scaling (cf. figure 7 and

1/2



(24)

is the effective compressive strain rate of turbulent eddies (e.g. Smyth, 1999) just prior to the arrest of motion at t = 0. This turbulent strain rate controls the rate of diffusion by controlling the Batchelor scale ( /e )1/2 : stronger strain generates fluctuations on smaller spatial scales, which then diffuse more rapidly when the strain is switched off. The important observation here is that the integral in (23) is independent of the molecular diffusivity. Therefore, if is independent of at t = 0, it will remain so as the scalar fluctuations diffuse. More generally, the ratio of the dissipation rates of two scalars having different molecular diffusivities will not change as the fluctuations diffuse. The foregoing argument rests on the assumption that the scalars in question are passive, i.e. that buoyancy effects are not important. In the present experiments, buoyancy effects are present and become increasingly dominant as turbulence decays. Buoyancy adds a new time scale, N -1 , to the problem, invalidating the Batchelor scaling and with it the above analysis. The mean shear

15


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

1 1 0.9 0.8 0.7 0.6 2 3 4

that characterizes the whole mixing event. For this reason, we define the cumulative diffusivity of in terms of the net change in the associated background potential energy, Pb . Because, the "end" of the event is chosen arbitrarily, we first let the cumulative diffusivity be a function of time: K
C

KSC/KTC

0.5 0.4 0.3 0.2 0.1 0 0

(t) =

Lz 0 u2 Pb (t) - Pb (0) - t . 2g t

(25)

1000

2000

3000

4000

5000

6000

7000

tf [s]

Figure 13: Evolution of the ratio of cumulative turbulent diffusivities for runs 1 and 2-4. Early times are not shown because the ratio KS C /KT C is undefined before the onset of turbulent mixing.

also adds a time scale, but mean shear is nearly proportional to N during the decay phase since the bulk Richardson number remains nearly constant (Smyth and Moum, 2000b). The two effects therefore become important at about the same time. The ability of buoyancy and shear to influence the dynamics depends on the ratio of the decay rate e to N , which is proportional to the square root of the buoyancy Reynolds number. Therefore, when Reb is small (as in runs 2-4), we expect that d will drift rapidly away from the value it had before turbulence began to decay. Conversely, d should remain close to its turbulent value for longer when the turbulent phase is characterized by larger Reb , as in run 1. The above argument pertains entirely to the dissipation rates. We know of no corresponding argument to explain the fact that the ratio of turbulent diffusivities remains high for longest when Reb is large, other than to note the evident fact that the two ratios were very similar during the decay phases of these simulations (figure 12).

The constant = 2g /Lz 0 u2 is the rate at which potential energy would increase if the fluid remained in the stable motionless state. This rate is determined entirely by the potential energy fluxes at the upper and lower boundaries (Winters et al., 1995; Winters and D'Asaro, 1996), and therefore remains steady as long as the mean densities on the upper and lower boundaries do not change appreciably, as is the case in the present simulations. By subtracting out this relatively small increase, we isolate potential energy changes due to fluid motions. We next define the ratio of the cumulative turbulent diffusivities: T PbS (t) - PbS (0) - S t KS C = . (26) KT C S PbT (t) - PbT (0) - T t This ratio as mixing turbulence provides a (figure 13) was undefined at early times, rose proceeded, then approached an asymptote as decayed. The asymptotic value of KS C /KT C useful metric for differential diffusion: d = lim
t

KS C (t) . KT C (t)

(27)

Note that, had we not subtracted t from the potential energies in (25) and (26), this asymptotic limit would not exist and the cumulative diffusivity ratio would not be well defined. (In contrast, we did not subtract when defining the instantaneous diffusivities in (18). Had we done so, the instantaneous diffusivities would have approached zero at early and late times, and their ratio d would then have been undefined.) That characteristic diffusivity ratio was about 0.82 for the high-Re case (figure 13, solid curve). The maximum ratio was significantly lower for the lower-Re cases, and varied by 14% over the range of R in runs 2-4 with the lower values corresponding to R < 1. The latter variation is in agreement with M04. In these experiments, the d. Cumulative diffusivities dependence on R was due mainly to differing time inFrom a parameterization perspective, we care less about tervals over which instantaneous KS /KT remained large the instantaneous diffusivity than about a net diffusivity (cf. figures 10, 12 and the accompanying discussion). 16


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

6. Comparison with previous work
0.8 0.75 0.7 0.65 0.6

In this section we survey results from previous laboratory, observational and DNS studies of differential diffusion. Figure 15 shows values of d and Reb from the present work along with representative results from previous studies.

d

a. Laboratory experiments
0.55 0.5 0.45 0 10
1 2

10

10

Reb

Figure 14: Ratio of cumulative turbulent diffusivities versus
buoyancy Reynolds number. The four larger runs 1-4 already discussed; smaller symbols runs included to give a more comprehensive determining d. Filled circles: R = 1. Open Asterisks: R = 5.0. symbols represent represent auxiliary view of the factors circles: R = 0.2.

The diffusivity ratio shows a close correlation with the maximum value of Reb , as illustrated in figure 14. This is in accord with GMH as well as with other studies as detailed below. Figure 14 includes results from all ten DNS runs listed in table 1. The runs covered about an order of magnitude of variation in Reb , and exhibited d values ranging from 0.51 to 0.82. Evident again is the tendency for d to increase with increasing density ratio. Runs 7, 8 and 9 had R = 0.2, 1.0 and 5.0, respectively. The buoyancy Reynolds number reached 7 in each case, so that the three cases line up vertically at that value on figure 14. The difference in d among these three cases was 20%. The difference between this result and the results of runs 2-4 described above indicates that the effect of R is b. Ocean observations most marked at low Reynolds number, consistent with the expectation that d should approach unity at high Reynolds Nash and Moum (2002) analyzed 350 turbulent patches number for all R . measured over the continental shelf off Oregon. Using a 17

The initial lab experiments of Turner (1968) have been re-analyzed by Nash and Moum (2002) in order to estimate the buoyancy Reynolds number. Equating the ratio of entrainment fluxes with the diffusivity ratio, Nash and Moum obtained the relation shown by the thick curve in figure 15. The diffusivity ratio increases with increasing Reb until it reaches a value near unity at Reb 102 . The thickness of the curve represents the uncertainty in the estimation of Reb from the original data. In the laboratory experiments of Jackson and Rehmann (2003), the work done on the fluid was measured in order to infer the kinetic energy dissipation rate and hence Reb . Beginning and ending profiles of temperature and salinity yielded the diffusivity ratio. The results, indicated by crosses on figure 15, fell into two broad groups based on buoyancy Reynolds number. In the upper group, Reb ranged from 500 to about 25000, and KS was generally larger than KT by a few percent. In the lower group, Reb was within a factor of three of 102 and Ks /KT ranged between 0.56 and 0.87. Hebert and Ruddick (2003) measured differential diffusion in internal gravity waves generated by a paddle in a uniformly stratified fluid. To avoid the effects of heat losses from the sidewalls, they measured differential diffusion of a pair of chemical dyes having different molecular diffusivities in place of heat and salt. The flows lay in a very different region of parameter space from the other results surveyed here (and are therefore not shown on figure 15), but the general trend of increasing KS /KT with increasing Reb was reproduced.


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

fast-response conductivity/temperature probe on a slowlyfalling profiler, they obtained the first in situ estimates of the salinity dissipation rate. Assuming a productiondissipation balance in the scalar variance budgets, they estimated the diffusivity ratio from the ratio of the thermal and saline dissipation rates and hence of d . Three unique aspects of these field observations are pertinent. First, in contrast with laboratory and numerical experiments, it is very difficult to follow a turbulent event through its cycle of growth and decay (in fact, this has yet to be accomplished in the field). Therefore, each observational estimate of d is based on the instantaneous characteristics of a turbulent event at some unknown stage in its evolution. This has been identified as a source of scatter in observational estimates of d , although the present results suggest that d may be relatively insensitive to the stage of the event at which it is measured (at least at high values of Reb after turbulence has been fully established). Second, each observation represents a single profile at a single horizontal location whose relationship to the horizontal extent of the turbulence is unknown. In contrast, laboratory and numerical estimates of d represent spatial averages. This is likely to be a major source of scatter in the observational results. Third, determination of KS and KT requires sufficiently strong signal to noise ratio in the raw data to resolve dissipation rates. This required that the Nash and Moum (2002) analyses be restricted to regions where values of Reb were higher than is common in the main thermocline (Moum, 1996a). For the bulk of the turbulent patches analyzed in Nash and Moum (2002), Reb was between O(102 ) and O(104 ), whereas measurements described in section 3 of the present study yield a median value of 29. The Nash and Moum (2002) analyses yielded a mean diffusivity ratio between 0.6 and 1.1. These limits, and the corresponding range of Reb , are represented on figure 15 by the lightly shaded ellipse. It was not possible to rigorously differentiate the mean value of d from unity, or to detect any trend with respect to Reb . This is due both to scatter and to the fact that the signal to noise requirement effectively restricted the analysis to a region of parameter space where d is close to unity and independent of Reb (according to laboratory and numerical results, e.g. figure

15). The determination of KS /KT in low Reb regimes is the object of another field study.

c. Direct numerical simulations
The DNS experiments of GMH were comparable to those described here, albeit with some significant differences. As in Merryfield et al. (1998), the initial condition consisted of uniform stratification and no mean shear. Turbulence was driven by the impulsive forcing of a finiteamplitude random velocity field (as opposed to the smallamplitude perturbations used here to catalyze the generation of turbulence via dynamic instability). As in the present case, turbulence intensity grew and decayed in time. In contrast to the present case, GMH's flows were spatially homogeneous in a triply periodic computational domain. (Both methods have advantages. GMH's approach offers efficient access to higher Reynolds numbers because the entire computational domain is occupied by turbulence, whereas the present approach is more realistic in the sense that turbulence is generated via a physically realizable flow instability known to be important in the ocean.) GMH quantified differential diffusion in various ways that did not include the ratio of turbulent diffusivities employed here. They did, however, compute the ratio of the time-integrated buoyancy fluxes due to heat and salt. Because the integrated buoyancy flux is equivalent to the potential energy gain,and taking account of GMH's buoyancy scaling, the ratio of fluxes should be equivalent to the ratio of diffusivities. For buoyancy Reynolds numbers ranging from 0.4 to approximately 103 , GMH's results give diffusivity ratios between 0.32 and 0.94, as shown by the shaded circles on figure 15. GMH found that differential diffusion becomes more pronounced at lower Reb (shaded circles on figure 15), as have we in the present study. For Reb less than about 100, both GMH and the present DNS study found levels of differential diffusion that were generally less pronounced than the comparable results from the laboratory experiments (shaded band and crosses). This could be due to the fact that the difference in molecular diffusivities was artificially reduced in both DNS studies. M04 has extended the computations of GMH to include variation of R and the limit of zero stratification. In the

18


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

1.1 1 0.9 0.8

d

0.7 0.6 0.5 0.4 0.3 0.2 -1 10
0 1 2 3 4

has been analyzed in terms of scalar gradient spectra, the buoyancy Reynolds number Reb , scalar variances, and total and background potential energies associated with each density component. We have calculated both instantaneous and net turbulent diffusivities for each scalar in order to test the hypothesis that the two scalars would mix at different rates. This anticipated difference is a consequence of the difference in molecular diffusivities, and conflicts with the predictions of high Reynolds number turbulence theory. The ratio d describing the relative degrees of mixing of the two scalars (specifically, the ratio of the cumulative turbulent diffusivity of salt to that of temperature) was tested for dependence on Reb and on the density ratio R . For all cases, d was less than unity, but that ratio increased toward unity with increasing Reb and also with increasing R . The dependence on Reb is largely due to differences in the duration over which instantaneous KS /KT is large. Even weakly turbulent flows attain KS /KT 1 at maximum Reb , but when the peak Reb is larger, KS /KT remains close to unity long after turbulence has subsided. This finding suggests that the role of nonstationarity in differential diffusion is more complex than was previously thought. We have proposed an explanation for the persistence of KS /KT in terms of the Batchelor (1959) scaling for scalar gradients. The results were compared with results from previous laboratory, observational and numerical studies. Considering the wide range of flow geometries, parameter values and experimental techniques, the results summarized in figure 15 present a remarkably consistent picture. The diffusivity ratio is near unity for Reb > 102 , and some studies have suggested that this ratio actually exceeds unity for high Reb (e.g. M04). For buoyancy Reynolds numbers below O(102 ), heat diffuses more rapidly than salt. For Reb < 10, a circumstance that is common in the thermocline (see figure 5), the difference is greater than a factor of two. The results of the present study are consistent with this picture, and we may therefore add KH billows to the list of turbulent flows exhibiting differential diffusion. As in GMH and M04, the present levels of differential diffusion represent an underestimate due to the artificially increased diffusivity of salt. These findings reinforce the impression that Reb

10

10

10

10

10

Re

b

Figure 15: Results from the present study (black bullets), along
with summary results from selected previous studies: Turner (1968) (thick curve), Nash and Moum (2002) (ellipse), Jackson and Rehmann (2003) (pluses), Gargett et al. (2003) (grey bullets).

latter limit, M04 has found diffusivity ratios in excess of unity, and has given a physical explanation for this result. As in the present study, Merryfield finds that d increases with increasing R . This work has also revealed a close correlation between differential diffusion and restratification, the latter being quantified in terms of Lagrangian particle displacements. Unfortunately, we cannot duplicate that calculation with the present data as it does not include particle trajectories. A similar calculation using potential energy evolution to quantify restratification revealed no consistent correlation.

7. Summary
We have described a sequence of direct computations of the growth and decay of turbulence driven by KH instability. The flows were stratified by a combination of two active scalars representing temperature and salinity. All parameter values were consistent with weak mixing events in the thermocline, except that the diffusivity of salt was increased to facilitate resolution of the smallest fluctuation scales with the available memory. Flow evolution

19


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

REFERENCES

102 is a useful estimate of the lower limit of the "high Reynolds number" regime of stratified turbulence, in which theoretical results valid in the limit of infinite Reynolds number remain accurate. Similar results have been found by Smyth (1999) for the Reb dependence of the Batchelor "constant" q , and by Itsweire et al. (1993) and Smyth and Moum (2000a) for various statistical relationships dependent on the assumption of local isotropy. If generally valid, this observation has important implications for the interpretation of microstructure measurements, as that science relies heavily on high Reynolds number theory. It is possible that the majority of ocean mixing is accomplished by high Reynolds number events. At any given time, however, large volumes of the ocean interior experience mixing weak enough (figure 5) that preferential diffusion of heat over salt should be anticipated (figure 15). The importance of such weak mixing events to the large scale circulation remains to be quantified. We are now working to extend the results of Nash and Moum (2002) via more extensive observations of microstructure in low Reynolds number mixing events, in combination with further analyses of the DNS experiments described here. We have so far confined our attention to the diffusively stable case in which both the thermal and saline components of density are stably stratified. In much of the ocean, one or the other of these components is unstably stratified, leading to the possibility of double diffusive instability and hence vastly more complex flow physics. The diffusively unstable case is now under investigation.

References
Altman, D. and A. Gargett, 1990: Differential property transport due to incomplete mixing in a stratified fluid. In E. List and G. Jirka, editors, Stratified Flows, 454­ 460, American Society of Civil Engineers. Batchelor, G. K., 1959: Small-scale variation of convected quantities like temperature in turbulent fluid. J. Fluid Mech., 5, 113­133. Collins, D. and S. Maslowe, 1988: Vortex pairing and resonant wave interactions in a stratified free shear layer. J. Fluid Mech., 191, 465­480. Corrsin, S., 1951: On the spectrum of isotropic temperature fluctuations in isotropic turbulence. J. Appl. Phys., 22, 469­473. Farmer, D. and L. Armi, 1999: Stratified flow over topography: the role of small-scale entrainment and mixing in flow establishment. Proc. Roy. Soc. A, 455, 3221­ 3258. Gargett, A., 2003: Differential diffusion: an oceanographic primer. Prog. Oceanogr., 56, 559­570. Gargett, A., W. Merryfield, and G. Holloway, 2003: Direct numerical simulation of differential scalar diffusion in three-dimensional stratified turbulence. J. Phys. Oceanogr., 33, 1758­1782. Gargett, A. E. and G. Holloway, 1992: Sensitivity of the GFDL ocean model to different diffusivities of heat and salt. J. Phys. Oceanogr., 22, 1158­1177.

Acknowledgements: This paper has benefited greatly from the comments of Ann Gargett and Bill Merryfield. Tim Pugh provided essential programming support. We Hebert, D. and B. Ruddick, 2003: Differential mixing by breaking internal waves. Geophys. Res. Lett., 30, 1758­ also thank Bill Merryfield for access to research results 1782. prior to publication and Rita Brown for editorial assistance. The work was supported by the National Science Foundation under grants OCE0095640 and OCE0136116. Howard, L., 1961: Note on a paper of John W. Miles. J. Fluid Mech., 10, 509­512. Computing resources were provided by the National Center for Atmospheric Research and by the National PartnerItsweire, E., J. Koseff, D. Briggs, and J. Ferziger, 1993: ship for Advanced Computational Infrastructure. Turbulence in stratified shear flows: Implications for interpreting shear-induced mixing in the ocean. J. Phys. Oceanogr., 23, 1508­1522. 20


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

REFERENCES

Jackson, P. and C. Rehmann, 2003: Laboratory measurements of differential diffusion in a diffusively stable, turbulent flow. J. Phys. Oceanogr., 33, 1592­1603.

over the continental shelf. J. Phys. Oceanogr., 33(10), 2093­2112.

Nash, J. and J. Moum, 2002: Microstructure estimates of Klaassen, G. and W. Peltier, 1985a: The effect of turbulent salinity flux and the dissipation spectrum of Prandtl number on the evolution and stability of Kelvinsalinity. J. Phys. Oceanogr., 32, 2312­2333. Helmholtz billows. Geophys. Astrophys. Fluid Dyn., Osborn, T. R. and C. S. Cox, 1972: Oceanic fine structure. 32, 23­60. Geophys. Fluid Dyn., 3, 321­345. Klaassen, G. and W. Peltier, 1985b: The onset of turbulence in finite-amplitude Kelvin-Helmholtz billows. J. Scinocca, J., 1995: The mixing of mass and momentum Fluid Mech., 155, 1­35. by Kelvin-Helmholtz billows. J. Atmos. Sci., 52, 2509­ 2530. Klaassen, G. and W. Peltier, 1991: The influence of stratification on secondary instability in free shear layers. J. Seim, H. and M. Gregg, 1994: Detailed observations of a Fluid Mech., 227, 71­106. naturally-occurring shear instability. J. Geophys. Res., 99(C5), 10,049­10,073. Kraichnan, R., 1968: Small-scale structure of a scalar field convected by turbulence. Phys. Fluids, 11, 945­ Smyth, W., 1999: Dissipation range geometry and scalar 953. spectra in sheared, stratified turbulence. J. Fluid Mech., 401, 209­242. Merryfield, W., 2004: Dependence of differential mixing Smyth, W. and J. Moum, 2000a: Anisotropy of turbuMerryfield, W., G. Holloway, and A. Gargett, 1998: Diflence in stably stratified mixing layers. Phys. Fluids, ferential vertical transport of heat and salt by weak 12, 1343­1362. stratified turbulence. Geophys. Res. Lett., 25, 2772­ Smyth, W. and J. Moum, 2000b: Length scales of tur2776. bulence in stably stratified mixing layers. Phys. Fluids, Merryfield, W., G. Holloway, and A. Gargett, 1999: A 12, 1327­1342. global ocean model with double diffusive mixing. J. Smyth, W., J. Moum, and D. Caldwell, 2001: The effiPhys. Oceanogr., 29, 1124­1142. ciency of mixing in turbulent patches: inferences from Metcalfe, R., S. Orszag, M. Brachet, S. Menon, and J. Ridirect simulations and microstructure observations. J. ley, 1987: Secondary instability of a temporally growPhys. Oceanogr., 31, 1969­1992. ing mixing layer. J. Fluid Mech., 184, 207­243. Thorpe, S., 1977: Turbulence and mixing in a Scottish Miles, J., 1961: On the stability of heterogeneous shear loch. Philos. Trans. Roy. Soc. London, A286, 125­181. flows. J. Fluid Mech., 10, 496­508. Moum, J., 1996a: Efficiency of mixing in the main ther- Turner, J., 1968: The influence of molecular diffusivity on turbulent entrainment across a density interface. J. mocline. J. Geophys. Res., 101(C5), 12,057­12,069. Fluid Mech., 33, 639­656. Moum, J., 1996b: Energy-containing scales of turbulence in the ocean thermocline. J. Geophys. Res., 101(C6), Winters, K. and E. A. D'Asaro, 1996: Diascalar flux and the rate of fluid mixing. J. Fluid Mech., 317, 179­193. 14095­14109. Moum, J., D. Farmer, W. Smyth, L. Armi, and S. Vagle, Winters, K., P. Lombard, J. Riley, and E. A. D'Asaro, 2003: Structure and generation of turbulence at inter1995: Available potential energy and mixing in densityfaces strained by solitary waves propagating shoreward stratified fluids. J. Fluid Mech., 289, 115­128. 21 on N and R . J. Phys. Oceanogr., (submitted).


S M Y T H et al., D I FF E R E N T I A L D I FF U S I O N I N K H B I L L OW S

REFERENCES

Winters, K., J. MacKinnon, and B. Mills, 2003: A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Oceanic Technol., 21(1), 69­94. Woods, J., 1968: Wave-induced shear instability in the summer thermocline. J. Fluid Mech., 32, 791­800.

22