Abstract

A large set of laboratory, direct numerical simulation (DNS), and large eddy simulation (LES) data indicates that in stably stratified flows turbulent mixing exists up to Ri ∼ O(100), meaning that there is practically no Ri(cr). On the other hand, traditional local second-order closure (SOC) models entail a critical Ri(cr) ∼ O(1) above which turbulence ceases to exist and are therefore unable to explain the above data. The authors suggest how to modify the recent SOC model of Cheng et al. to reproduce the above data for arbitrary Ri.

1. Introduction

Second-order closure (SOC) turbulence models have been widely employed to describe geophysical flows (e.g., Launder et al. 1975; Lumley 1978; Zeman and Lumley 1979; André et al. 1978; Mellor and Yamada 1974, hereafter MY74; Mellor and Yamada 1982, hereafter MY82; Canuto 1992; Shih and Shabbir 1992; Cheng et al. 2002, hereafter C02).

SOC models are derived from the basic dynamic equations for the velocity, temperature, moisture, salinity, etc. Flows in geophysical settings come in two “flavors”: unstably and stably stratified. In the first case, a local approximation of the SOC models [level 2 in Mellor–Yamada (MY) terminology] is not physically justifiable since eddy sizes are commensurable with that of the “container,” thus making large eddies the rule rather than the exception; under those conditions, the computation of heat fluxes with a local model is hard to justify. To account for nonlocality is not an easy task, but some recent results have appeared that combine turbulence and plume models (Canuto et al. 2007).

Stably stratified flows are characterized by smaller eddies, and thus a local approximation may be more justified. The difficulty resides somewhere else. In such flows, one has two competing factors: a shear S that acts as a source of mixing and a stable temperature gradient N2 that acts as a sink. The two combine into a single parameter, the Richardson number Ri, defined as

 
formula

Here, U, V, and Θ are the mean velocity components and mean potential temperature, respectively; z is the vertical height; g is the gravitational acceleration; α is the coefficient of thermal expansion; S is the shear; and N is the Brunt–Väisälä frequency.

The physical question to be answered is the following: Is there a critical value of Ri at which the effects of stable stratification overpower the mixing action of shear so as to lead to zero mixing? In principle, the SOC models should be able to yield an answer. Equilibrium SOC models exhibit a finite critical Richardson number Ri(cr) beyond which the SOC equations have no real solutions. For example, the MY82 model yields Ri(cr) = 0.19, while the most recent model by C02 increases the value to Ri(cr) = O(1), thus overcoming the difficulty first raised by Martin (1985) that the MY model yielded too shallow ocean mixed layers. A stability analysis including nonlinearities (Abarbanel et al. 1984) also arrived at Ri(cr) = O(1).

However, even an Ri(cr) ≈ 1 is not satisfactory. In fact, several data show that mixing persists for Ri > 1. These data include meteorological observations (Kondo et al. 1978; Bertin et al. 1997; Mahrt and Vickers 2005; Uttal et al. 2002; Poulos et al. 2002; Banta et al. 2002), laboratory experiments (Strang and Fernando 2001; Rehmann and Koseff 2004; Ohya 2001), large eddy simulations (LESs; Zilitinkevich et al. 2007, 2008; Zilitinkevich and Esau 2007), direct numerical simulations (DNSs; Stretch et al. 2001), oceanic measurements (Mack and Schoeberlein 2004), and theoretical modeling (Sukoriansky et al. 2005).

Zilitinkevich et al. (2007) were the first to try to construct a model, which they called an “energy and flux-budget turbulence closure model,” that has no critical Richardson number. They then proceeded to show how the model can explain a large variety of data. Their model differs significantly from the “mainstream” closure models used in traditional SOCs. In particular, their parameterization of the pressure correlations differs from the “standard” ones developed by numerous authors and supported by theoretical and experimental evidences (Rotta 1951; Launder et al. 1975; Lumley 1978; Zeman and Lumley 1979; André et al. 1978; MY82; Canuto 1992; C02). In order not to spoil the abundant low-Ri data [e.g., the behavior of several variables in the planetary boundary layer (PBL)], the model required the introduction of additional parameterizations.

Given this situation, it was our goal to answer the following questions: Can the traditional SOC models encompass an arbitrary large Ri(cr), and what changes are needed to do so without spoiling their performance for small and medium Richardson numbers?

2. New Ri dependence

