## Abstract

Stably stratified shear-driven turbulence is analyzed using the gradient Richardson number, Ri, as the stability parameter. The method overcomes the statistical problems associated with the widely used Monin–Obukhov stability parameter. The results of the Ri-based scaling confirm the presence of three regimes: the weakly and the very stable regimes and the transition in between them. In the weakly stable regime, fluxes scale in proportion with variance, while in the very stable regime, stress and scalar fluxes behave differently. At large Ri, the velocity field becomes highly anisotropic and the turbulent potential energy becomes approximately equal to half of the turbulent kinetic energy. It appears that even in the strongly stable regime, beyond what is known as the critical gradient Richardson number, turbulent motions are present.

## 1. Introduction

Near the earth’s surface turbulence dominates atmospheric flows in the planetary boundary layer on time scales shorter than about 1 h. Turbulent motions, however, occur in most of the atmosphere (Nastrom et al. 1984) in baroclinic flows and in conjunction with mesoscale phenomena. Flows are stably stratified when dense air is found below potentially lighter air, as is the case for most parts of the atmosphere. Stably stratified boundary layers (SBLs) occur, for instance, during night over land or sea ice when heat is lost at the surface due to imbalanced radiation, or when warm air is advected over a colder surface. In the SBL, turbulence is generated by wind shear.

Close to the ground, turbulence is the main factor controlling the vertical exchange of momentum, heat, and trace substances. This holds true outside the viscous sublayer, where molecular diffusion dominates. Strong stable stratification serves to suppress turbulent motion, hereby reducing the magnitude of the turbulent fluxes. This reduction of the fluxes lowers the SBL height, which may be as shallow as a few meters. These very shallow SBLs can accommodate very low temperatures during clear skies, permitting the formation of fogs and the rapid buildup of contaminants, and they manifest a wealth of variability and phenomena (e.g., Mahrt 1999; Poulos et al. 2002).

Observational studies of turbulence have used different stability parameters to determine the extent to which the flow is dominated by stratification. Most studies are based on the Monin–Obukhov length, *L*, which is a function of the turbulent fluxes:

where *g* is gravity, *k* is the von Kármán constant, *θ* is the mean potential temperature, *wθ* is the vertical turbulent potential temperature flux, and *τ* = −*uw* is the turbulent stress, here taken in the direction of the wind. This means that observing negative stress is possible. Note that often the turbulent stress is defined by −*uw* multiplied by the density of the air or the magnitude of the total stress vector, −, regardless of direction. The *L* was derived by Obukhov (1946) and the significance of the related stability parameter *z*/*L*, where *z* is the height above ground, was verified from observations by Monin and Obukhov (1954). The *z*/*L* can either be evaluated close to the surface (e.g., Businger et al. 1971) or locally throughout the SBL (Nieuwstadt 1984; Mahrt et al. 1998; Pahlow et al. 2001; Klipp and Mahrt 2004; Grachev et al. 2005). Other investigators use the mean profile of wind and potential temperature to form stability parameters (Kondo et al. 1978; Schumann and Gerz 1995; Poulos and Burns 2003; Klipp and Mahrt 2004). Yet other studies simply make a classification by dividing between cases of strong and weak turbulence (Sorbjan 2006; Mahrt and Vickers 2006), or between continuous, intermittent, and radiatively dominated turbulent regimes (Van de Wiel et al. 2003), deliberately avoiding the use of stability parameters when examining the stably stratified turbulence.

Studies by Businger et al. (1971) and Nieuwstadt (1984) have made significant contributions to our current understanding of the weakly stratified SBL (Sorbjan 2006). Monin–Obukhov-based studies result in semiempirical relations between mean profiles of wind, temperature and tracers, and their respective turbulent fluxes. These relations are widely used when modeling the effects of turbulence in the atmospheric boundary layer (e.g., Louis 1979; Delage 1997). The use of *z*/*L* as a stability parameter is, however, complicated by the fact that the method suffers from self-correlations. When using the same variables on both axes, random errors can generate an unphysical self-correlation. In unfavorable cases this may lead to artificial power-law relations obscuring the results (Hicks 1978; Mahrt et al. 1998; Andreas and Hicks 2002; Klipp and Mahrt 2004; Baas et al. 2006). For instance, a plot intending to show the stability dependence of *τ* will be affected because *τ* itself appears in the denominator of *z*/*L*. In this case the turbulent stress will be overestimated at small *z*/*L* and underestimated at large *z*/*L*.

