## 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 *N* ^{2} that acts as a sink. The two combine into a single parameter, the Richardson number Ri, defined as

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 *U _{i}* fields, are as follows (see C02):

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)]:

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

C02 and many other previous SOC models further assume that

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 *τ _{pθ}* that enters the buoyancy spectrum

*B*(

*k*) = −

*τ*(

_{pθ}*k*)

*N*

^{2}

*E*(

*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

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

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

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 *h*_{3} = *h*

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; *τ _{pθ}* is a relaxation time scale, while

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

_{J}*E*(

*k*),

*E*(

_{z}*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

where *P*_{33}(*k*) = 1 − *k*^{2}_{z}/*k*^{2} 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

we derive the following expressions for *τ _{pθ}* and

*τ*

_{J}:the ratio of which is given by

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

We proceed to model the dissipation time scale *τ _{J}* as proportional to the turbulent kinetic energy dissipation time scale

*τ*:

which is of the form (7c).

We have concentrated on *τ _{pθ}* 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 .

_{pυ}## 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 *τ _{pθ}*. In the model, the time scale

*τ*normalized (as are all other time scales) to the dynamical time scale

_{pθ}*τ*, was written as

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

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:

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

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

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.

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

where

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 *G _{M}* versus Ri. For large Ri, we have obtained the asymptotic value

Thus, as Ri ≫ 1, we have

## 6. Can an SOC model yield *σ*_{t} ∼ Ri for Ri > 1?

_{t}

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* = *λ*^{−1}_{5} (see appendix B), the asymptotic behavior of the turbulent Prandtl number *σ _{t}* can be found to be

where *G _{H}* → ∞ corresponds to the limit in which turbulence ceases to exist and Ri

_{cr}≡ Ri(cr) is the corresponding critical Richardson number. On the other hand, we have a second relation,

where *R _{f}* is the flux Richardson number with a well-known critical value

*R*

_{f}_{cr}≃ 0.25. Unlike Ri, which is a large-scale variable,

*R*is a turbulence variable, and its finite critical value does not cause conceptual or mathematical problems. Using C02’s constant

_{f}*s*values (their Table 2) in (15a), we can simultaneously solve for Ri

_{cr}and

*σ*from (15a) and (15b) to obtain

_{t}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), (13a–e), and (14d), yields

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 S_{M} and S_{H} in the new model as Ri → ∞:

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 *S _{M}* in which both the numerator and the denominator are linear in

*G*, and an

_{H}*S*in which the numerator is independent of

_{H}*G*while the denominator is linear in

_{H}*G*. Thus,

_{H}*σ*=

_{t}*S*/

_{M}*S*→ ∞, as

_{H}*G*→ ∞, accommodating an infinitely large Ri

_{H}_{cr}. 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

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

Here ℓ is such that for small *z*, ℓ = *κ**z*, where *κ* = 0.4 is the von Kármán constant and *B*_{1} = 19.3.

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

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

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

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.

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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

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

where

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

where

### APPENDIX B

#### Flux Component Equations of the New Model

The new model simply replaces the time-scale ratio *τ*_{pθ}/*τ* = *λ*^{−1}_{5} as in the C02 model by *r* as defined in (11).

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)

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.