SOC models have been widely used in the past since they constitute a manageable set of equations that are able to incorporate different physical processes (buoyancy, shear, rotation, etc.) following relatively well structured rules. One major weak point of these models is the modeling of the pressure correlations (third-order moments). A brief historical discussion of the closures proposed for such functions is given in C02.

To quantify the problem, let us consider the dynamic equations for the second-order moments representing the Reynolds stresses, turbulent kinetic energy, heat flux, and the temperature variance. The closure-independent equations, including both potential temperature Θ and mean velocity Ui fields, are as follows (see C02):

  • traceless Reynolds stresses bij = − ⅔Kδij: 
    formula
  • kinetic energy K = ½: 
    formula
  • heat flux hi = : 
    formula
  • temperature variance : 
    formula

The explicit forms of the diffusion terms (D) can be found in C02 [their Eqs. (8b), (3b), and (4b)], but since it is not pertinent to the present discussion, we do not need to write their form explicitly. Consider instead the two pressure correlation terms in Eqs. (2) and (4), which are defined as follows [see C02, their Eqs. (2b) and (3b)]:

 
formula

As discussed in some detail in C02, the pressure (p) correlations in (6a) have been modeled as follows (considering only the slow terms):

 
formula

C02 and many other previous SOC models further assume that

 
formula

While the closures (6b) and (6c) have been widely used for more than 30 yr [see C02, their Eqs. (6) and (7)], they may not necessarily apply to the stably stratified flows that we are discussing here. The best way to assess their validity would be to compare them against LES data. The authors have searched for such data, and, while more than one group is carrying out such computations, no results are available as yet in the literature. Lacking that information, we suggest the following arguments.

First, many studies of stably stratified flows, beginning with the pioneering theoretical work of Weinstock (1978), pointed out that one must include a buoyancy damping factor to the time scale τ that enters the buoyancy spectrum B(k) = −τ(k)N2E(k). His derivation, which is presented in a simplified form in appendix A of Cheng and Canuto (1994), shows that in wavenumber space the damping factor has the form

 
formula

The physical reason behind (7a) is that in stably stratified flows, eddies work against gravity and lose kinetic energy, which is converted to potential energy. The presence of the factor in the denominator is analogous to the “energy absorbed from the particle motion” process, a damping factor, common in similar situations in plasma physics (Weinstock 1978). Canuto et al. (1994) found that (7a), once translated from k space in the form

 
formula

significantly improved the PBL simulation results. Canuto et al. (2001) and Cheng et al. (2005) also used (7b) in the PBL simulations with higher-order turbulence models. Later on in section 5, we show that when both shear and buoyancy are included, (7b) can be generalized to

 
formula

A second argument can also be formulated as follows. Let us substitute the second of (6b) in (4). For the case i = 3, one obtains for h3 = h

 
formula

which we shall employ shortly. In (8a) the last term is the dissipation term of h, which is usually neglected or modeled together with the pressure correlation term in front of it; τ is a relaxation time scale, while τJ is a dissipation time scale. On the other hand, consider the two-point closure model developed and tested against a large variety of turbulent flows (Canuto and Dubovikov 1996a, b; Canuto and Dubovikov 1997). In that work, the dynamic equations for the spectra E(k), Ez(k), J(k), and Eθ(k) were derived, and it was shown how integrating those over all wavenumbers yields the one-point dynamic equations given in Eqs. (2)(5) above (Canuto and Dubovikov 1998). The main advantage of using a two-point closure is that the pressure terms are no longer present since they are eliminated from the very beginning using the incompressibility relation, which makes the pressure field part of the nonlinear terms. The explicit form of the four dynamic equations is given in Canuto and Dubovikov (1998) and need not be repeated here. The derivation of the one-point equations is also discussed. For our purposes, we only need to consider the equation for the heat flux J(k), which reads

 
formula

where P33(k) = 1 − k2z/k2 is one component of the projection operator that appears when the pressure term is eliminated. The functions νt(k), χt(k) are the turbulent viscosity and diffusivity, the ratio of which can be taken to be nearly independent of k and identified with the Prandtl number σt, where ν and χ are molecular viscosity and diffusivity. Integrating (8b) over all wavenumbers to recover (8a), we note that with

 
formula

we derive the following expressions for τ and τJ:

 
formula

the ratio of which is given by

 
formula