When determining to which extent the flow is dominated by stratification, an alternative to using *z*/*L* is to evaluate the nondimensional gradient Richardson number (Richardson 1920). This was proposed for instance by Klipp and Mahrt (2004) as well as Sorbjan (2006) using similarity arguments. The gradient Richardson number is defined as

where *g* is gravity, *θ* is potential temperature, *z* is the height above ground, *U* and *V* are the mean wind components, *N* is the Brunt–Väisälä frequency and *S* is the total vertical shear. In practice, Ri is calculated from finite differences between instrument levels, as to be described later.

In a theoretical study of the turbulent kinetic energy budget equation, Richardson (1920) found that the limit Ri_{cr} = 1, determines the onset of turbulence from laminar flow. The Ri_{cr} is known as the critical Richardson number. Chandrasekhar (1961) used a package exchange method to arrive at a similar result (see Miles 1986), while Miles (1961) and Howard (1961) employed linear stability analysis of infinitesimal disturbances in the Navier–Stokes equations to arrive at Ri_{cr} = ¼. Observations of atmospheric flows, however, indicate that turbulence exists above Ri_{cr}, both in the SBL (Kondo et al. 1978) and in the stably stratified free atmosphere (Tjernström 1993).

In the present observational study of stably stratified shear-driven atmospheric turbulence, we use estimates of Ri as the basic stability parameter. The method avoids the self-correlation problems inherent to studies based on *z*/*L*. From the analysis, we find that we can separate data into three regimes, similar to what Mahrt et al. (1998) found using *z*/*L*:

In the “weakly stable regime,” where Ri <

*O*(0.1), the fluxes are proportional to the variance (e.g., enhanced temperature fluctuations are associated with a stronger downward heat flux). The turbulent energy is dominated by kinetic energy.The “transition regime,”

*O*(0.1) < Ri <*O*(1), is characterized by a rapid decrease in normalized fluxes with increasing stability.The “very stable regime,” when Ri >

*O*(1), has a very different character. Here the momentum flux is small but finite, while it appears that the heat flux is very small, possibly zero. A large fraction of the turbulent energy is potential.

The goal of this study is to investigate the effects of stable stratification on selected properties of atmospheric turbulence. We analyze data from six datasets, collected in very different physical environments. The strength in using several different datasets is that we can identify common features where correspondence occurs, while disagreement points to the need for further investigation. The results of the analysis may be useful in a wide range of applications. For example, they could be used as closure assumptions for parameterizations of the stably stratified atmospheric boundary layer as an alternative to Monin–Obukhov-based parameterizations.

## 2. Data analysis

The datasets used in this study were made available to the community by the investigators responsible for conducting the measurements. The choice of datasets was made primarily on the basis of practical considerations, such as the possibility of estimating Ri and to provide a large portion of stably stratified data, including cases with very strong stability. Furthermore, the intention was to include a variety of physical environments. The data include observations from grassland [the Cooperative Atmosphere–Surface Exchange Study-99 (CASES-99)], snow-covered land [the Fluxes over Snow Surfaces, phases 1 and 2 (FLOSS) and (FLOSS-II), respectively], sea ice [the Surface Heat Budget of the Arctic Ocean (SHEBA)], mountain–forest canopy [the Carbon in the Mountains Experiment (CME)], and sea as observed from an island in the Baltic Sea [Östergarnsholm (OGH)]. Details of the experiments are summarized in Table 1 and schematically shown in Fig. 1.

The datasets were prepared by the investigators who conducted the measurements. This is a complicated process, which requires specialized expertise and good knowledge of the instrumentation. A crucial step is the separation of the raw signals into a mean and a turbulent part, *U* = *U* + *u*. This is in practice accomplished by assuming stationarity and time averaging the signal to estimate the mean, *U*, whereafter the mean is subtracted from the raw data to yield the turbulent part. Then the second moments (i.e., fluxes and variances) can be calculated (e.g., *wθ*). Alternative methods include Fourier and wavelet decompositions (e.g., Howell and Mahrt 1997; Vickers and Mahrt 2003; Cornish et al. 2006).

The choice of the averaging window is a compromise. If it is chosen too short, not all the relevant scales are included. On the other hand, if it is too long the assumption of stationarity is likely to be violated. The fluxes can be divided into contributions from both microscale turbulence and mesoscale phenomena (e.g., Finnigan and Einaudi 1981). Vickers and Mahrt (2003, 2006) found that it is sometimes necessary to make initial averages as short as 1–2 min in order to separate these two scales. This was done by assuming the presence of a spectral gap at the time scales where the heat flux is close to zero. However, Tjernström (2005) analyzed velocity spectra and found that clear spectral gaps were only present 22% of the time during the period under consideration. If no spectral gap is present the calculated turbulent fluxes and variances will depend on the averaging window (Duynkerke 1998; Mahrt et al. 1998). In these cases the Reynolds averaging method allows contamination by leakage from scales larger than the averaging window. Furthermore, different averaging windows could lead to varying results, even in the absence of leakage, depending on the shape of the spectra and cospectra for the variables.

In the present study, we choose to use the flux estimates provided by the investigators responsible for the experiments, including their instrumental calibrations and quality checks. For CASES-99, CME, and FLOSS the averaging time was 5 min. SHEBA had an averaging time of 10 min for estimating the turbulent moments, whereafter these values were averaged over an hour, while for OGH, a 10-min running mean was used over 1-h intervals. Finally, for FLOSS-II the variable window method of Vickers and Mahrt (2006) was applied to calculate the turbulence statistics. While the use of varying analysis methods in the present study is considered a weakness, we have no other practical possibilities for the time being. There is a strong need for high-quality datasets, which are analyzed in the same way.

The stability parameter, Ri, is here estimated as a finite-difference approximation from wind and temperature measurements above (subscript *a*) and below (subscript *b*) the sensor where the turbulent fluxes and variances are observed, leading to the so-called bulk- or finite-difference Richardson number:

Figure 1 shows the instrumentation on the towers in the different experiments. For example, for SHEBA fluxes measured at 3 m above ground, gradients were calculated between the 2- and 5-m levels. The difference *z _{a}* −

*z*varied between 3 and 10 m, unless otherwise stated.

_{b}As noted by Arya (1991), taking a discrete difference to approximate the gradient is questionable when *z _{a}*/

