## Abstract

It is argued that the mixing efficiency of naturally occurring stratified shear flows, *γ*=*Rf*/(1−*Rf*), where *Rf* is the flux Richardson number, is dependent on at least two governing parameters: the gradient Richardson number *Ri* and the buoyancy Reynolds number *Re*_{b}=*ε*/*vN*^{2}. It is found that, in the range approximately 0.03<*Ri*<0.4, which spans 10^{4}<*Re*_{b}<10^{6}, the mixing efficiency obtained via direct measurements of fluxes and property gradients in the stable atmospheric boundary layer and homogeneous/stationary balance equations of turbulent kinetic energy (TKE) is nominally similar to that evaluated using the scalar balance equations. Outside these *Ri* and *Re*_{b} ranges, the commonly used flux-estimation methodology based on homogeneity and stationarity of TKE equations breaks down (e.g. buoyancy effects are unimportant, energy flux divergence is significant or flow is non-stationary). In a wide range, 0.002<*Ri*<1, the mixing efficiency increases with *Ri*, but decreases with *Re*_{b}. When *Ri* is in the proximity of *Ri*_{cr}∼0.1–0.25, *γ* can be considered a constant *γ*≈0.16–0.2. The results shed light on the wide variability of *γ* noted in previous studies.

## 1. Introduction

The specification of eddy viscosity and diffusivities remains a central problem in developing predictive models for oceans, lakes and the atmosphere. For stably stratified flows, the vertical diffusivities can be defined in terms of vertical momentum (, ) and buoyancy ( or temperature () fluxes as
1.1
and
1.2
where *K*_{M}≡*K*_{U}=*K*_{V} and *K*_{T} are the eddy viscosity and diffusivity, respectively, and their ratio *Pr*_{tr}=*K*_{M}/*K*_{T} is the turbulent Prandtl number [1]. When temperature *T* is the major contributor to density *ρ*, as in dry atmosphere, freshwater lakes and some marine environs, the buoyancy fluctuations become *b*′=−*g*(*ρ*−*ρ*_{o})/*ρ*_{o}=−*gρ*′/*ρ*_{o}≈*αgT*′, and *K*_{T} can be interpreted as the mass diffusivity. Here, *g* is the gravity, *ρ*_{o} the reference density and *α* the thermal expansion coefficient. The instantaneous velocities (*u*,*v*,*w*) in the (*x*,*y*,*z*) directions are written in the usual forms *u*=*U*+*u*′, *v*=*V* +*v*′, *w*=*W*+*w*′, the temperature and the density where *U*,*V*,*W*, and represent the mean and the primes are the fluctuations. In the atmosphere, sonic anemometers and gradient measurements allow direct evaluation of *K*_{M} and *K*_{T} from (1.1) and (1.2) by turbulent fluxes (using covariance calculation) and the gradients of temperature and velocity components at a suitable vertical separation. Here, the linear overbar stands for averaging over specific time segments (assuming Taylor frozen turbulence hypothesis). In oceans and other water bodies, however, the measurement of fluxes still remains a technical challenge, and eddy diffusivities are obtained either indirectly via local microstructure measurements or directly by tracer dispersion observations over large space–time domains. In the former, the simplified equations for the budget of turbulent kinetic energy (TKE) and turbulent thermal variance ,
1.3
and
1.4
are used with the assumption of homogeneity (neglecting advection and diffusion of *q*^{2} and *Θ*^{2}), where is the shear production of TKE and is the buoyancy flux [1]. Here, mixing is assumed to be internal; that is, TKE is produced and approximately balanced locally *vis-à-vis* external mixing where the energy flux divergence and advective terms transport turbulence, which is generated elsewhere, to the mixing location [2]. Here,
1.5
are the dissipation rates of *q*^{2}and *Θ*^{2} in isotropic approximation, and *ν* and *D* are the molecular viscosity and diffusivity, respectively. As mentioned previously, *ε* and *χ* in natural waters are estimated using microstructure profiling measurements [3–8] as well as acoustic Doppler current profiler and acoustic Doppler velocimeter records [9–12].

Defining traditional flux Richardson number as [1]
1.6
or using a generalized definition of *Rf* for non-stationary and inhomogeneous turbulence [13],
1.6a
where the denominator accounts for all local and non-local sources of the TKE production, and introducing the mixing efficiency [14]
1.7
it is possible to write the stationary balance of TKE using equations (1.3) and (1.2) as
1.8
or in non-dimensional form,
1.9
where
1.10
is the mixing Reynolds number. Note that (1.6*a*) is equivalent to (1.6) if pure internal mixing is occurring, that is shear-produced turbulence is dissipating locally and consumed for buoyancy flux. For ocean mixing, Osborn [14] suggested *Rf*≤0.17, which gives an upper bound for *K*_{T} in (1.8) with *γ*≈0.2.

For the case where the temperature flux and gradients are measured directly *K*_{f}≡*K*_{T}, the mixing Reynolds number (1.10) becomes
1.10a
and the buoyancy Reynolds number is given by
1.11
where is the squared buoyancy frequency.

Conversely, following (1.2) and (1.4), the temperature diffusivity *K*_{T} for stationary turbulence can be estimated as
1.12
where can be directly measured, and *χ* can be estimated from temperature microstructure data. The corresponding mixing Reynolds number is (cf. equation (1.10*a*))
1.10b
and the mixing efficiency *γ*_{χ}≡*γ* can be estimated using (1.8) and (1.12) [4], provided that the assumptions of stationarity and homogeneity are satisfied simultaneously [15]. Another approach of estimating *γ*≡*γ*_{f} is the flux-based measurement of *K*_{T}≡*K*_{f} using (1.2) and (1.8) or *γ*≡*γ*_{o} using (1.7) with (1.6) or (1.6*a*), where the temperature flux is measured directly. For convenience, these are listed in table 1.

The objectives of this paper are (i) to compare the mixing efficiencies *γ*_{χ} and *γ*_{f} calculated by the two methods and (ii) to discuss the dependence of mixing efficiency on *Ri* and *Re*_{b}. The canonical flow type of interest here is the stratified shear flow that is common in oceans, lakes and the atmosphere. We will use data taken in the atmospheric nocturnal boundary layer at very high Reynolds numbers.

The present work is particularly helpful in delving into the variability of *γ* and therefore *Rf* observed in different studies that have employed a variety of techniques by analysing mixing mechanisms [15–17]. In the laboratory, Linden [18], Rohr *et al.* [19], Rohr & Van Atta [20] and Monti *et al.* [21], among others, have used non-stationary turbulence, whereas Strang & Fernando [22] and McEwan [23] used stationary or quasi-stationary flows, yielding somewhat different results. Numerical studies with evolving [17,24,25] and stationary flows [26] show disparate *Rf* and *γ* behaviour. In all, past results illustrate the dependence of *Rf* and *γ* on the mixing mechanism [27], the state of evolution of mixing [21,24,25,28], the nature of forcing (internal versus external [2,29]) and even the effect of topography [30]. Results of several experimental and theoretical studies that showed a clear dependence of *γ* on *Ri* were summarized in fig. 10 of Lozovatsky *et al.* [31], and the repercussions of variable mixing efficiency are discussed in Balmforth *et al.* [32]. The notion of Phillips [33] that *γ* grows with *Ri*, reaches a maximum at some Richardson number (0<*Ri*<∼1), and then decreases has been supported by laboratory experiments [18,34], measurements in a stably stratified atmospheric boundary layer [35] and used in the Goddard Institute for Space Studies (GISS) oceanographic numerical model of Canuto *et al.* [36]. Initial attempts to estimate *γ* from oceanic flux measurements and balance equations of *q*^{2} and *Θ*^{2} yielded a *γ*_{χ} that is significantly smaller than *γ*_{f} (a single depth towing measurements in a tidal front by Gargett & Moum [37]), but the mixing efficiency was confined to a narrow range (0.13–0.17) on the basis of vertical profiling through turbulent patches [38]. The authors noted a large uncertainty for *γ* estimates owing to technical difficulties of flux measurements in the ocean interior and a limited length of observational records.

It was our hope that the evaluations of *γ* from a comprehensive dataset that allows direct measurement (using *γ*_{o} and *γ*_{f}) and indirect evaluation (via *γ*_{χ}) of *γ* will help clarify the wide variability of mixing efficiencies reported in the literature. Furthermore, our dataset allows evaluation of *γ* as a function of governing variables of natural stratified sheared flows.

## 2. Governing parameters

Consider a stratified shear flow away from the boundaries, which is characterized by the buoyancy frequency *N* and mean shear , with their length scales of variation *L*_{N} and *L*_{Sh}, respectively. The turbulence in the flow is specified by a characteristic velocity *q*, length scale *l*_{tr} and time scale *τ*. The molecular parameters involved are *v* and *D*. Using *ε*≈*q*^{3}/*l*_{tr}, the governing dimensional variables become *N*, *Sh*, *ε*, *q*, *v*, *D*, *L*_{N}, *L*_{Sh} and *τ*, yielding the governing non-dimensional variables *Ri*=*N*^{2}/*Sh*^{2}, *Re*_{b}=*ε*/*vN*^{2}, *SN*=*Sh*×(*q*^{2}/*ε*), *τ*×*Sh*, *l*_{tr}/*L*_{Sh}, *l*_{tr}/*L*_{N} and *Pr*=*v*/*D*, where *l*_{tr}∼*q*^{3}/*ε* has been reintroduced. For slowly spatially varying flows or turbulent regions that are much smaller than *L*_{N} and *L*_{Sh}, the governing parameters become the gradient Richardson number *Ri*, buoyancy Reynolds number *Re*_{b} and Prandtl number *Pr*. The additional constraint is the evolution of turbulence in equilibrium with shear (whence the shear number *SN* becomes a constant [39]). At high Reynolds numbers, or at geophysical scales, it is possible to invoke the Reynolds number similarity, and assume that the flows are *Pr* independent.

If local shear production of turbulence is the dominant source of (internal) mixing with *l*_{tr}=(*L*_{Sh},*L*_{N}), the simplest stationary TKE balance *K*_{M}*Sh*^{2}≈*ε* gives the following relationship between *Re*_{b}, *Ri* and *R*_{m}:
2.1
where the turbulent Prandtl number *Pr*_{tr}=*K*_{M}/*K*_{T} itself can be a function of *Ri* whence *Ri* exceeds a critical Richardson number *Ri*_{cr}≈0.1–0.25 [40,41]. Conversely, *Pr*_{tr} tends to be close to unity for *Ri*≪*Ri*_{cr} [1,31,40,42]. Shih *et al.* [25] obtained an inverse dependence between *Re*_{b} and *Ri* similar to equation (2.1) (their eqn (5.1)), and also showed a relatively weak decrease in *Pr*_{tr} (from approx. 1 to 0.85) when *Re*_{b} increases from approximately 10 to 100. For *Re*_{b}>∼100, *Pr*_{tr} remained approximately constant. Although *Ri* and *Re*_{b} are formally independent governing parameters, but for shear-generated turbulence, they are related to each other in a specific parameter range. Nevertheless, it is important to note the external (the mean shear for *Ri*) and internal (the dissipation rate for *Re*_{b}) nature of these parameters in turbulent stratified flows.

## 3. Observations and data processing

The data (wind velocity and sonic temperature) used in this study were obtained during the vertical transport and mixing experiment (VTMX) in Salt Lake City, Utah conducted from 30 September to 7 October 2000 (see Monti *et al.* [43] and Princevac *et al.* [44] for a general description of the experiment). A meteorological mast at a gentle (approx. 4^{°}) mountain slope was equipped with two three-cup anemometers at heights *h*=2.0 (CA-1) and *h*=7.3 (CA-2) m above ground level (agl) and with two thermistor sensors at *h*=1.8 (T-1) and *h*=6.9 (T-2) m agl. The threshold speed of the anemometers was 0.5 m s^{−1}, with an accuracy of 1.5 per cent. The 5 min averaged data were collected to characterize the mean wind speed and air temperature. The high-frequency (10 Hz sampling rate) records of velocity components *u* (downslope), *v* (cross-slope), and *w* (upward) and temperature *T* were obtained at *h*=4.5 m agl using a sonic anemometer/thermometer (Applied Technologies, Inc. and Metek GmbH). The resolution and accuracy of data were 0.01 and 0.05 m s^{−1}, and 0.01^{°}C and 0.05^{°}*C*, respectively. A sketch showing positions of the instruments at the mast is given in figure 1*c*. Clear skies and light synoptic winds characterized the weather conditions from 30 September until 6 October. During the last night, wind increased up to 12 m s^{−1} owing to synoptic influence. Because of the very dry atmosphere (relative humidity during the experiment did not exceed 5–8%), we did not apply a moisture correction to the sonic temperature.

The 5 min averaged records of temperature, along-slope wind velocity and vertical components of temperature and momentum fluxes at *h*=4.5 m agl are shown in figure 1*a*,*b*. The night-time were predominantly negative as a result of the stable stratification near the ground (note 30 September 2000 is day 274 of 2006; this differs by 1 day from Monti *et al.* [43], where the leap year adjustment was not made). To obtain accurate estimates of *γ*, specific 20 min segments of data with almost constant fluxes and winds were selected (see a number of such segments for each night in figure 1*a*). Note that the homogeneity and stationarity of flow at specified segments were confirmed by comparing (1.6) and (1.6*a*). An example of , |*Sh*| and *Ri* records during a typical night (3–4 October) is given in figure 2, with seven segments chosen for the calculation of *γ*. The flow was usually weakly stable at the beginning of the night (*Ri*<0.1 in this case). At about midnight (figure 2), a local maximum of *Ri*≈2 was observed owing to a sharp drop of shear and a continuous increase in wherein the flow was dominated by quasi-periodic internal-wave oscillations [44]. Oscillations with longer periods of approximately 2–3 h were also observed [35] that could be attributed to global intermittency associated with periods of intense turbulence production and enhance stability. Continuous cooling of the surface and general reduction of mean shear towards the end of the night led to an increasing Richardson number. The natural meteorological variability during the observational period of seven nights allowed estimates of *γ* to be obtained over a wide range of *Ri* with high statistical confidence.

The mixing efficiency, as mentioned in table 1, was calculated directly using equations (1.6) or (1.6*a*) and (1.7) and equations (1.9) and (1.2) to obtain *γ*_{o} and *γ*_{f}, respectively, as well as indirectly to yield *γ*≡*γ*_{χ} (equations (1.9) and (1.12)). In order to estimate the mixing efficiency directly (i.e. *γ*_{o} in table 1), we first compare the flux Richardson numbers *Rf* and *Rf*_{II} defined differently for stationary (equation (1.6)) and non-stationary (equation (1.6*a*)) balances of TKE. Strong linear correlation between the two variables can be seen in figure 3 with *r*^{2}=0.97 for the best least-square linear regression *Rf*_{II}=1.1*Rf*. This result suggests that statistically, *Rf*_{II} and *Rf* are almost identical, supporting the assumption of homogeneity and stationarity for the dataset considered. Thus, we can confidently use equation (1.6) paired with equation (1.7) to calculate mixing efficiency *γ*_{o}≡*γ*_{f}.

The covariance and Reynolds stresses were computed for selected 20 min segments using the direct covariance method and integrating the corresponding co-spectra (see several examples of in figure 4*a*). A characteristic relative difference between the two temperature flux estimates was 11 per cent by amplitude, with 21 per cent standard deviation. The gradients of the mean temperature were taken as the finite differences between the 20 min averaged temperatures and located *Δh*=5.1 m apart. The TKE and temperature dissipation rates, *ε* and *χ*, were obtained from the inertial and inertial-convective subranges of the velocity and temperature spectra
3.1
at each *i* segment, and the 20 min averaged wind velocities were used to convert the frequency spectra of *w*′ and *T*′ to the corresponding wavenumber spectra *E*_{w}(*κ*) and *E*_{T}(*κ*), where *κ* is the horizontal wavenumber. To ensure high quality of the dissipation estimates, only those spectra that exhibited a near-perfect −5/3 subrange were used with the canonical values of spectral constants [45]; for the longitudinal flow component (which is *u* in our case), *c*_{ku}=0.52 and for transversal components (*w* or *v*), *c*_{kw}=*c*_{kv}=4/3*c*_{ku}=0.67 (an example is given in figure 4*b*). A concern was the significant variation of the Obukhov–Corrsin constant *c*_{T} reported in previous studies. Sreenivasan [46] suggested 0.3<*c*_{T}<0.5 with a tendency for lower *c*_{T} at higher Reynolds numbers; for that reason, *c*_{T}=0.3 was used in this study. The spectra for *u*′,*v*′ and *w*′ exhibited clear inertial subranges, which, as expected, were wider for *E*_{u}(*κ*) and *E*_{v}(*κ*) compared with *E*_{w}(*κ*). At many segments, the flow was not horizontally unidirectional, thus posing problems of selecting the longitudinal and transverse directions. Therefore, the estimates of *ε* were made using *E*_{w}(*κ*) with *c*_{kw}=0.67 because *w* is the unequivocal transversal component. Two methods were used for the *R*_{m} (equation (1.10)) calculation, based on *K*_{T}≡*K*_{f} using the temperature flux and gradient measurement in (1.2), i.e. *R*_{mf}, and on *K*_{T}≡*K*_{χ} obtained via (1.12) by estimating *χ* through *E*_{T}(*κ*), i.e. *R*_{mχ}.

## 4. Dependence of mixing Reynolds numbers *R*_{mχ} and *R*_{mf} on *Ri*

On the basis of arguments of §2, the diffusivities must be dependent on the gradient Richardson number *Ri* as well as the buoyancy Reynolds number *Re*_{b}. To explore the *R*_{m}(*Ri*) dependence, a combined plot of *R*_{mf} and *R*_{mχ} versus *Ri* is shown in figure 5. Both diffusivities therein are almost identically affected by *Ri* following approximately *R*_{m}∼*Ri*^{−3/2} [47]. This simple parametrization has been previously used in several numerical models [48]. The *Ri*^{−3/2} fit is, however, applicable for a relatively narrow range of *Ri*.

In figure 5, *R*_{mf} flattens at a value of *R*_{mf}≈2.2×10^{3} in the range of *Ri* between *Ri*_{cr}=0.25 and *Ri*=1, and may even show a slightly rising tendency at larger *Ri*, although the latter cannot be substantiated owing to availability of few experimental points. It is, however, a possibility that the *R*_{mf} samples at *Ri*>1 represent intermittent turbulent patches that were advected to the observational site rather than generated locally by weak vertical shear. In this case, for *Ri*>1, the local *Ri* does not represent the state of turbulence in the same way as that for continuous turbulence.

On the contrary, *R*_{mχ} continuously decreases beyond *Ri*_{cr}, but much slower than in the *Ri*_{0}−*Ri*_{cr} range, where *Ri*_{0}≈0.025, being about an order of magnitude smaller than *Ri*_{cr}. For *Ri*>1, *R*_{mχ} tends to a background value *R*_{mb}=600. In the original stably stratified layer, shear-induced mixing onsets at *Ri*=*Ri*_{cr}, which rapidly grows until the Richardson number decreases to *Ri*=*Ri*_{0}, and then it reaches a saturation level typical of non-stratified shear flows. At *Ri*<*Ri*_{0}, *R*_{mf} flattens, deviating from the −3/2 power law towards *R*_{mf}∼3×10^{5}, but *R*_{mχ} continuously increases at low *Ri*, reaching approximately 10^{7} at *Ri*≈10^{−3}.

The dependences of *R*_{mf} and *R*_{mχ} on *Ri* evident from figure 5 can be approximated by a scaling formula
4.1
with *s*=2. For *R*_{mχ}, the fit is shown by the dashed line (*b*), with *R*_{mn}=2×10^{7} and *R*_{0}=2.5×10^{−3}, whereas the heavy line (*c*) is drawn for *R*_{mf} with *R*_{mn}=3×10^{5} and *R*_{0}=2.5×10^{−2}. The formula (4.1) belongs to the family of negative *Ri* power-law parametrizations of eddy viscosity and diffusivity in stably stratified flows [49–55]. Different *s* values have been proposed, ranging between 1 and 2.5. The value of *Ri*_{0} in (4.1) is usually taken from 0.1 to 0.3; however, *Ri*_{0} as low as 0.02–0.05 has been used sometimes to satisfy experimental data [31]. The fitted *Ri*_{0}=2.5×10^{−3} is much smaller than the previously suggested values, and the corresponding *R*_{mn}=2×10^{7} is also unusually large. Perhaps, the presence of wall-induced turbulence in the present case may explain the anomalies; previous comparisons of (4.1) have been conducted with data taken from the thermocline or free shear flows.

This can be checked by applying the law-of-the-wall [45] test for vertical diffusivity *K*_{z}=*κu*_{*}*z*, where *u*_{*} and *κ*=0.4 are the friction velocity and von Karman constant, respectively. All the highest values of *R*_{mf} and *R*_{mχ} in figure 5 that correspond to the lowest *Ri*<10^{−2} represent seven segments of data obtained on 7 October under strong katabatic winds (the along-slope wind component varied between −10 and −12 m s^{−1} (figure 1). A characteristic estimate of *u*_{*} for these data is approximately 1 m s^{−1} (the quadratic law formula). Thus, at *h*=4.5 m agl, , and hence *R*_{m}≈1.5×10^{5}. This value matches well with the normalized diffusivity *R*_{mf} in figure 5 at the lower end of the *Ri* axis, but it is more than an order of magnitude smaller than the corresponding *R*_{mχ} for the same values of *Ri*. The test implies that the balance-based calculation of *R*_{mχ} produces unreliable estimates of the normalized diffusivity for very low Richardson numbers (*Ri*<0.03). Probably, we can make the same conclusion about the flux-based estimates *R*_{mf} for high *Ri*>1. Hence, the dependence of mixing Reynolds number *R*_{m} on *Ri* is best represented by equation (4.1) with *R*_{mn}=3×10^{5} and *R*_{0}=2.5×10^{−2}, which is shown by line (*c*) in figure 5.

## 5. Mixing efficiency

### (a) Calculations of *γ* and its dependence on *Re*_{b}

The estimates of *γ*_{χ} and *γ*_{f} can be obtained from the regression plots *R*_{mχ}(*Re*_{b}) and *R*_{mf}(*Re*_{b}) shown in figure 6*a*,*b*. According to (1.8), *γ* must be constant when *R*_{m} and *Re*_{b} are linearly dependent, which is satisfied for the regression *R*_{mχ}(*Re*_{b}) shown in figure 6*a* as
5.1
For *R*_{mf}(*Re*_{b}), however, the least-squared fit yields
5.2
(figure 6*b*), where *ξ*≈33 is a non-dimensional regression coefficient and the exponent *p*=1/2. The above leads to the *Re*_{b} dependence of
5.3
In addition, *γ*_{o}, which is equivalent to *γ*_{mf}, was evaluated (table 1) using independent measurements of *B* and *P*, and the results are approximated in figure 7 as
5.4a
where *c*_{o}≈1.5. Note that the approximations (5.1), (5.2) and (5.4*a*) are obtained with high statistical confidence, with coefficients of determination *r*^{2}=0.96, 0.92 and 0.93, respectively.

At first glance, *γ*_{χ}≈0.16 obtained in (5.1) for the nocturnal stably stratified atmospheric boundary layer is in good agreement with many previous estimates of *γ* for natural waters. In ocean mixing studies, *γ*=0.2 is frequently used [4,14,54,56–63]; however, higher (up to 0.4) and lower (approx. 0.1) values of *γ* have also been suggested (e.g. [64] and [38,65,66]). Limnologists [13,67,68] usually prefer *γ*=0.15–0.17, or in some cases, even smaller values, *γ*=0.04–0.06 [69], at very low stabilities. Recent analysis of DNS data [15,25] showed an approximately linear regression between *R*_{mχ} and *Re*_{b}, also supporting *γ*=0.17, but only in a relatively narrow range of *Re*_{b}=7–10^{2}. The authors identified this range as a stationary transition turbulent regime sandwiched between decaying and developing turbulence.

Further, the direct numerical simulation (DNS) data [25] produced *R*_{mf} that appeared to be proportional to 2*PrRe*^{1/2}_{b} for *Re*_{b}=10^{2}–10^{3}. This is identical to the functional dependence (5.3) shown in figure 6*b* based on our atmospheric dataset for much larger buoyancy Reynolds numbers, *Re*_{b}>10^{3}, and therefore the dimensionless constant *ξ* is different. Results of several laboratory experiments [70,71,72] compiled in Shih *et al.* [25], together with DNS data, suggest that *R*_{m}∼*Re*^{1/3}_{b} (in our notation) when *Re*_{b} increases from 10^{2} to 10^{5}. Note that the diffusivity in the laboratory experiments of Barry *et al.* [70] with grid-generated turbulence was calculated using the change of system's background potential energy before and after mixing events, and hence is an integral measure of different turbulent regimes. A slightly weaker than *Re*^{1/2}_{b} dependence can be suggested for the eight most energetic () samples shown in figure 6*b* that can be approximated by a power function *R*_{mf}∼*Re*^{0.4}_{b}. It is, however, reasonable to conclude that for the atmospheric nocturnal boundary layer *R*_{mf}∼*Re*^{1/2}_{b} for *Re*_{b}>≈10^{4}. The difference points to the sensitivity of *γ* to the calculation methods used as well as to the nature of the flow.

### (b) Interdependence between *R*_{mf} and *R*_{mχ}

In order to examine the disparity between the constant (equation (5.1)) and *Re*_{b}-dependent (equation (5.3)) mixing efficiencies evaluated using different methods, we analysed the relationship between mixing Reynolds numbers *R*_{mχ} and *R*_{mf} (figure 8). It appeared that 68 per cent of data in the plot occupy a relatively narrow range of intermediate values of *R*_{mχ} and *R*_{mf} (within the rectangle), wherein the linear regression *R*_{mf}=*c*_{R}*R*_{mχ} is valid with a relatively high coefficient of determination *r*^{2}=0.64 and a regression coefficient *c*_{R}=0.65, having 95 per cent confidence bounds from 0.53 to 0.79. The deviation of *c*_{R} from the ‘perfect agreement’ case of *c*_{R}=1 can be attributed to the uncertainties of evaluating *ε* and *χ* using Kolmogorov and Obukhov–Corrsin spectra, especially those associated with spectral constants *c*_{kw} and *c*_{T} in equation (3.1) and flux measurements.

The difference between *R*_{mf} versus *R*_{mχ} is substantial at high and low ends of the *R*_{m} diagram, where approximately 10^{6}<*R*_{mχ}<2×10^{7} and *R*_{mf}<(2–3)×10^{5}. This is equivalent to high and low values of *Re*_{b} (figure 6*a*,*b*). Note, however, that large and small values of turbulent variables are usually subjected to highest uncertainties.

The striking loss of the parity between *R*_{mf} and *R*_{mχ} (or *K*_{f} and *K*_{χ}) at small and large *Re*_{b} invites explanation, which has a bearing on the estimation of fluxes in natural waters. Because *Re*_{b}∼(*L*_{O}/*L*_{K})^{4/3}, where *L*_{O}=(*ε*/*N*^{3})^{1/2} and *L*_{K}=(*ε*/*ν*^{3})^{1/4} are the buoyancy (or Ozmidov) and Kolmogorov scales, respectively, *Re*_{b}≫1 implies very weakly stratified turbulence, where the TKE production essentially balances the dissipation, and fluxes are determined by r.m.s. velocity and temperature fluctuations (, *c*_{θ} being a correlation coefficient [73]). This flux saturation, a reflection of a weak gradient or entraining fluid from non-turbulent regions, also implies non-stationarity and hence unsuitability of indirect methods of flux evaluation. On the other hand, small *Re*_{b} implies lack of the inertial subrange, which is essential for indirect flux estimation. A stationary balance between production of TKE, buoyancy flux and the rate of dissipation can be established only in a specific (figure 8) intermediate range of *Re*_{b}.

In other words, in fully developed turbulence at high *Re*_{b}=10^{7}–10^{8}, a constant rate of mixing sustains until the density/temperature gradient almost completely erodes to a level that cannot uphold continuous growth of buoyancy (temperature) flux, whence the normalized diffusivity *R*_{mf} tends to saturation (compare the highest *R*_{mf} values in figures 5 and 6*b*). Less energetic turbulence (*Re*_{b}<10^{4}) confounded by stratification produces weak mixing, which is characterized by the smallest *R*_{mf}≈(2–3)×10^{3}. The results indicate that highly energetic (high *Re*_{b}) as well as underdeveloped (low *Re*_{b}) turbulence does not support stationary, non-diffusive and non-advective balance of *Θ*^{2} in stratified flows. These findings are consistent with the interpretation of DNS data by Shih *et al.* [25]. The difference is the actual range of *Re*_{b} (a narrow one in both studies) where stationary turbulence prevails. The results can also be cast in terms of *Ri*, which is discussed below.

### (c) Dependence of *γ* on *Ri*

The parametrization of mixing Reynolds number *R*_{mf} as a function of *Ri* given in §4 is a good representation for Richardson numbers below approximately 1. We have directly evaluated *γ*_{f}≡*γ*_{o} (equation (1.6*a*) paired with equation (1.7)) and plotted it in figure 9 as a function of *Ri*. The growing trend of *γ*_{o}(*Ri*) up to *Ri*≈1 seen here has been reported in several laboratory experiments, numerical and field studies [20,32,69,74–76]. Lozovatsky *et al.* [31] proposed an approximation,
5.5
to the GISS modelling results of Canuto *et al.* [36], in which the eddy coefficients are parametrized using specific damping functions. If we discard in figure 9 the two largest samples of *γ*_{o} corresponding to *Ri*≫1 as outliers (per discussion on two largest *R*_{mf} samples shown in figure 5 for *Ri*>4), then the data trend broadly mimics equation (5.4), although the scatter is high. A better fit to this specific dataset is
5.6
which nicely captures the main trend shown in figure 9 by a bold line. The squared box in figure 9 encompasses the same range of *Ri* as that in figure 5 (0.03<*Ri*<0.4) and the corresponding range of *R*_{m} boxed in figures 6 and 8. The box-averaged mixing efficiency appeared to be equal 〈*γ*_{o}〉=0.2, with the box-median value [*γ*_{o}]=0.165. This is close to *γ*=0.16 obtained using equations (1.8) and (1.12) in the range 3.7–5.7 (figure 6*a*). The box-averaged Richardson number 〈*Ri*〉=0.1 is smaller than the conventional critical value *Ri*_{cr}=0.25, but still is in the range where shear-induced turbulence and internal mixing are dominant.

The above result can explain why a majority of microstructures in stratified oceans and lakes, where mixing is generated mainly by shear instabilities with *Ri* slightly below critical, leads to *γ*=0.16–0.2 when evaluated using indirect methods (table 1). Note that atmospheric and laboratory measurements have spanned a wide range of *Ri* and *Re*_{b}, and hence exhibited significant variations of *γ*. A monotonic increase in *γ*_{o} with *Ri* can be seen from about 0.01 at *Ri*≈(2–3)×10^{−3} to 0.4–0.5 when *Ri* is approaching unity. However, in a limited range around *Ri* approximately 0.1, the mixing efficiency can be considered approximately constant close to 0.2. Note that in a weakly stratified upper oceanic layer, the median of *Ri* can be as low as 0.1 [31]. The cumulative distribution of the Richardson number is often well approximated by a lognormal probability law, showing that the probability of *Ri*<0.25 can be above 60 per cent and for *Ri*<1, can approach 80 per cent.

Considering that two variables *Ri* and *Re*_{b} are involved in determining fluxes (§2), we plot contours of *γ*_{o}(*Ri*,*Re*_{b}) in figure 10, which show that *γ*_{o} in the boxed area (*Ri* and *Re*_{b} range associated with shear-generated turbulence) is between 0.1 and approximately 0.3. On the basis of our data, the upper cut-off *Re*^{up}_{b} was identified as approximately 5×10^{5} (*Ri*<0.4), beyond which buoyancy effects are insignificant, whereas the lower limit was *Re*^{lw}_{b}≈10^{4} (*Ri*>0.03), signifying the absence of a clear inertial subrange due to suppression of local production of TKE. The variability of *γ*_{o} outside this *Ri*−*Re*_{b} range is possibly due to other types of TKE sources unrelated to mean shear and/or to non-stationary turbulence (developing at *Ri*<0.03, *Re*_{b}>∼10^{6} or decaying at *Ri*>0.3, *Re*_{b}<10^{4}).

## 6. Discussion and conclusions

Vertical heat, mass and momentum fluxes have not been directly measured in oceans and lakes until recently, and thus the eddy viscosity and diffusivities that are central to predictive modelling are inferred indirectly. One method relies on the assumption that the mixing efficiency *γ* is a constant in the expression *K*_{T}*N*^{2}=*γε*. Alternatively, a simplified scalar balance equation is used. The possible variability of *γ* was investigated in this paper based on sonic anemometer data obtained from the stable atmospheric boundary layer. The diffusivity *K*_{T}≡*K*_{f} and the corresponding mixing Reynolds number *R*_{mf} (equation (1.10*a*) and table 1) were evaluated using temperature flux and gradient measurements and *ε* via the kinetic energy spectra; *γ* so obtained was designated as *γ*=*γ*_{f}. The results were compared with *γ*=*γ*_{x}, which was evaluated using *K*_{T}≡*K*_{χ} and the corresponding *R*_{mχ} (equation (10*b*)) based on the commonly used (in physical oceanography) scalar dissipation technique (equation (1.12)). Dimensional arguments suggest that, at high buoyancy Reynolds numbers *Re*_{b},*γ* is a function of *Ri* and *Re*_{b}, but the dependence on molecular parameters such as *Pr* is negligibly weak [25].

It was found that *γ*≈0.16 is nominally a constant in the range 10^{4}<*Re*_{b}<10^{6}, and for our data, this corresponds to approximately 0.03<*Ri*<0.4. Both normalized diffusivities *R*_{mf} and *R*_{mχ} coincide in this regime (figures 5 and 8), leading to *R*_{m}∼*Ri*^{−3/2}, indicating an approximate equivalence between *γ*_{f} and *γ*_{x}. Outside the above parameter ranges, *R*_{m} can be parametrized as a function of *Ri* using equation (4.1). For very low and high *Re*_{b} and *Ri*, the basic assumptions underlying *K*_{T}*N*^{2}=*γε* were untenable, and should be used with circumspection. Furthermore, the accuracy of measurements, calibration of sensors as well as the usage of Obukhov–Corrsin spectral forms all affect the numerical value of *γ*.

The mixing efficiency was found to be a growing function of *Ri* (equation (5.5)) in a wide range, 0.002<*Ri*<1. It was a decreasing function of *Re*_{b}, according to *γ*∼*Re*^{−1/2}_{b} (equation (5.3)), in the range 3×10^{4}<*Re*_{b}<3×10^{7}, which is in agreement with the DNS results of Shih *et al.* [25] for 10^{2}<*Re*_{b}<10^{3}, but at odds with Gargett's [77] suggestion of *γ*∼*Re*^{−1}_{b} for approximately the same range of *Re*_{b} (based on laboratory data of Itsweire *et al*. [78]). If turbulence is anisotropic at smaller Reynolds numbers, as argued by Gargett [77], then it impacts the dissipation estimates and hence mixing efficiency calculations. Detailed studies on the dependence of *γ* with *Ri* and *Re*_{b} were not possible, given the difficulty of obtaining *γ* versus *Ri* and *Re*_{b} in nature when one or the other parameter is constant (equation (2.1)). The mixing efficiency may vary from approximately 0.01 at *Ri*≈(2–3)×10^{−3} to approximately 0.5–0.6 at *Ri*∼1. Phillips [33] suggested an increase of *Ri* up to a critical value *Ri*_{0}, following a *γ* decrease at higher *Ri*. This trend was implemented by Canuto *et al.* [36] in modelling, and was experimentally observed by Strang & Fernando [22], Guyez *et al.* [79] and others. A decrease in *γ* at high *Ri* was evident (figure 9 and equation (5.5)), but could not be confirmed using the present dataset, as in the case of stably stratified Arctic boundary-layer data presented by Grachev *et al.* [80,81]. In our case, the number of high *Ri* data points are only a handful, thus precluding any inferences.

At Richardson numbers close to or below a critical value *Ri*_{cr}∼0.1–0.25 (viz. approx. 0.03<*Ri*<0.4), direct measurements of *γ*=*γ*_{o} via fluxes and gradients could be approximately treated as a constant with a characteristic value between 0.16 and 0.2. This *Ri* range is most pertinent to shear-generated stratified turbulent layers of oceans and lakes, which may explain why *γ* is often observed to be close to 0.2 and treated as such [82]. In this range, the microstructure-based estimates *γ*_{x} and *γ*_{f} are consistent with *γ*_{o}, showing *γ*≈0.16. Comparison of different observations, nonetheless, is stymied by the dependence of *Ri* on measurement resolution, in particular, the separation with which the gradients are estimated. DeSilva *et al.* [83] showed that when this separation is greater than the buoyancy scale, the measurement may not be representative of local *Ri*, and our data were on the verge of this limit.

Our results help shed light on the differences between mixing efficiencies often encountered in oceanographic and limnological studies. Both constant and variable values of mixing efficiency are used, which forms the basis of closure in numerical models [84]. The present results show that a constant mixing efficiency can be used only in a limited range of governing parameters (*Ri*, *Re*_{b}), and many oceanographic measurements appear to be in this range [85]; therein turbulence is internally generated and approximately satisfies conditions of stationarity and homogeneity.

## Acknowledgements

We are thankful to the participants of the VTMX experiment who conducted the field measurements and to Andrea Dato (a former student of University of Rome ‘La Sapienza’, Italy) who actively participated in the data processing during his visit to the USA. This work was supported by the National Science Foundation (CMG) and ONR grants nos N00014-10-1-0738 and N00014-11-1-0709 (MATERHORN Programme).

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales I’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.