In the last ratio in (8e), the numerator and the denominator are affected by Ri in similar fashion, and the variations are largely canceled out, so the ratio has only a relatively weak Ri dependence, while the turbulent Prandtl number σt is well known to increase nearly linearly with Ri (see Fig. 7). For example, Schumann (1991) proposed

 
formula

We proceed to model the dissipation time scale τJ as proportional to the turbulent kinetic energy dissipation time scale τ :

 
formula

Using (9a, b) in (8e), we obtain

 
formula

which is of the form (7c).

Fig. 7.

Turbulent Prandtl number σt [Eq. (17c)] from the present model (solid line) vs Richardson number, Ri, compared to meteorological observations (Kondo et al. 1978, slanting black triangles; Bertin et al. 1997, snowflakes), laboratory experiments (Strang and Fernando 2001, black circles; Rehmann and Koseff 2004, slanting crosses; Ohya 2001, diamonds), LES (Zilitinkevich et al. 2007, 2008, triangles), DNS (Stretch et al. 2001, five-pointed stars). The previous SOC model (C02) result is plotted as a dashed line.

Fig. 7.

Turbulent Prandtl number σt [Eq. (17c)] from the present model (solid line) vs Richardson number, Ri, compared to meteorological observations (Kondo et al. 1978, slanting black triangles; Bertin et al. 1997, snowflakes), laboratory experiments (Strang and Fernando 2001, black circles; Rehmann and Koseff 2004, slanting crosses; Ohya 2001, diamonds), LES (Zilitinkevich et al. 2007, 2008, triangles), DNS (Stretch et al. 2001, five-pointed stars). The previous SOC model (C02) result is plotted as a dashed line.

We have concentrated on τ rather than on the other time scale τ for a variety of reasons. A first, quite practical, reason was that a Ri dependence of the latter did not seem to be able to fit any of the available data. A second reason is that it has been known for many years that stable stratification reduces the heat flux more than the momentum flux (Gerz et al. 1989). A third reason is that a dependence on Ri means an interplay between temperature (acting like a sink) and mean velocity fields (acting like a source), and the combination most likely to be directly affected is the one that entails both those fields, and that is the heat flux .

3. Modified SOC model

For ease of reference, the algebraic SOC model by C02 is summarized in appendix A. The model employed the state-of-the-art turbulence closure schemes for the pressure–velocity and pressure–temperature correlations. Due to its more complete pressure correlations, the model improved on the Mellor and Yamada model (MY74; MY82) in several aspects and thus provided a base for further modeling development to include both the buoyancy and shear effects on the time scale τ. In the model, the time scale τ normalized (as are all other time scales) to the dynamical time scale τ, was written as

 
formula

where λ5 is a model constant whose value was 11.04. In our new model, (10) is extended to include the Ri effect,

 
formula

The algebraic expressions for the various fluxes that exhibit the factor r are discussed in the next section.

4. Model solutions

In appendices A and B we present the general equations of the SOC model. As one can see, the equations in appendix B are still implicit. Solving them yields the following results in terms of the large-scale fields only:

  • heat flux: 
    formula
  • momentum fluxes: 
    formula

where the expressions for the momentum and heat diffusivities are given by

 
formula
 
formula
 
formula
 
formula

It is then an easy matter to obtain the expressions for all the remaining turbulence variables. For computational efficiency, we write the coefficients in (13b, c) so as to exhibit their dependence on r:

 
formula

Table 1 lists the coefficients of these polynomials, which are derived in appendix C from the basic model constants λ’s, which are the same as in C02 and are listed in Table 2.

Table 1.

Constant coefficients in Eq. (13e).

Constant coefficients in Eq. (13e).
Constant coefficients in Eq. (13e).
Table 2.

Basic model constants.

Basic model constants.
Basic model constants.

5. Consistency between (7b) and (7c)

We begin with the so-called equilibrium model, which is based on the local assumption that production equals dissipation. In this case, Eq. (3) becomes the algebraic relation

 
formula

Substituting (13b) and (13c) into (14a) yields an equation for GM as a function of Ri:

 
formula

where

 
formula

Using the Ri dependence of the coefficients c’s as from Eqs. (14c) and (13e), as well as the values given in Table 1, Eq. (14b) can be solved to yield GM versus Ri. For large Ri, we have obtained the asymptotic value

 
formula

Thus, as Ri ≫ 1, we have

 
formula