*z*is large, because real profiles are not linear close to the ground. It is therefore important to use values of wind and temperature at some distance above ground (e.g., Vogelezang and Holtslag 1996). On the other hand Arya (1991) showed that linear finite differences are superior over logarithmic differences even for large

_{b}*z*/

_{a}*z*(2 and 4 were tested) in the very stable regime, which is the focus in the current study. Further, evaluating finite differences over shorter distances than what we have done here may be questioned, because errors in the measurements may have a larger impact when two nearly equal mean values are subtracted.

_{b}In general the datasets were of high quality. However, some selection of data was necessary. The wind speed was required to increase with height. This is to avoid complex wind profiles, for instance, when a low-level jet is present below the sensor in which case the stress is expected to reverse, which was the case in 8% of the observations. Including these cases in the analysis increased the scatter among the datasets, but did not affect the overall conclusions of the present study. For the heat flux, we had to impose the criteria that *N* > 0.02 s^{−1} in order to avoid locally stable, but nonlocally unstable cases. These can for example occur in a growing convective boundary layer, when the ground is heating rapidly. Though this is not a perfect criterion, it is the only generally applicable one we could find. In individual datasets, these cases could be identified with independent instruments placed close to the ground, though such an approach would give different criteria for different datasets. In the analysis, we removed data points where the variance was identically zero, since they were statistically overrepresented, probably due to the analysis algorithms. For OGH, we only included data when the wind came from the open sea as described in Smedman et al. (2003).

We divided all data in bins of Ri, equally spaced on a logarithmic axis. For each bin, we calculated the 95% confidence interval using a two-sided *t* test. The mean value in this test is in the middle of the confidence interval, and hence, we only plotted the confidence intervals in most cases in order to avoid cluttered graphs. Data is most abundant at Ri ∼ 0.1, while for very small and very large Ri fewer data are available. This gives the confidence intervals a “trumpetlike” shape. We included only bins with more than 10 values.

The proposed approach is inherently “*z* less,” as the height above ground is not included in the stability parameter. While the individual sensors are fixed on the masts, the SBL may vary dramatically in depth, such that on different occasions the sensors are in the surface layer, while at other times may be situated above the boundary layer.

