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:

 
formula

where g is gravity, k is the von Kármán constant, θ is the mean potential temperature, 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

 
formula

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 Ricr = 1, determines the onset of turbulence from laminar flow. The Ricr 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 Ricr = ¼. Observations of atmospheric flows, however, indicate that turbulence exists above Ricr, 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:

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

  2. The “transition regime,” O(0.1) < Ri < O(1), is characterized by a rapid decrease in normalized fluxes with increasing stability.

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

Table 1.

Dataset overview.

Dataset overview.
Dataset overview.
Fig. 1.

Schematic overview of the instrumentation at the experimental sites. The vertical lines are towers and the numbers indicate the height of each instrument level. Horizontal lines are levels with sonic anemometers used for turbulence measurements, while dots are levels with instruments used in the calculation of gradients of mean temperature and wind. In all cases, gradients are calculated from instrumentation above and below the level of interest. For instance, CASES-99 gradients at the 5-m level is calculated as finite differences between 1.5 and 10 m. SHEBA heights are nominal only. The actual measured heights were used in the analysis. See Table 1. 

Fig. 1.

Schematic overview of the instrumentation at the experimental sites. The vertical lines are towers and the numbers indicate the height of each instrument level. Horizontal lines are levels with sonic anemometers used for turbulence measurements, while dots are levels with instruments used in the calculation of gradients of mean temperature and wind. In all cases, gradients are calculated from instrumentation above and below the level of interest. For instance, CASES-99 gradients at the 5-m level is calculated as finite differences between 1.5 and 10 m. SHEBA heights are nominal only. The actual measured heights were used in the analysis. See Table 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., ). 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:

 
formula

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 zazb varied between 3 and 10 m, unless otherwise stated.