which makes (7c) compatible with (7b).

6. Can an SOC model yield σt ∼ Ri for Ri > 1?

A great deal of data indicate a linear behavior of the turbulent Prandtl number versus Ri for Ri > 1, and the derivation of such a behavior was the central goal of recent theoretical work (Galperin et al. 2007). Here we show analytically that simply by using relation (11), such a behavior is indeed recovered from a SOC model.

a. Old C02 model

We begin by considering that it was not clear to us why the standard SOC models such as C02 lead to a critical Richardson number of order 1 or less. Here, we suggest an explanation. In the old C02 model with r = λ−15 (see appendix B), the asymptotic behavior of the turbulent Prandtl number σt can be found to be

 
formula

where GH → ∞ corresponds to the limit in which turbulence ceases to exist and Ricr ≡ Ri(cr) is the corresponding critical Richardson number. On the other hand, we have a second relation,

 
formula

where Rf is the flux Richardson number with a well-known critical value Rfcr ≃ 0.25. Unlike Ri, which is a large-scale variable, Rf is a turbulence variable, and its finite critical value does not cause conceptual or mathematical problems. Using C02’s constant s values (their Table 2) in (15a), we can simultaneously solve for Ricr and σt from (15a) and (15b) to obtain

 
formula

Thus, the C02 model is not consistent with the experimental fact that σt increases with Ri linearly for Ri up to 100.

b. New model

Instead of the old C02 model result (15a), the present new model, with (11), (13ae), and (14d), yields

 
formula

The value of c = 4 in (16a), predicted by the new model, falls between the value c = 5 obtained by Zilitinkevich et al. (2007) and c = 3, preferred by Esau and Grachev (2007). The function σt(Ri) produced by the present model matches the data quite well for all Ri up to Ri = 100 (Fig. 7).

For completeness, we also derive the asymptotic behavior of SM and SH in the new model as Ri → ∞:

 
formula
 
formula
 
formula

We notice that Zilitinkevich et al. (2007) parameterize and terms in the and equations as part of the “effective dissipation rate,” thus effectively canceling the and terms for all Ri, resulting in an SM in which both the numerator and the denominator are linear in GH, and an SH in which the numerator is independent of GH while the denominator is linear in GH. Thus, σt = SM/SH → ∞, as GH → ∞, accommodating an infinitely large Ricr. In comparison, in our model these and terms vanish only for very large Ri. For small Ri, our new model behaves essentially the same as the original C02 model, and thus we do not need to introduce the additional changes in the model that Zilitinkevich et al. (2007) made.

7. Results

In Figs. 1 –6, we plot the dimensionless variables:

 
formula

The dimensionless function AM,H originates from writing the diffusivities (13a) in terms of a dissipation length scale ℓ as follows:

 
formula

Here ℓ is such that for small z, ℓ = κz, where κ = 0.4 is the von Kármán constant and B1 = 19.3.

Fig. 1.

The new dimensionless structure function for momentum, SM [Eqs. (13b) and (13c), solid line], vs Ri compared to the previous SOC model (C02, dashed line).

Fig. 1.

The new dimensionless structure function for momentum, SM [Eqs. (13b) and (13c), solid line], vs Ri compared to the previous SOC model (C02, dashed line).

Fig. 6.

Same as in Fig. 1 but for the dimensionless function AH defined in the last relation of Eq. (17a).

Fig. 6.

Same as in Fig. 1 but for the dimensionless function AH defined in the last relation of Eq. (17a).

8. Comparison with the data

Following Zilitinkevich et al. (2007, 2008), we compare our new model results with experimental, DNS, and LES data. In Figs. 7 –11 we plot the quantities

 
formula

against a variety of data using the present model (solid line) and the C02 model (dashed line).

Fig. 11.

Same as in Fig. 9 but for the dimensionless w-variance defined in Eq. (17c). The DNS data of Stretch et al. (2001) are shown by five-pointed stars. The model results and the data are not normalized to their values at Ri = 0.

Fig. 11.

Same as in Fig. 9 but for the dimensionless w-variance defined in Eq. (17c). The DNS data of Stretch et al. (2001) are shown by five-pointed stars. The model results and the data are not normalized to their values at Ri = 0.

9. Comments about the data at Ri > 1