The large amount of data allows us to show more confidence in the results, than we could have with a single dataset. The six datasets provide us with a total equivalent of 5.5 yr of usable data, excluding data failing quality control and convectively unstable data. The Ri was larger than unity for 8% of the data, corresponding to more than five equivalent months. The use of results from different environments, different instruments, and data prepared with different algorithms allow us to identify common features and differences.

## 3. Results

Here we present several properties of stably stratified shear-driven turbulence as a function of Ri. We emphasize that all of these properties can be displayed while avoiding self-correlations between the stability parameter and the entity of interest. Furthermore, common practice in studies of this kind is to nondimensionalize variances by fluxes (Nieuwstadt 1984; Mahrt et al. 1998; Pahlow et al. 2001). However, as we shall see, this may be problematic for the very stable regime, where the fluxes are small or possibly zero. On the contrary we suggest to normalize the turbulent fluxes of momentum, heat, and tracers by variances.

### a. Turbulent stress

First we are interested in the turbulent stress. The turbulent stress is the momentum flux in the opposite direction of the wind. We normalize the stress either by the turbulent kinetic energy to form the nondimensional stress:

or by horizontal and vertical variances to form the stress correlation coefficient:

where *τ* = −*uw* is the stress, and the turbulent kinetic energy *E* = ½(*σ*^{2}_{u} + *σ*^{2}_{υ} + *σ*^{2}_{w}), is half the sum of the squared variances of the three wind components. Here cor* _{τ}* is an approximation to the correlation coefficient, since we used the average horizontal variance, rather than variance in the direction of the wind. The denominator of (5) is dependent on the zero-plane rotation (e.g., Wilczak et al. 2001), while for (4) it is not. Therefore the latter is more trustworthy when signals are weak. However, the vertical flux estimates are still dependent on the zero-plane rotation.

The normalized stresses are shown in Fig. 2. The general behavior of the stress with the two normalizations is very similar. We see that most datasets exhibit a nearly constant value of *f _{τ}* in the near-neutral regime (Ri < 0.1) of 0.15–0.19 and 0.35–0.4 for cor

*. We note that Nieuwstadt (1984) found 2*

_{τ}*E*/

*τ*≈ 2.7 − 3.0 corresponding to

*τ*/

*E*≈ 0.22–0.27 for near-neutral conditions, which differs considerably from the values found in the present study. Mahrt et al. (1998) found (

*σ*

^{2}

_{u}+

*σ*

^{2}

_{υ})/

*τ*≈ 3 and

*σ*/

_{w}*τ*≈ 1 implying

*τ*/

*E*≈ 0.20 (i.e., closer to the values found here). Possible explanations for this discrepancy are, the use of inverse stresses and a different definition of

*τ*. In the former case it is clear that the mean is not conserved under inversion, while in the latter case however, it is not clearly stated in the mentioned papers how the calculation was done. We also note that in both these studies data was plotted in such a manner that self-correlation cannot be avoided. Self-correlation will lead to larger estimates since random positive stress errors permits lower values of

*z*/

*L*(i.e., overestimated near-neutral values). We found self-correlation to explain the difference up to the level found by Mahrt et al. (1998), by making a plot (not shown) with

*z*/

*L*instead of Ri with the same datasets.

Scaled stress decreases in the transition regime (0.1 < Ri < 1) to a new level in the very stable regime (Ri > 1) which varies among the datasets. For all datasets, we see that there is a mean nonzero stress even in the very stable regime. FLOSS data have a 95% confidence interval, which includes zero for *f _{τ}*, while SHEBA and FLOSS confidence intervals for cor

*include zero. Therefore, we cannot conclude that they are nonzero, though the mean values in all bins are positive. As noted above,*

_{τ}*f*is preferable over cor

_{τ}*when signals are weak. Since all binned means are positive it is very likely that further sampling would lead to more narrow confidence intervals that does not include zero. OGH exhibits larger values, of normalized stress than the other datasets at all Ri and SHEBA likewise for Ri < 0.1. This behavior could be attributed to the longer averaging window used, as more scales are included. However, CME contradicts this hypothesis in the very stable regime where it has a larger normalized stress than SHEBA. The different terrains, consisting of comparably rough mountain forest–canopy in CME and of relatively smooth sea ice in the case of SHEBA, could explain the behavior of CME at large Ri. A simple interpolation between two regimes was chosen, which is centered at Ri = 0.25. The empirical curves shown in Fig. 2 are*