As noted by Arya (1991), taking a discrete difference to approximate the gradient is questionable when za/zb 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 za/zb (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.

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:

 
formula

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

 
formula

where τ = −uw is the stress, and the turbulent kinetic energy E = ½(σ2u + σ2υ + σ2w), 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 2E/τ ≈ 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 (σ2u + σ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.

Fig. 2.

Two-sided t test 95% confidence intervals on the binned mean of the nondimensionalized stresses (a) fτ and (b) corτ, as functions of Ri. Colors indicate dataset, while gray shading indicates 0.1 < Ri < 1 and the vertical dashed line is Ri = ¼. The thick dashed curves are empirical.

Fig. 2.

Two-sided t test 95% confidence intervals on the binned mean of the nondimensionalized stresses (a) fτ and (b) corτ, as functions of Ri. Colors indicate dataset, while gray shading indicates 0.1 < Ri < 1 and the vertical dashed line is Ri = ¼. The thick dashed curves are empirical.

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.

Fig. 3.

Probability distributions of the measure fτ for the three Ri regimes including all datasets. The curves are normalized. For each distribution a corresponding Gaussian distribution is plotted with thin lines.

Fig. 3.

Probability distributions of the measure fτ for the three Ri regimes including all datasets. The curves are normalized. For each distribution a corresponding Gaussian distribution is plotted with thin lines.

b. Turbulent heat flux

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

 
formula
 
formula

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

Fig. 4.

Same as in Fig. 2 but for the nondimensional potential temperature flux measures (a) f and (b) cor.

Fig. 4.

Same as in Fig. 2 but for the nondimensional potential temperature flux measures (a) f and (b) cor.

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 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(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σθ. It follows that for positive Ri tending to zero, cor should approach a constant.

c. Turbulent water vapor flux

The generalized nondimensional tracer flux and correlation coefficients are

 
formula
 
formula

where wc is the vertical tracer flux and σ2c 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.

Fig. 5.

(a) Nondimensional water vapor flux and (b) the correlation coefficient for the 5-m level of CASES-99. The plots show binned mean and median values, 95% confidence intervals on the mean as shading, and standard deviations as vertical solid lines. The thick dashed curves are empirical.

Fig. 5.

(a) Nondimensional water vapor flux and (b) the correlation coefficient for the 5-m level of CASES-99. The plots show binned mean and median values, 95% confidence intervals on the mean as shading, and standard deviations as vertical solid lines. The thick dashed curves are empirical.

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 fc(Ri) ≈ 0.145/ (1 + 4 × Ri) and corc (Ri) ≈ 0.3/(1 + 4 × Ri), which have the same functional shapes as the ones suggested for f and cor. The limited data supports a somewhat sharper transition between the regimes.

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:

 
formula

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.

Fig. 6.

Median of anisotropy as a function of Ri for the same data as in Fig. 2. The horizontal dashed black line is for isotropic turbulence.

Fig. 6.

Median of anisotropy as a function of Ri for the same data as in Fig. 2. The horizontal dashed black line is for isotropic turbulence.

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

 
formula

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 N2. Here, we avoided this issue by using different estimates of N2 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 σ2w. 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/σ2w plotted using the same inner estimate of N2 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 N2, 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/σ2w at small Ri.

Fig. 7.

Ratio of turbulent potential energy (P) to squared vertical wind variance from the CASES-99 dataset. The shaded area is the 95% confidence intervals and the lines are mean values. See the text for further details.

Fig. 7.

Ratio of turbulent potential energy (P) to squared vertical wind variance from the CASES-99 dataset. The shaded area is the 95% confidence intervals and the lines are mean values. See the text for further details.

In stationary, horizontally homogeneous, and near-neutral conditions P tends to zero. As we can see in Fig. 7, P/σ2w 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 σ2w 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 Pwave = Ewave (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 Ricr 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

REFERENCES
Andreas
,
E. L.
, and
B. B.
Hicks
,
2002
:
Comments on “Critical test of the validity of Monin–Obukhov similarity during convective conditions.”.
J. Atmos. Sci.
,
59
,
2605
2607
.
Arya
,
S. P.
,
1991
:
Finite-difference errors in estimation of gradients in the atmospheric surface layer.
J. Appl. Meteor.
,
30
,
251
253
.
Baas
,
P.
,
G. J.
Steeneveld
,
B. J. H.
van de Wiel
, and
A. A. M.
Holtslag
,
2006
:
Exploring self-correlation in flux–gradient relationships for stably stratified conditions.
J. Atmos. Sci.
,
63
,
3045
3054
.
Businger
,
J. A.
,
J. C.
Wyngaard
,
Y.
Izumi
, and
E. F.
Bradley
,
1971
:
Flux-profile relationships in the atmospheric surface layer.
J. Atmos. Sci.
,
28
,
181
189
.
Chandrasekhar
,
S.
,
1961
:
Hydrodynamic and Hydromagnetic Stability.
Clarendon Press, 652 pp
.
Cornish
,
C. R.
,
C. S.
Bretherton
, and
D. B.
Percival
,
2006
:
Maximal overlap wavelet statistical analysis with application to atmospheric turbulence.
Bound.-Layer Meteor.
,
119
,
339
374
.
Delage
,
Y.
,
1997
:
Parameterising sub-grid scale vertical transport in atmospheric models under statically stable conditions.
Bound.-Layer Meteor.
,
82
,
23
48
.
Duynkerke
,
P. G.
,
1998
:
Dynamics of cloudy boundary layers.
Clear and Cloudy Boundary Layers, A. A. M. Holtslag and P. G. Duynkerke, Eds., Royal Netherlands Academy of Arts and Sciences, 199–218
.
Finnigan
,
J. J.
, and
F.
Einaudi
,
1981
:
The interaction between an internal gravity wave and the planetary boundary layer. Part II: Effect of the wave on the turbulence structure.
Quart. J. Roy. Meteor. Soc.
,
107
,
807
832
.
Grachev
,
A. A.
,
C. W.
Fairall
,
P. O. G.
Persson
,
E. L.
Andreas
, and
P.
Guest
,
2005
:
Stable boundary-layer scaling regimes: The SHEBA data.
Bound.-Layer Meteor.
,
116
,
201
235
.
Hicks
,
B. B.
,
1978
:
Some limitations of dimensional analysis and power laws.
Bound.-Layer Meteor.
,
14
,
567
569
.
Howard
,
L. N.
,
1961
:
Note on a paper of John W. Miles.
J. Fluid Mech.
,
10
,
509
512
.
Howell
,
J. F.
, and
L.
Mahrt
,
1997
:
Multiresolution flux decomposition.
Bound.-Layer Meteor.
,
83
,
117
137
.
Klipp
,
C. L.
, and
L.
Mahrt
,
2004
:
Flux gradient relationship, self-correlation and intermittency in the stable boundary layer.
Quart. J. Roy. Meteor. Soc.
,
130
,
2087
2103
.
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
1102
.
Louis
,
J. F.
,
1979
:
A parametric model of vertical eddy fluxes in the atmosphere.
Bound.-Layer Meteor.
,
17
,
187
202
.
Mahrt
,
L.
,
1999
:
Stratified atmospheric boundary layers.
Bound.-Layer Meteor.
,
90
,
375
396
.
Mahrt
,
L.
, and
D.
Vickers
,
2005
:
Boundary layer adjustment over small-scale changes of surface heat flux.
Bound.-Layer Meteor.
,
116
,
313
330
.
Mahrt
,
L.
, and
D.
Vickers
,
2006
:
Extremely weak mixing in stable conditions.
Bound.-Layer Meteor.
,
119
,
19
39
.
Mahrt
,
L.
,
J.
Sun
,
W.
Blumen
,
T.
Delany
, and
S.
Oncley
,
1998
:
Nocturnal boundary-layer regimes.
Bound.-Layer Meteor.
,
88
,
255
278
.
Miles
,
J.
,
1961
:
On the stability of heterogeneous shear flows.
J. Fluid Mech.
,
10
,
496
508
.
Miles
,
J.
,
1986
:
Richardson’s criterion for the stability of stratified shear flow.
Phys. Fluids
,
29
,
3470
3471
.
Monin
,
A. S.
, and
A. M.
Obukhov
,
1954
:
Basic laws of turbulent mixing in the ground layer of the atmosphere.
Akad. Nauk SSSR Geofiz. Inst. Tr.
,
151
,
163
187
.
Nappo
,
C. J.
,
2002
:
An Introduction to Atmospheric Gravity Waves.
Academic Press, 260 pp
.
Nastrom
,
G.
,
K.
Gage
, and
W. H.
Jasperson
,
1984
:
Kinetic energy spectrum of large- and mesoscale atmospheric processes.
Nature
,
310
,
36
38
.
Nieuwstadt
,
F. T. M.
,
1984
:
The turbulent structure of the stable, nocturnal boundary layer.
J. Atmos. Sci.
,
41
,
2202
2216
.
Obukhov
,
A. M.
,
1946
:
Turbulence in an atmosphere with inhomogeneous temperature.
Tr. Inst. Teor. Geofiz. Akad. Nauk SSSR
,
1
,
95
115
.
Pahlow
,
M.
,
M. B.
Parlange
, and
F.
Porte-Agel
,
2001
:
On Monin–Obukhov similarity in the stable atmospheric boundary layer.
Bound.-Layer Meteor.
,
99
,
225
248
.
Poulos
,
G. S.
, and
S. P.
Burns
,
2003
:
An evaluation of bulk Ri-based surface layer flux formulas for stable and very stable conditions with intermittent turbulence.
J. Atmos. Sci.
,
60
,
2523
2537
.
Poulos
,
G. S.
, and
Coauthors
,
2002
:
CASES-99: A comprehensive investigation of the stable nocturnal boundary layer.
Bull. Amer. Meteor. Soc.
,
83
,
555
581
.
Richardson
,
L. F.
,
1920
:
The supply of energy from and to atmospheric eddies.
Proc. Roy. Soc. London
,
A97
,
354
373
.
Schumann
,
U.
, and
T.
Gerz
,
1995
:
Turbulent mixing in stably stratified shear flows.
J. Appl. Meteor.
,
34
,
33
48
.
Smedman
,
A-S.
,
X. G.
Larsen
,
U.
Högstro
,
K. K.
Kahma
, and
H.
Pettersson
,
2003
:
Effect of sea state on the momentum exchange over the sea during neutral conditions.
J. Geophys. Res.
,
108
.
3367, doi:10.1029/2002JC001526
.
Sorbjan
,
Z.
,
2006
:
Local structure of turbulence in stably-stratified boundary layers.
J. Atmos. Sci.
,
63
,
1526
1537
.
Tjernström
,
M.
,
1993
:
Turbulence length scales in stably stratified free shear flow analyzed from slant aircraft profiles.
J. Appl. Meteor.
,
32
,
948
963
.
Tjernström
,
M.
,
2005
:
The summer arctic boundary layer during the Arctic Ocean Experiment 2001 (AOE-2001).
Bound.-Layer Meteor.
,
117
,
5
36
.
Uttal
,
T.
, and
Coauthors
,
2002
:
Surface heat budget of the Arctic Ocean.
Bull. Amer. Meteor. Soc.
,
83
,
255
275
.
Van de Wiel
,
B. J. H.
,
A. F.
Moene
,
O. K.
Hartogensis
,
H. A. R.
De Bruin
, and
A. A. M.
Holtslag
,
2003
:
Intermittent turbulence in the stable boundary layer over land. Part III: A classification for observations during CASES-99.
J. Atmos. Sci.
,
60
,
2509
2522
.
Vickers
,
D.
, and
L.
Mahrt
,
2003
:
The cospectral gap and turbulent flux calculations.
J. Atmos. Oceanic Technol.
,
20
,
660
672
.
Vickers
,
D.
, and
L.
Mahrt
,
2006
:
A solution for flux contamination by mesoscale motions with very weak turbulence.
Bound.-Layer Meteor.
,
118
,
431
447
.
Vogelezang
,
D. H. P.
, and
A. A. M.
Holtslag
,
1996
:
Evaluation and model impacts of alternative boundary-layer height formulations.
Bound.-Layer Meteor.
,
81
,
245
269
.
Wilczak
,
J. M.
,
S. P.
Oncley
, and
S. A.
Stage
,
2001
:
Sonic anemometer tilt correction algorithms.
Bound.-Layer Meteor.
,
99
,
127
150
.

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