The data for Ri > 1 require some comments. First, one may ask whether the large Ri data are observed frequently and continually to exclude contamination with spuriously large Ri due to, for example, occasional short-term reduction of local velocity gradients. Figure 12 gives an idea of the data distribution in DATABASE64 (LES), the Ohya (2001) wind tunnel dataset, and the Surface Heat Budget of the Arctic (SHEBA) and the Cooperative Atmosphere–Surface Exchange Study-1999 (CASES-99) field databases. In the LES and CASES-99 data, the fraction with Ri > 1 is significant, but in the SHEBA data only a small fraction exhibits large Ri. This fraction, however, represents relatively long periods of supercritical stratified turbulence, as reported by Grachev et al. (2005). Thus, large Ri cannot be attributed to immanent effects of data sampling.

Fig. 12.

Data distribution over stability ranges in numerical DATABASE64 (triangles), atmospheric SHEBA (circles) and CASES-99 (downward-facing triangles), and wind tunnel Ohya (diamonds) datasets. The total number of data is different in each database.

Fig. 12.

Data distribution over stability ranges in numerical DATABASE64 (triangles), atmospheric SHEBA (circles) and CASES-99 (downward-facing triangles), and wind tunnel Ohya (diamonds) datasets. The total number of data is different in each database.

Second, one may further ask whether the large Ri data represent nonturbulent features of the flow, for example, internal waves that may not be modeled with a turbulence closure approach. This is a difficult question since we lack a model for the wave–turbulence transition that could help us differentiate the appropriate modeling tools. In such a situation, we can only reverse the problem and try to see if the data we use are contaminated by internal wave processes.

It is relatively easy to estimate the contamination by internal waves in LES data. The maximum wavelength allowed in LES is defined by the size of the computational domain. Thus, by varying the domain size, it is easy to quantify changes in turbulence statistics induced by internal wave activity. No significant changes have been found in large domain experiments set up for a few runs from DATABASE64. Waves indeed developed in larger LES domains but their energy, and hence the role in turbulent statistics, remained relatively small. This conclusion is also consistent with Fedorovich et al.’s (2004) study, in which vertical transport due to irradiance of gravity waves was found negligible.

In SHEBA, orographically generated gravity waves and katabatic flows could be excluded, as the site was located on a flat sea ice surface a few hundred kilometers from the land. Thermal generation by openings in the sea ice is possible, but reasonable agreement between SHEBA and LES statistics at large Ri suggests that the role of such contamination is not important in this dataset.

For the wind tunnel experiment of Ohya (2001), the effects of gravity waves and external, transitional disturbances can also be safely excluded, because of the size of the apparatus and the inflow conditions.

In any case, the use of normalized values, where both numerator and denominator would be equally contaminated, would probably alleviate the problem. In spite of the above disparities, these heterogeneous data yield surprisingly similar results as to the shape and even the coefficients of the Ri dependences of the variables in Figs. 7 –11. The ideal situation would be to rely on good experiments—simulations of small temperature fluctuations in strongly stratified, sheared flows—to obtain more convincing data, and these are under way.

10. Conclusions

We have presented a simple generalization of the C02 model. Although the generalization is not derived from first principles, we have tried to justify it using various arguments. It is shown that with a single change, the new model can accommodate arbitrary Richardson numbers while preserving the good qualities of the original C02 model. The new model for the buoyancy time scale is the key ingredient that allows the new model to perform well for small, medium, and large Ri. The increase of the computational expenses is minimal.

Fig. 2.

Same as in Fig. 1 but for SH [Eqs. (13b) and (13c)].

Fig. 2.

Same as in Fig. 1 but for SH [Eqs. (13b) and (13c)].

Fig. 3.

Same as in Fig. 1 but for the dimensionless function GM = (τS)2.

Fig. 3.

Same as in Fig. 1 but for the dimensionless function GM = (τS)2.

Fig. 4.

Same as in Fig. 1 but for the dimensionless function GH = (τN)2.

Fig. 4.

Same as in Fig. 1 but for the dimensionless function GH = (τN)2.

Fig. 5.

Same as in Fig. 1 but for the dimensionless function AM defined in the last relation of Eq. (17a).

Fig. 5.

Same as in Fig. 1 but for the dimensionless function AM defined in the last relation of Eq. (17a).

Fig. 8.

Same as in Fig. 7 but for flux Richardson number Rf , defined in Eq. (17c).

Fig. 8.

Same as in Fig. 7 but for flux Richardson number Rf , defined in Eq. (17c).