_{τ}*f*(Ri) ≈ 0.17 [0.75/(1 + 4 Ri) + 0.25] and cor

_{τ}*(Ri) ≈ 0.37 [0.75/(1 + 4Ri) + 0.25].*

_{τ}When calculating the 95% confidence interval it is assumed that the probability distribution of the data follows the Student’s *t* distribution. When there are many data points this distribution approaches the normal distribution according to the central limit theorem. Figure 3 shows the probability distributions of *τ*/*E* for the three stability regimes including all datasets. For each regime, we graphed the corresponding normal distribution with the same mean and standard deviation. Figure 3 shows that the probability distributions are nearly normally distributed. The stronger the stability, the wider the probability distribution. As a consequence, a large part of the distribution in the very stable regime is negative. We believe this is related to random errors made when estimating the fluxes, which with the presented approach is simply handled as an error. In Monin–Obukhov scaling it is assumed a priori that the stress should be positive.

### b. Turbulent heat flux

Next, we examine the nondimensional heat flux and correlation coefficient, which we define in a way analogous to the stress:

where *wθ* is the vertical potential temperature flux and *σ*^{2}_{θ} is squared potential temperature variance. The results are shown in Fig. 4. The probability distributions (not shown) are multimodal due to the lack of agreement among the datasets. Again the results for the two normalizations does not differ significantly.

The weakly stable (0.01 < Ri < 0.1) regime values of normalized heat flux varies considerably among the datasets. This could be due to the dataset preparation, since three of the datasets that agree well with each other (CASES-99, CME, and FLOSS) all happen to be prepared using the same algorithm based on a 5-min averaging window. Another possibility is bias in the instruments determining the potential temperature gradient, causing different degrees of contamination of data with cases of unstable stratification. Further, the signal-to-noise ratio from the fast-response temperature sensors decreases as neutral stratification is approached. This is both due to instrumental noise as well as horizontal temperature advection. Noise contributes to the variance but is unlikely to correlate with *w*, causing *f _{θ}* and cor

*to decrease. This gives a bias toward zero as Ri approaches zero, making the larger values from SHEBA more trustworthy. Since we required that*

_{wθ}*N*> 0.02 s

^{−1}fewer values are available in the weakly stable regime, making the curves end earlier.

In the very stable regime five of the datasets have confidence intervals that include zero, whereas the FLOSS-II dataset exhibits positive heat fluxes in this regime. This could possibly be due to horizontal heteorogeneity, generated by partial snow cover, causing the fluxes to be out of equilibrium with the mean flow gradients (Mahrt and Vickers 2005). Another possibility is problems with measuring the temperature in weak-wind situations, due to the lack of ventilation (L. Mahrt 2006, personal communication). The results from the five other datasets indicate that we cannot reject the hypothesis that the heat flux is zero in the very stable regime. As a consequence, Monin–Obukhov similarity scaling fails in the very stable regime (cf. Sorbjan 2006).

Empirical curves are suggested to be *f _{wθ}*(Ri) ≈ −0.145/(1 + 4Ri) and cor

*(Ri) ≈ −0.3/(1 + 4Ri), where the values −0.145 and −0.3 are rather uncertain. Both are taken from the normalized water vapor flux to be determined in section 3c. In the weakly stable regime these curves do not fit the data very well. However, as noted above, this could be due to contamination of data subject to unstable stratification and low signal-to-noise ratios. We may consider potential temperature as a passive tracer in the near-neutral limit, which follows the turbulent flow without altering it. Then simple reasoning gives*

_{wθ}*wθ*∝

*σ*. It follows that for positive Ri tending to zero, cor

_{w}σ_{θ}*should approach a constant.*

_{wθ}### c. Turbulent water vapor flux

The generalized nondimensional tracer flux and correlation coefficients are

where *wc* is the vertical tracer flux and *σ*^{2}_{c} is the squared tracer variance. Figure 5 shows, as an example, an application to the water vapor flux measured at 5 m above ground during CASES-99. After data selection, an accumulated total of 16 equivalent days was included in the analysis. In addition to the 95% confidence intervals, both mean and median values, along with standard deviations are shown. The median values are less sensitive to outliers, and in a skewed distribution it will be closer to the mode than the mean value.

It is important to keep track of the mean scalar gradient when applying (8) and (9). This was done by inspecting the difference between values observed at 5- and 20-m heights. When both measurements were available and Ri was positive the water vapor concentration was largest close to the ground. As a result the normalized fluxes for the weakly stable regime were positive and almost constant over a wide range. This was because the campaign occurred in very dry conditions. Under more moist conditions dew deposition may occur.

The very stable regime exhibits quite wide confidence intervals, due to the small amount of data. As for the case of heat flux (section 3b), the null hypothesis of zero water vapor flux at large Ri cannot be rejected. Furthermore, median values are closer to zero than the mean values, indicating that the mean values are affected by a few positive spikes in the data. The method could be extended to include more datasets and other variables, which can be measured with good precision at high rates. The empirical curves are *f _{c}*(Ri) ≈ 0.145/ (1 + 4 × Ri) and cor

*(Ri) ≈ 0.3/(1 + 4 × Ri), which have the same functional shapes as the ones suggested for*

_{c}*f*and cor

_{wθ}*. The limited data supports a somewhat sharper transition between the regimes.*

_{wθ}### d. Anisotropy

Turbulence is isotropic when the statistical properties of the flow are invariant under rotations of the coordinate system. We already know that this is not the case for shear-driven turbulence, since there is a positive stress for all Ri (Fig. 2). The degree to which the flow is anisotropic can be estimated by the ratio of the horizontal velocity variance to the vertical:

The stability dependency of this anisotropy measure is presented in Fig. 6. By definition, the error distribution of *A* is skewed because it is limited to positive values. We therefore chose to use medians, rather than means in Fig. 6. These are seen to have values in the weakly stable regime of 5–8, which is close to what Mahrt et al. (1998) and others found. This is considerably larger than for isotropic turbulence, which has *A* = 2.

In the very stable regime values of *A* are in the range of 8–30, indicating that vertical motion is suppressed by the increasing stability. Four datasets agree well with each other in the very stable regime with values of 12–15, whereas OGH and FLOSS exhibit more extreme values. Possible errors in the zero-plane rotation influences *A*, which could explain part of the spread. Note that it is not the same subset of datasets that agreed in the weakly stable limit of the nondimensional heat flux (Fig. 4). Across all datasets *A* exhibits a global minimum close to Ri = 0.1; however, we have no physical explanation for this.

### e. Turbulent potential energy

An interesting feature of stably stratified turbulence is that the energy of the turbulent part of the flow comprises both turbulent kinetic energy, *E*, and turbulent potential energy, *P*. The latter is proportional to the square of the density variations. Here *P* is related to the potential temperature variance (e.g., Schumann and Gerz 1995):

As we are interested in the Ri dependence of *P*, the problem of self-correlation appears since the results on both axes have a functional dependency on *N* ^{2}. Here, we avoided this issue by using different estimates of *N* ^{2} on each axis. This was accomplished by calculating two independent gradients over vertical distances of 10 and 20 m, respectively, centered around the flux level. We used the CASES-99 dataset (Fig. 1). Due to this restriction only flux data from 20-, 30-, and 40-m levels were included. A similar method could be used to avoid self-correlation in *z*/*L*-based studies, by placing two or more sonic anemometers at the same level.