Fig. 10.

Same as in Fig. 9 but for the squared dimensionless heat fluxes as defined in Eq. (17c).

Fig. 10.

Same as in Fig. 9 but for the squared dimensionless heat fluxes as defined in Eq. (17c).

Fig. 9.

Same as in Fig. 7 but for the squared dimensionless turbulent momentum flux defined in Eq. (17c) compared with laboratory experiments (Ohya 2001, diamonds), LES (Zilitinkevich et al. 2007, 2008, triangles), and meteorological observations (Mahrt and Vickers 2005, squares; Uttal et al. 2002, circles; Poulos et al. 2002; Banta et al. 2002, overturned triangles). The model results and all the data have been normalized to their values at Ri = 0.

Fig. 9.

Same as in Fig. 7 but for the squared dimensionless turbulent momentum flux defined in Eq. (17c) compared with laboratory experiments (Ohya 2001, diamonds), LES (Zilitinkevich et al. 2007, 2008, triangles), and meteorological observations (Mahrt and Vickers 2005, squares; Uttal et al. 2002, circles; Poulos et al. 2002; Banta et al. 2002, overturned triangles). The model results and all the data have been normalized to their values at Ri = 0.

Acknowledgments

VMC would like to express his sincere thanks to Prof. S. S. Zilitinkevich for valuable correspondence. INE acknowledges the financial support from Norwegian Research Council.

REFERENCES