Figure 7 shows three curves derived from the CASES-99 dataset. As in Schumann and Gerz (1995) we normalize *P* by *σ*^{2}_{w}. The outer gradients, situated 20 m apart, were used to calculate Ri, for the “outer Ri, inner *P*” curve and vice versa for the “inner Ri, outer P” curve. To demonstrate a potential effect of self-correlation, Fig. 7 also include *P*/*σ*^{2}_{w} plotted using the same inner estimate of *N* ^{2} for both Ri and *P*. The behavior of the self-correlated curve differs from the other two at small Ri. This can easily be understood, since a measuring error in *N* ^{2}, causing Ri to be small, will automatically act to increase *P*. Thus simple random scatter in the data will cause the unrealistic increase of *P*/*σ*^{2}_{w} at small Ri.

In stationary, horizontally homogeneous, and near-neutral conditions *P* tends to zero. As we can see in Fig. 7, *P*/*σ*^{2}_{w} tends to a finite value above zero, due to instrumental noise and horizontal advection, see section 3b. It is noteworthy that the transition between the two regimes is not centered around Ri = 0.25, as has been the case in the other plots, except for the self-correlated dataset. This is because the inner gradients can be larger than the outer gradients, giving a positive bias on *P* for the inner Ri, outer P curve, due to the fact that *N* is in the denominator of *P*. The opposite is the case for the outer Ri, inner *P* curve. At increasing values of Ri, the proportion of *P* to *σ*^{2}_{w} increases, and at large Ri the data indicates a value about 3–4, which is somewhat larger than the value 2.5 found by Schumann and Gerz (1995). Combining this with the anisotropy curve for CASES-99, Fig. 6, we find that *P* is about half of *E* in the very stable regime. For comparison a field of linear internal gravity waves have *P*_{wave} = *E*_{wave} (e.g., Nappo 2002).

## 4. Conclusions

We have analyzed tower-based turbulence observations from six different datasets, including grassland, snow-covered land, mountain–forest canopy, ice-covered ocean, and ocean as observed from an island in the Baltic Sea, accumulating a total of 5.5 yr worth of data of stably stratified situations. The gradient Richardson number, Ri, was used as the stability parameter, rather than the widely used Monin–Obukhov parameter, *z*/*L*. By doing this, we avoid the effects of self-correlation for the properties we are interested in. We confirm that the division of the stability range, defined by Mahrt et al. (1998), into weakly stable, transition stability, and very stable regimes, also applies within Ri-based studies.

In the weakly stable regime nondimensionalized stress, potential temperature flux, and water vapor fluxes appear to have almost constant values. All datasets exhibit a monotonic transition from the weakly stable to the very stable regime, where scaled fluxes attain new, reduced levels. While we are statistically confident that stress is nonzero in most of the six datasets, heat and scalar fluxes cannot be distinguished from zero. At the same time the turbulence becomes very anisotropic with a large dominance of the horizontal modes compared to those in the vertical and the turbulent potential energy increases to a level close to half the turbulent kinetic energy. The observed properties of atmospheric shear-driven turbulence at large Ri, with nonzero stress, but zero heat flux, are not inconsistent with internal gravity wave activity.

Both in the weakly and very stable regimes some datasets disagree with the others. There seems to be no systematic behavior of the aberrant datasets, though OGH seems to diverge more often than the rest. The reasons for this can be the 1) different physical environments, 2) different instrumentation, and 3) data treatment. The latter including zero-plane rotation and time averaging. Here OGH differ from the rest in all points, being from the marine boundary layer, using Gill sonic anemometers and utilizing a 10-min running mean. From the present study we are not able to judge the relative importance of these factors. During the preparation of this paper we had the opportunity to work with two versions of the FLOSS-II dataset. One used 5-min averages and the other variable averaging time. There were some statistically significant differences between the two versions; the rather strange upward heat flux at large Ri was stronger and in the near-neutral *τ*/*E* was closer to 0.16 in the former case. This points at the importance of data treatment, though the overall conclusions were the same independent of dataset.

The results of the present analysis have a wide range of applicability to turbulence parameterizations in atmospheric models, it may further be utilized in dispersion models of chemicals and small aerosols. It could be used to validate laboratory experiments, large eddy simulations, and direct numerical simulations of stably stratified flows. The results pose new challenges to theoretical development, due to the observed presence of turbulent motions, far beyond any Ri* _{cr}* predicted by linear theory. One might speculate that the transition to the very stable regime is nature’s way of upholding turbulence even at strong flow stability.

## Acknowledgments

We wish to thank H. Bergström, T. Horst, G. Maclean, S. Oncley, O. Persson, G. Poulos, A.-S. Smedman, J. Sun, T. Uttal, D. Vickers, and their collegues for supplying the data used in the present study and S. Zilitinkevich, B. Grisogono, G.-J. Steeneveld, I. Esau, M. Tjernström, P. Lundberg, L. Mahrt, and several others for inspiring discussions and help in preparing the manuscript. Support for this work was from the Swedish Research Council Contract 621-2001-1846.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address*: Thorsten Mauritsen, Department of Meteorology, Stockholm University, Svante Arrhenius Väg 12, Stockholm 106 91, Sweden. Email: thorsten@misu.su.se