REFERENCES
Abarbanel
,
H. D.
,
D.
Holm
,
J. E.
Marsden
, and
T.
Ratiu
,
1984
:
Richardson number criterion for the nonlinear stability of 3D stratified flow.
Phys. Rev. Lett.
,
52
,
2352
2355
.
André
,
J. C.
,
G.
DeMoor
,
P.
Lacarrére
, and
R.
du Vachat
,
1978
:
Modeling the 24-hour evolution of the mean and turbulent structures of the planetary boundary layer.
J. Atmos. Sci.
,
35
,
1861
1883
.
Banta
,
R. M.
,
R. K.
Newsom
,
J. K.
Lundquist
,
Y. L.
Pichugina
,
R. L.
Coulter
, and
L.
Mahrt
,
2002
:
Nocturnal low-level jet characteristics over Kansas during CASES-99.
Bound.-Layer Meteor.
,
105
,
221
252
.
Bertin
,
F.
,
J.
Barat
, and
R.
Wilson
,
1997
:
Energy dissipation rates, eddy diffusivity, and the Prandtl number: An in situ experimental approach and its consequences on radar estimate of turbulent parameters.
Radio Sci.
,
32
,
791
804
.
Canuto
,
V. M.
,
1992
:
Turbulent convection with overshooting: Reynolds stress approach.
Astrophys. J.
,
392
,
218
232
.
Canuto
,
V. M.
, and
M. S.
Dubovikov
,
1996a
:
A dynamical model for turbulence. I: General formalism.
Phys. Fluids
,
8
,
571
586
.
Canuto
,
V. M.
, and
M. S.
Dubovikov
,
1996b
:
A dynamical model for turbulence. II: Shear-driven flows.
Phys. Fluids
,
8
,
587
598
.
Canuto
,
V. M.
, and
M. S.
Dubovikov
,
1997
:
A dynamical model for turbulence. IV: Buoyancy driven flows.
Phys. Fluids
,
9
,
2118
2131
.
Canuto
,
V. M.
, and
M. S.
Dubovikov
,
1998
:
Stellar turbulent convection. I: Theory.
Astrophys. J.
,
493
,
834
847
.
Canuto
,
V. M.
,
F.
Minotti
,
C.
Ronchi
,
R. M.
Ypma
, and
O.
Zeman
,
1994
:
Second-order closure PBL model with new third-order moments: Comparison with LES data.
J. Atmos. Sci.
,
51
,
1605
1618
.
Canuto
,
V. M.
,
Y.
Cheng
, and
A. M.
Howard
,
2001
:
New third-order moments for the convective boundary layer.
J. Atmos. Sci.
,
58
,
1169
1172
.
Canuto
,
V. M.
,
Y.
Cheng
, and
A. M.
Howard
,
2007
:
Non-local ocean mixing model and a new plume model for deep convection.
Ocean Modell.
,
16
,
28
46
.
Cheng
,
Y.
, and
V. M.
Canuto
,
1994
:
Stably stratified shear turbulence: A new model for the energy dissipation length scale.
J. Atmos. Sci.
,
51
,
2384
2396
.
Cheng
,
Y.
,
V. M.
Canuto
, and
A. M.
Howard
,
2002
:
An improved model for the turbulent PBL.
J. Atmos. Sci.
,
59
,
1550
1565
.
Cheng
,
Y.
,
V. M.
Canuto
, and
A. M.
Howard
,
2005
:
Nonlocal convective PBL model based on new third- and fourth-order moments.
J. Atmos. Sci.
,
62
,
2189
2204
.
Esau
,
I. N.
, and
A. A.
Grachev
,
cited
.
2007
:
Turbulent Prandtl number in stably stratified atmospheric boundary layer: Intercomparison between LES and SHEBA data. [Available online at http://ejournal.windeng.net/16/01/Esau_Grachev_manuscript_published.pdf.]
.
Fedorovich
,
E.
,
R.
Conzemius
, and
D.
Mironov
,
2004
:
Convective entrainment into a shear-free, linearly stratified atmosphere: Bulk models reevaluated through large eddy simulations.
J. Atmos. Sci.
,
61
,
281
295
.
Galperin
,
B.
,
S.
Sukoriansky
, and
P. S.
Anderson
,
2007
:
On the critical Richardson number in stably stratified turbulence.
Atmos. Sci. Lett.
,
8
,
65
69
.
Gerz
,
T.
,
U.
Schumann
, and
S. E.
Elghobashi
,
1989
:
Direct numerical simulation of stratified homogeneous turbulent shear flows.
J. Fluid Mech.
,
200
,
563
594
.
Grachev
,
A. A.
,
C. W.
Fairall
,
P. O. G.
Persson
,
E. L.
Andreas
, and
P. S.
Guest
,
2005
:
Stable boundary-layer scaling regimes: The Sheba data.
Bound.-Layer Meteor.
,
116
,
201
235
.
Kondo
,
J.
,
O.
Kanechika
, and
N.
Yasuda
,
1978
:
Heat and momentum transfers under strong stability in the atmospheric surface layer.
J. Atmos. Sci.
,
35
,
1012
1021
.
Launder
,
B. E.
,
G.
Reece
, and
W.
Rodi
,
1975
:
Progress in the development of a Reynolds-stress turbulent closure.
J. Fluid Mech.
,
68
,
537
566
.
Lumley
,
J. L.
,
1978
:
Computational modeling of turbulent flows.
Adv. Appl. Mech.
,
18
,
123
176
.
Mack
,
S. A.
, and
H. C.
Schoeberlein
,
2004
:
Richardson number and ocean mixing: Towed chain observations.
J. Phys. Oceanogr.
,
34
,
736
754
.
Mahrt
,
L.
, and
D.
Vickers
,
2005
:
Boundary layer adjustment over small-scale changes of surface heat flux.
Bound.-Layer Meteor.
,
116
,
313
330
.
Martin
,
P. J.
,
1985
:
Simulation of the mixed layer at OWS November and Papa with several models.
J. Geophys. Res.
,
90
,
903
916
.
Mellor
,
G. L.
, and
T.
Yamada
,
1974
:
A hierarchy of turbulence closure models for planetary boundary layers.
J. Atmos. Sci.
,
31
,
1791
1806
.
Mellor
,
G. L.
, and
T.
Yamada
,
1982
:
Development of a turbulence closure model for geophysical fluid problems.
Rev. Geophys. Space Phys.
,
20
,
851
875
.
Ohya
,
Y.
,
2001
:
Wind-tunnel study of atmospheric stable boundary layers over a rough surface.
Bound.-Layer Meteor.
,
98
,
57
82
.
Poulos
,
G. S.
, and
Coauthors
,
2002
:
CASES-99: A comprehensive investigation of the stable nocturnal boundary layer.
Bull. Amer. Meteor. Soc.
,
83
,
555
581
.
Rehmann
,
C. R.
, and
J. R.
Koseff
,
2004
:
Mean potential energy change in stratified grid turbulence.
Dyn. Atmos. Oceans
,
37
,
271
294
.
Rotta
,
J. C.
,
1951
:
Statistische Theorie nichthomogener Turbulenz.
Z. Phys.
,
129
,
547
572
.
Schumann
,
U.
,
1991
:
Subgrid length-scales for large-eddy simulation of stratified turbulence.
Theor. Comput. Fluid Dyn.
,
2
,
279
290
.
Shih
,
T-H.
, and
A.
Shabbir
,
1992
:
Advances in modeling the pressure correlation terms in the second moment equations.
Studies in Turbulence, T. B. Gatski, S. Sarkar, and C. G. Speziale, Eds., Springer-Verlag, 91–128
.
Strang
,
E. J.
, and
H. J. S.
Fernando
,
2001
:
Vertical mixing and transports through a stratified shear layer.
J. Phys. Oceanogr.
,
31
,
2026
2048
.
Stretch
,
D. D.
,
J. W.
Rottman
,
K. K.
Nomura
, and
S. K.
Venayagamoorthy
,
2001
:
Transient mixing events in stably stratified turbulence. Proc. 14th Australasian Fluid Mechanics Conf., Adelaide, Australia, 612–628 pp
.
Sukoriansky
,
S.
,
B.
Galperin
, and
I.
Staroselsky
,
2005
:
A quasinormal scale elimination model of turbulent flow with stable stratification.
Phys. Fluids
,
17
,
085107
.
1
28
.
Uttal
,
T.
, and
Coauthors
,
2002
:
Surface heat budget of the Arctic Ocean.
Bull. Amer. Meteor. Soc.
,
83
,
255
275
.
Weinstock
,
J.
,
1978
:
On the theory of turbulence in the buoyancy subrange of stably stratified flows.
J. Atmos. Sci.
,
35
,
634
649
.
Zeman
,
O.
, and
J. L.
Lumley
,
1979
:
Buoyancy effects in entraining turbulent boundary layers: A second-order closure study.
Turbulent Shear Flows, F. Durst et al., Eds, Vol. 1, Turbulent Shear Flows, Springer-Verlag, 295–306
.
Zilitinkevich
,
S.
, and
I. N.
Esau
,
2007
:
Similarity theory and calculation of turbulent fluxes at the surface for the stably stratified atmospheric boundary layers.
Bound.-Layer Meteor.
,
125
,
193
206
.
Zilitinkevich
,
S. S.
,
T.
Elperin
,
N.
Kleeorin
, and
I.
Rogachevskii
,
2007
:
Energy- and flux-budget (EFB) turbulence closure model for the stably stratified flows. Part I: Steady-state, homogeneous regimes.
Bound.-Layer Meteor.
,
125
,
167
192
.
Zilitinkevich
,
S. S.
,
T.
Elperin
,
N.
Kleeorin
,
I.
Rogachevskii
,
I. N.
Esau
,
T.
Mauritsen
, and
M. W.
Miles
,
2008
:
Turbulence energetics in stably stratified geophysical flows: Strong and weak mixing regimes.
Quart. J. Roy. Meteor. Soc.
,
in press
.

APPENDIX A

Algebraic Equations of the C02 Model

As stated in C02, the state-of-the-art pressure–velocity and pressure–temperature correlations are as follows:

 
formula
 
formula

where

 
formula

Using (A1) and (A2) and neglecting the lhs of the dynamic Eqs. (2), (4), and (5), C02 obtained the algebraic model

 
formula
 
formula
 
formula

where

 
formula

The values of the λ’s are given in Table 2. In deriving (A4)(A6), the following expressions for the dissipation terms have been used:

 
formula

APPENDIX B

Flux Component Equations of the New Model

The new model simply replaces the time-scale ratio τpθ/τ = λ−15 as in the C02 model by r as defined in (11).

  1. Heat fluxes: 
    formula
     
    formula
     
    formula
  2. Temperature variance: 
    formula
  3. Momentum fluxes: 
    formula
     
    formula
     
    formula
  4. Kinetic energy components: 
    formula
     
    formula
     
    formula

The above equations still require the specification of the turbulence kinetic energy K, which is the solution of Eq. (3).

APPENDIX C

The Coefficients in Eq. (13e)

The constants in Eq. (13e) are related to λ’s in Table 2 as follows:

 
formula

These derived constants are calculated only once, and their values are listed in Table 1.

Footnotes

Corresponding author address: V. M. Canuto, NASA GISS, 2880 Broadway, New York, NY 10025. Email: vcanuto@giss.nasa.gov

* This paper is dedicated to Prof. Sergej S. Zilitinkevich on the occasion of his 70th birthday.