Non-stationary Generation of Weak Turbulence for Very Stable and Weak-Wind Conditions

Turbulence measurements for very stable conditions near the surface are contrasted among three sites: a high altitude basin during winter with grass or snow-covered grass, a broad valley with complex agricultural land use, and a more narrow valley that is influenced by a valley cold pool and cold air drainage. In contrast to previous studies, this investigation emphasizes the very weak turbulence with large bulk Richardson number occurring during extensive periods between brief mixing events. The relationship of the turbulence to the non-stationary wind and stratification is examined along with the impact of short-term flow accelerations, directional shear and downward diffusion of turbulence from higher levels. The failure of the turbulence for strong stratification to decrease with further increase of stratification is explored. Additional analyses are applied to weak-wind cases for the entire range of stratification, including weak stratification associated with cloudy conditions.

2006; Yagüe et al. 2006;Ohya et al. 2008;Conangla et al. 2008;Vindel and Yagüe 2011). Since the frequency distribution of the turbulence for a given period is normally not bimodal but strongly skewed, the definition of intermittent mixing events can be sensitive to the criteria for defining mixing events (Nakamura and Mahrt 2005). Depending on these criteria, some nights have few or no events (van de Wiel et al. 2003;Steeneveld et al. 2005). The surface turbulent heat flux may become only a minor term in the surface energy balance (van de Wiel et al. 2003).
Turbulence for the most stable conditions does not obey similarity theory (Grachev et al. 2005;Sorbjan 2010), including the lack of an inertial subrange (Grachev et al. 2012a). Traditional turbulence concepts (Tennekes and Lumley 1972) may not apply to the fluctuating flow. For example, such fluctuations can be characterized by very small correlations and vertical velocity fluctuations and much larger horizontal velocity and temperature variations, compared to more traditional turbulence in weakly stable conditions . For strong stratification, the very stable regime can be defined as winds weaker than a threshold wind speed (Sun et al. 2012). Intermittency occurs when the flow switches back and forth across this threshold. For the very weak-wind regime, turbulent eddies may not be interacting with the surface even as low as 2 m above the ground.
Turbulence in the very stable regime might be generated primarily by wave-like motions and other small submeso motions on time scales of minutes or tens of minutes (Conangla et al. 2008) such that equilibrium between the turbulence and non-turbulent flow is not established (Mahrt 2011). In such non-stationary conditions, the surface turbulence may be poorly related to the "mean" flow computed from traditional time averages ranging from 5 min to 1 h. The poor relationship is due to the exclusion of turbulence-generating non-turbulent motions on time scales between the turbulence and the averaging time.
Here, "non-stationary" refers to conditions when the amplitude of these non-turbulent submeso motions is greater than the speed of the large-scale flow. For sufficient stratification, such submeso non-turbulent motions might include the buoyancy subrange (Sukorianski and Galperin 2012). Wind profiles on submeso time scales can become severely distorted compared to those satisfying similarity theory (Mahrt 2008b). Low-level wind maxima occur frequently.
Little emphasis has been placed on the extensive periods of very weak turbulence between mixing events, which is the subject of this study. Since very weak turbulence accounts for most of the record for very stable conditions, it can contribute significantly to the timeaveraged turbulence and to the heat flux . The degree of "weakness" of this turbulence partly determines the rapid cooling rate, build-up of contaminants and the potential for ground fog. Sun et al. (2012) finds that for the very stable regime, the weak turbulence is relatively insensitive to the magnitude of the stratification. This unexpected result also occurs in the three datasets analyzed in this study and is examined in some detail. We also evaluate the variation of turbulence relationships among the three sites.

Observations
We analyze observations from three sites for the nocturnal period between 2100 to 0600 local time. We have excluded the commonly analyzed CASES-99 dataset since the number of very stable cases was significantly less than in the datasets analyzed in this study, listed in Table 1. The total number of nights regardless of stability. Height (m) of the upper-level sonic (u.s.), mid-level sonic (m.s.) and lower-level sonic (l.s.), height of the upper-level temperature (u.T.) and lower-level temperature (l.T.). Parenthetical values refer to the heights of the thermocouples in the "thin-layer" calculations

Sonic Anemometer Measurement Errors
For representation of the turbulence, we analyze sonic anemometer observations at a height of 1 m for all three sites in spite of possible significant loss of fluctuation variance due to path-length averaging. Because the very stable boundary layer is often shallow, with depths <10 m (Smedman 1988;King 1990), the turbulence measured at levels >1 m, where pathlength averaging is less of a problem, is less representative of the surface conditions. The vertical path length of 0.1 m for both the CSAT3 and RM Young 81000 sonic anemometers is of greatest concern because the stratification reduces the vertical length scale of the eddies. Extrapolation of cospectra to smaller unresolved scales indicates only minor loss of the heat flux. However, for very stable conditions, the assumption of constant spectral slopes or spectral similarity theory to estimate such loss may not be valid. In addition, the use of Taylor's hypothesis with very weak winds to convert from the time domain of the measurements to the space domain of the path-length averaging is problematic because the wind direction is constantly changing. Our videos of striated fog indicate that during periods of weakest turbulence, significant fine-scale diffusion would not be captured by the sonic anemometers. Consultation of the manufacturer specifications for the Campbell CSAT and RM Young sonic anemometers used in this study indicate that periods with the weakest turbulence (vertical velocity fluctuations of a few tens of mm s −1 ) are probably contaminated by significant measurement errors. Our wind-tunnel comparisons between the two types of sonic anemometers indicate close comparisons even with such weak winds. Nonetheless, we concede the inability to quantitatively estimate errors for very weak turbulence and that such errors may be important.
Rotation of the data to correct for sonic misalignment was not applied since directiondependent mean vertical motions for strongly stratified flows may be real. Nonetheless, this is an open and important issue for strong stratification. While fog releases reveal that towers and support structures lead to only minor disturbances in very weak flows, an examination of the directional dependence of σ w /V reveals some augmentation of turbulence for flow through the towers prior to reaching the sonic anemometers. However, no observation was excluded since eliminating specific wind directions with meandering wind directions would lead to fragmented time series.
The data were laboriously quality controlled by hand as automated methods can eliminate some of the unusual characteristics found in very stable conditions. Various averaged quantities for each site are shown in Tables 2 and 3 (Sect. 3). Roughness lengths (Table 3) were computed from windy near-neutral conditions. Because similarity theory becomes invalid in very stable conditions, the quantitative application of such numerical estimates of the roughness length to very stable conditions is without formal theoretical support but can be used as a general indicator of the surface roughness for very stable conditions. Vertical gradients of  The columns are the site and level (if different from 1 m), σ w (m s −1 ), σ w (m s −1 ) scaled by the wind speed, V , the correlation between vertical velocity and temperature fluctuations (R), σ w (m s −1 ) adjusted by the correlation coefficient corresponding to Rσ w (m s −1 ), the adjusted σ w scaled by the wind speed, the vertical length scale defined as σ w scaled by the Brunt-Väisälä frequency N (s −1 ) where N is computed from the difference between the top and bottom temperature measurements (Table 1). Products and ratios are computed by first averaging individual quantities and then performing multiplications and divisions Table 3 Averages for very stable conditions (R b > 1) Site wind speed and potential temperature are computed from simple finite differencing between the upper and lower levels of the wind speed and temperature measurements for each site (Table 1).

FLOSS
The first dataset was collected in North Park, Colorado, USA, during the Fluxes Over a Snow Surface II (FLOSSII) experiment within a broad deep valley (Mahrt 2007). Very stable conditions are less frequent at this site (Table 3), but the long observational period of four months leads to a significant number of very stable cases. The surface consists of matted grass, sometimes with a shallow snow cover; the aerodynamic roughness length for this site is quite small, less than 0.001 m with snow cover. This study analyses CSAT sonic anemometer data from the 1, 2, 5 and 10-m levels; the computation of the vertical temperature gradient required recalibration of the data (Mahrt 2007).

BPP
The BPP site (Thomas et al. 2012) is located on the Botany and Plant Pathology farm of Oregon State University, Corvallis, Oregon, USA. The observations analyzed here were collected from late August until mid-October 2011. The site consists of a network of instruments located within a grass-covered shallow depression 20 m across and roughly 50 m long and about 1 m deep, framed by gradual side slopes generally less than 10 %. This observational domain is embedded within a larger flat region of mixed agricultural use that includes vineyards, orchards and isolated buildings. The turbulence is measured with R. M. Young 81000 VRE sonic anemometers. The roughness length is wind-direction dependent but always greater than a few tens of mm. A sonic anemometer in the centre of the domain is selected for analysis in this study, though the conclusions are not sensitive to this choice. The shear is evaluated from the wind speed at the top of a 12-m tower and a 4-station vector average of the wind speed at 1-m height at stations surrounding the tower. The temperature profile at the BPP site is evaluated from 10 aspirated Omega TMTSS-020G thermistors deployed on the 12-m tower.

Rock Springs
The third dataset was collected within the Nittany valley of central Pennsylvania, USA, with sidewalls of 300-350 m elevation gain . We analyze observations from Station 9, located on the valley floor but near the bottom of the sidewall slope. This station was generally located within the valley cold pool on the frequent weak-wind nights. The surface is grass-covered or bare. This station also includes 10 naturally-ventilated Omega TMTSS-020G thermistors for evaluation of the stratification in the lowest 10 m.

Choice of Scales
Based on results below, we nominally define "very stable" as those cases where the bulk Richardson number is >1. Högström (1990) and Basu et al. (2006) have shown that using an averaging window that is too large can lead to false conclusions on the behaviour of the turbulence. Unfortunately, the characteristics of the flow for very stable conditions vary only gradually with scale without a clear cut-off scale that separates turbulence from non-turbulent motions. Thus, any choice of averaging scale is somewhat arbitrary. Fortunately, the statistics of the flow are not sensitive to the exact choice of averaging time. Deviations from 36-s averages, which provide 100 samples per hour, account for most of the systematic vertical flux of heat and is consistent with the analysis of Viana et al. (2010) for very stable conditions. At the same time, deviations from a 36-s average include a substantial contribution from non-turbulent motions leading to low correlations between the vertical velocity fluctuations and other quantities such as temperature. Since turbulence is diffusive by nature and therefore transports scalars, a very low correlation between vertical velocity and temperature fluctuations is viewed as a mix of turbulent and non-turbulent motions. The choice of 36 s can lead to an underestimation of the turbulence covariance for near-neutral windy conditions, but we forego the complexities of a record-dependent averaging length.
This study relates the turbulence to the "mean flow" computed with the same average as is used to compute the fluctuations without the additional averaging of turbulence quantities. Thus, the mean flow includes non-stationarity on scales just larger than the turbulent scales. Our study relates the turbulence, represented by deviations from the 36-s averages, to the 36-s averaged mean flow. Averaging quantities over a longer period used to define the mean flow omits non-turbulent motions between time scales of 36 s and the larger averaging period. Such omitted motions may dominate the shear generation of turbulence for conditions of very weak large-scale flow. As a result, the use of a larger averaging time for the mean flow degrades the relationship between the turbulence and the mean flow.

Interval Averaging
Presenting observations for very stable conditions is problematic because of large scatter, partly due to averaging the turbulence quantities over only 36 s. Consequently, we collect samples of turbulence quantities into subsets based on either different intervals (bins) of the bulk Richardson number or different intervals of the wind speed, and then average turbulence quantities over all of the values within a given interval. A minimum of 10 samples is required for each interval before averaging. Ratios of quantities are computed from interval-averaged quantities. Ratios are not directly averaged since they are often characterized by large skewness and kurtosis and averages can be strongly affected by extreme values. Averaging quantities before taking the ratio in the bulk Richardson number reduces some of the problems with interval (bin) averaging discussed by Grachev et al. (2012b). However, interval averages cannot be interpreted as an estimate of an ensemble average since observations within a given interval do not constitute a uniform population. The flow depends on more than just the bulk Richardson number. The standard error for a given interval is generally small because of the large number of samples. However, because the flow is strongly non-stationary and distributions are strongly skewed for large bulk Richardson numbers, the meaning of the standard error is uncertain.

Stability
The goal here is to predict turbulence exclusively in terms of mean variables without information on the turbulence, as in Mauritsen and Svensson (2007). Therefore, we use a bulk where δθ is the difference of potential temperature between the upper and lower temperature level (Table 1). The wind speed, V , is computed at the upper level, z, while the temperature difference δθ is computed from the potential temperature at the lower level (1 m) instead of using the surface radiation temperature. The footprint and quality of the measurement of the surface radiation temperature vary among sites. Since the wind at the lower level in Eq. 1 is taken as zero at the ground surface, the above form of the bulk Richardson number is asymmetric with respect to treatment of temperature and wind. We have also evaluated the gradient Richardson number based on the shear between the upper and lower levels. Although physically more direct, the shear is more vulnerable to observational errors compared to the wind speed. The shear becomes sensitive to the exact choice of levels with very stable conditions where the wind profiles include a variety of complex structures. Evaluation of the Richardson number over a deeper layer than used here leads to larger values, consistent with Balsley et al. (2008) and others. Thus, numerical comparison of the turbulence-Richardson number relationship among studies is inadvisable. Although, the scatter between the turbulence and the gradient Richardson number is greater than that between the turbulence and the above bulk Richardson number, the turbulence depends on both Richardson numbers in the same qualitative manner. One exception is the most stable cases where directional shear can be systematically important (Sect. 5). In spite of the influence of directional shear, the gradient Richardson number is often an order of magnitude larger than the bulk Richardson number. That is, the magnitude of the velocity difference is generally significantly smaller than the wind speed at the upper level.
A companion ratio is defined as where z is again the height of the upper level where V is evaluated. The ratio V dT can be related to the prediction of the surface heat flux by the bulk relation, which is equal to C h V δθ, where C h is the bulk transfer coefficient for heat. The same value of the bulk Richardson number can correspond to large V and δθ (large V dT ) when winds are significant and stratification is strong, or small V and δθ (small V dT ) when winds are especially weak but the stratification is also weak due to cloud cover. Large V dT can be forced by strong surface radiational cooling or cold-air advection at the surface. Although V dT contains useful information, it lacks direct physical roots such as in the derivation of the Richardson number from the turbulence energy equation.

Turbulence Quantities
The standard deviation of the vertical velocity fluctuations based on deviations from 36-s averages, σ w , is used as a general indicator of the intensity of turbulence since horizontal velocity fluctuations in a stratified flow are more dominated by non-turbulent motions compared to σ w . However, even σ w can be contaminated by non-turbulent submeso motions. Because the turbulent and non-turbulent scales appear to overlap, the turbulence cannot be isolated from non-turbulent motions. As an alternative approach, we define a modified velocity scale where R is the correlation coefficient between vertical velocity and temperature fluctuations.
Here, Rσ w better represents the mixing efficiency and vanishes with vanishing correlations, and it is more related to the transporting part of the vertical motions that are correlated with temperature. Analogous versions of Eq. 3 could be constructed for other scalars. Although Rσ w is expected to reduce the influence of non-turbulent motions, it does not eliminate them. The ratio of the fluctuation potential energy to the vertical component of the kinetic energy (Mahrt 2011) provides an additional characteristic of the fluctuating flow. This ratio is large with strong stratification when a given vertical velocity fluctuation produces larger temperature fluctuations.
A measure of the directional shear is constructed as the difference between the wind vector and speed shear, written as where δ refers to differences between the top and bottom sonic levels (Table 1) and V is the wind speed. This measure of wind-direction shear avoids issues with the cyclic behaviour of the wind direction but includes no information on the sign of the wind-direction change.

Wind Speed
The frequency distribution of σ w at 1 m for 36-s windows is bimodal for FLOSS (Fig. 1), a distribution that is caused by the bimodal distribution of synoptic wind speed. The frequency maximum at 0.04 m s −1 occurs almost completely for very stable conditions (R b > 1), while the broader frequency maximum centered around 0.23 m s −1 corresponds to the weakly stable regime (R b < 1). The weakly stable regime includes both windy and cloudy nights and events embedded within the very stable nights. For the other two sites where the wind speed for most of the nights is low, the frequency maximum for the stronger turbulence is missing, although the distribution is skewed toward larger values. For all three sites, σ w depends approximately linearly on wind speed for a broad range of wind speeds (Fig. 2a). This dependence of σ w on wind speed becomes a little less than linear for stronger winds and a little more complex for very weak winds (Fig. 2b). σ w becomes almost independent of wind speed for the very weakest winds and asymptotes to a non-zero value as the wind speed vanishes. This asymptotic value is approximately 0.035 m s −1 for FLOSS and about 0.02 m s −1 for the BPP and Rock Springs sites. Non-zero σ w for vanishing wind speed is partly due to the inclusion of non-turbulent motions. The importance of this non-turbulent contribution is suggested by the rapid decrease of correlation between vertical velocity and temperature fluctuations, R (Eq. 3), with decreasing low wind speeds at all three sites (Fig. 3). The correlations become near zero for the weakest wind-speed categories. As a result, Rσ w (not shown) decreases faster than σ w with decreasing low wind speed. Non-zero σ w for vanishing wind speed is also related to the transient nature of weak winds and finite decay time scale of the turbulence. Even if the wind vanishes, decaying turbulence continues. Additional influences include: the vertical diffusion of turbulence from adjacent layers and non-zero shear even when the wind speed vanishes at a given level. The role of noise in the data is not known.
For the weakest wind speed, the slope of σ w as a function of V increases to a larger value for wind speed greater than a transition wind speed (Fig. 2), as in Sun et al. (2012). The transition appears to be gradual with the stretched x-axis of Fig. 2b. This transition occurs at  Fig. 4 The dependence of the interval-averaged σ w as a function of the wind speed for FLOSS for 1 m (black), 2 m (red), 5 m (green) and 10 m (green) for weak winds only three sites, the asymptotic value of σ w increases with increasing roughness length while the transition value of the wind speed decreases with increasing roughness length. The transition wind speed at the FLOSS site increases with height ( Fig. 4), as in Sun et al. (2012). In contrast to Sun et al. (2012), this transition becomes more obscure with height, possibly due to the influence of surface heterogeneity at the FLOSS site and increasing footprint with height. This trend continues to the top of the 30-m tower (not shown). For a given wind speed below the transition value, σ w increases with height ( Fig. 4), pointing to a stronger source of turbulence at higher levels, sometimes referred to as the "upside down" boundary layer. For a given wind speed above the transition value, σ w decreases with height, as in a classical atmospheric boundary layer.
The standard deviation of σ w within each wind-speed category is small compared to the mean value for significant wind speeds but increases to about half of the mean value for the weakest wind-speed category. This large scatter is due to non-stationarity and other influences on σ w , discussed above. Due to the large number of samples, the standard error is small compared to the mean value of σ w and increases only to a few percent for the weakest wind categories. However, the concept of the standard error is obscure for strongly non-stationary flow with highly skewed frequency distributions.

Stability
The effects of wind speed and stratification are traditionally combined by using the Richardson number as a stability parameter. For all three sites, σ w decreases rapidly with increasing R b for small R b , as found in previous studies and shown in Fig. 5a (black) for the FLOSS site. For larger bulk Richardson number, σ w decreases much more slowly with increasing bulk Richardson number. σ w becomes almost independent of R b for R b > 10, visually obvious when accounting for the logarithmic scale of the x-axis in Fig. 5a. In fact, for the largest bulk Richardson numbers, σ w increases slightly. If significant, this increase is probably due to increasing non-turbulent contributions to σ w associated with an extremely small correlation between vertical velocity and temperature fluctuations.
Subject to the qualification that the standard error has uncertain meaning in non-stationary conditions with strongly skewed frequency distributions, the standard error is quite small, partly due to the large number of samples. The standard error divided by the interval average of σ w does not exceed 10 % until the interval of largest Richardson number where the sample size is small. The within-interval standard deviation of σ w scaled by V increases steadily with increasing R b and reaches 80 % for the last three intervals of the largest bulk Richardson number.
Relating scaled quantities to each other provides the possibility of universal similarity theory at the risk of introducing significant self-correlation. Since we are attempting to isolate the influence of the non-turbulent flow on the turbulence, the usual choice of friction velocity u * (Pahlow et al. 2001) is not considered here as a scaling variable for σ w . In addition, u * varies erratically for very stable conditions (Acevedo et al. 2009) where the stress, wind and shear vectors may assume significantly different directions. With very large stability, the surface stress may no longer be a significant influence on the turbulence (Grachev et al. 2005). Instead, we scale the turbulence by the speed of the "non-turbulent" flow V defined in terms of the 36-s averages. The resulting ratio σ w /V is roughly independent of R b for the entire range of R b (red, Fig. 5a); the Rock Springs and BPP datasets also show near independence of σ w /V from R b .
Although difficult to estimate, positive self-correlation due to V in the denominator of both σ w /V and R b is expected to be small (Mahrt 2008a). Since σ w is nearly independent of the bulk Richardson number for the very stable regime, significant positive self-correlation would have produced an increase of σ w /V with increasing R b . This is not observed. The near independence of σ w /V from R b implies that most of the variance of σ w is explained by variations of V and stratification is less important. This speculation is explored further in Sect. 5. Similar results are obtained using the bulk Richardson number evaluated over a thinner layer defined in Table 1.

Site Differences
While the qualitative dependence of σ w and σ w /V on R b are similar for all three sites, quantitative differences occur. For the very stable regime, σ w /V at the BPP site is about twice as large compared to the other two sites ( Table 2). The BPP site is rougher and more heterogeneous than the other two sites. At the same time, σ w is greatest for the FLOSS site, because the winds in FLOSS are not as weak compared to the other two sites for the very stable regime. That is, for a given bulk Richardson number, both the wind speed and stratification are greater at the FLOSS site compared to the other two sites, which can be quantified in terms of a larger V dT (Table 3). The very strong stratification at the FLOSS site is partly due to strong radiational cooling associated with partial snow cover and high-altitude (2,400 m) winter conditions, where concentrations of water vapour and particulate matter in the overlying atmosphere are small. This leads to strong radiational cooling. The Rock Springs site is characterized by the smallest values of σ w /V and V dT , due to very weak winds and weaker (but still strong) stratification compared to the other two sites.
As a possible consequence of the larger value of V dT at the FLOSS site, the averaged kinematic heat flux is significantly larger at the FLOSS site: −1.5 × 10 −3 K m s −1 compared to −0.5 × 10 −3 K m s −1 for BPP and −0.7 × 10 −3 K m s −1 for Rock Springs. At the FLOSS site, P E/K E is smaller than at the other two sites (Table 3) because of more significant wind and larger σ w .
Isolation of the influences of V dT and surface roughness is not possible with just three sites. Neither R b nor V dT contain information about the vertical structure of the flow, which, for very stable conditions, varies substantially. For the FLOSS and Rock Spring sites, a nearsurface wind-speed maximum occurs about half of the time for the 36-s averages (Table 3). This percentage is only about 10 % at the BPP site. This percentage is smaller than that for the other two sites, possibly due to greater roughness, to the absence of significant slopes in the region and to the location of the upper sonic level at a height of 12 m, in contrast to 10 m and 8.4 m at the other two sites. The impact of near-surface wind maxima on the turbulence is difficult to isolate.

Rσ w
For very stable conditions, σ w cannot be interpreted purely as turbulence. For example, the correlation between the vertical velocity and temperature fluctuations for very stable conditions decreases in magnitude to near −0.1 (Table 2), similar to that found by Sorbjan (2006), as compared to about −0.25 for weak stability in the present datasets and −0.35 for nearneutral conditions in Sorbjan (2006). Turbulence is normally not associated with correlation coefficients with magnitudes as small as 0.1. Shrinking the averaging window by a factor of 5 or 10 only modestly increases the correlation coefficient, but excludes significant downward turbulent heat flux.
Because both σ w and the correlation coefficient decrease with increasing R b , Rσ w decreases faster than σ w with increasing R b . In this sense, Rσ w better isolates the very stable regime. Rσ w becomes extremely small and relatively constant for R b greater than unity (black, Fig. 5b). From another point of view, significant correlation is confined to increasingly smaller time scales with increasing stability, thus allowing greater influence of non-turbulent motions within a fixed averaging time . In contrast to σ w /V (red, Fig. 5a), Rσ w /V (red, Fig. 5b) is not independent of stability but decreases modestly with increasing bulk Richardson number until R b = O(1) and then becomes relatively independent of the bulk Richardson number. Thus, use of Rσ w /V instead of σ w /V restores some dependence on stability by emphasizing the heat-transporting part of the flow. Similar results are found for the other two sites. The smallest values of R b are off scale in Fig. 5 where the correlation between the vertical velocity fluctuations and the small temperature fluctuations becomes small, causing Rσ w to decrease. This regime is outside the scope of our study and is characterized by large relative errors in the small values of δθ.

Very Stable Conditions; Dependence on δθ
Since σ w and σ w /V both become relatively independent of the bulk Richardson number for very stable conditions (R b > 1), we now investigate the behaviour of σ w for very stable conditions by collecting all of the observations with R b > 1 into one data subset. This condition partitions the flow according to the bulk Richardson number, in contrast to intermittency studies that directly divide the time series into turbulent events and quieter periods between the events. Relating the turbulence to the non-stationary mean flow is the goal of this investigation.
In general, a bulk Richardson number >1 corresponds to quite weak turbulence with only a few cases of more significant turbulence. For example, at the FLOSS site, only 0.2 % of the 36-s values of σ w exceed 0.2 m s −1 . Some of these cases are probably decaying mixing events that have not reached equilibrium with the increasing stability. However, the vast majority of the events with significant turbulence are associated with R b < 1. Care must be taken when examining inter-variable relationships for R b > 1 because the imposed lower limit on the bulk Richardson number corresponds to a maximum allowed wind speed that increases according to the square root of the stratification, δθ. Scaling σ w by the wind speed reduces the impact of this stratification-dependent restriction on wind speed.
The general independence of σ w /V from δθ (Fig. 6b) again supports the inference that the turbulence for very stable conditions depends mainly on the wind speed (Fig. 6a) and much less on stratification. The ratio Rσ w /V also fails to systematically decrease with increasing stratification for very stable conditions (not shown). Evaluating the stratification using different combinations of levels of temperature observations did not alter this conclusion. The reduction of σ w due to vertical path-length averaging should become more important with increasing stratification (reduced eddy size) so that path-length averaging would reduce the measured σ w /V with increasing stratification. This is not observed.
The failure of σ w /V and Rσ w /V to decrease with increasing stratification for R b > 1 may be due to increasing directional shear. We now examine the impact of wind-directional shear, defined as the magnitude of the vector shear minus the speed shear (Eq. 5). σ w increases systematically with wind-directional shear for the Rock Springs and BPP sites (Fig. 7), but not for the FLOSS site. The wind-directional shear increases systematically with increasing stratification for the FLOSS and Rock Springs sites but not for the BPP site. The site dependence of these tendencies might be related to the different causes of directional shear among the three sites. Stronger directional shear at the Rock Springs site is often associated with very weak north-easterly along-valley flow at the surface (1 m) within the valley cold pool and "less cold" south-south-easterly flow at 8.4 m that is roughly aligned with the slope to the southeast of the station. The latter could be drainage flow that overrides the colder air in the valley cold pool . The strong directional shear associated with very strong stratification at the BPP site is often confined to the lowest few metres but occurs in a variety of situations. The increase of directional shear for strong stratification at the FLOSS site is often due to cold northerly flow from the direction of a small cold pool, undermining prevailing southerly winds (Mahrt 2010). Since the cold pool begins only a few tens of metres north-north-west of the tower and turbulence in the cold pool is extremely weak, the flow may pass the tower before significant turbulence is generated by the directional shear. These results suggest that the increase of directional shear with increasing stratification might, at least partially, offset the influence of increasing buoyancy destruction with increasing stratification, but examination of the important site dependence requires data from more sites.
For very stable conditions in the stationary polar night, Ekman turning also contributes to directional shear close to the surface, which increases with increasing stability (Grachev et al. 2005). Although the very stable nocturnal boundary layer is generally non-stationary as well as heterogeneous, Ekman effects could contribute to the directional shear even if not a dominating factor. In addition to the influence of wind speed, stratification, and wind-direction shear, a variety of other influences potentially governs the behaviour of σ w and σ w /V .

Non-stationarity
From a Lagrangian point of view, the flow is constantly adjusting to changing surface conditions. Even minor surface features can become important with strong stability. In addition, the turbulence is always changing due to non-stationary external forcing by wave-like and other submeso motions on time scales that are not large compared to those of the largest turbulent eddies. Internal intermittency due to the interplay between the "mean flow" and the turbulence may also occur (Pardyjak et al. 2002;Fernando 2003), but is difficult to isolate from observations at fixed points. Internal intermittency is viewed as a sequence of events beginning with increasing shear (decreasing Richardson number), enhanced generation of turbulence, reduction of shear by the mixing, subsequent decay of the turbulence, rebuilding of shear, and so forth.
Acceleration was computed as V (i) − V (i − 1) where i represents the current 36-s time interval used to evaluate σ w , and i − 1 is the previous 36-s interval. During periods of flow acceleration (increasing speed between successive intervals), the development of turbulence lags the increasing wind speed because shear instability requires a finite amount of time to generate turbulence. As a result, the turbulence relative to the current value of the wind speed (σ w /V ) is smaller compared to that for vanishing acceleration (Fig. 8). During periods of deceleration, the decaying turbulence does not decay sufficiently rapidly to maintain equilibrium with the weakening wind field due to the finite time scale of the turbulence decay. As a result, σ w /V is larger than its value for vanishing acceleration (Fig. 8). In very stable conditions, vanishing acceleration generally corresponds to a time change in the sign of the acceleration and does not necessarily imply equilibrium between the turbulence and mean flow. The cause of the much greater impact of accelerations for the BPP site compared to the marginal impact for the FLOSS site (Fig. 8) is not known. As an aside, the important influence of the acceleration on σ w /V extends to the weakly stable regime as well.

Vertical Transport
For very stable conditions, σ w , on average, increases with the height at all three sites. σ w /V at all three sites increases substantially with the downward transport of turbulence toward the surface ( w 3 < 0), as shown in Fig. 9. Apparently, downward bursts of turbulence augment σ w more than they increase the wind speed through downward momentum transport. Figure 9 also shows that σ w /V increases with increasing upward transport (w 3 > 0). The upward transport of turbulent energy is associated with periods of greater σ w and σ w /V due to turbulence generated by surface-based disturbances.

Sensitivity and Weak Winds
Increasing the averaging window from the default value of 36 s to several minutes increases σ w by a small fraction and modestly decreases the correlation coefficient. While σ w increases with height for very stable conditions, the relationship of σ w to other variables does not qualitatively change with height, at least in the lowest 10 m. The above results are not sensitive to the exact choice of the specified threshold value of the bulk Richardson number for defining very stable conditions (here, unity). The main impact is through increased scatter with larger values of the specified threshold Richardson number due to fewer samples.
It is not clear what results for R b > 1 require very strong stratification or require weak winds or require both conditions. For example, weak winds generally correspond to strong non-stationarity on submeso time scales, regardless of the stratification. Weak-wind cases are now examined for all values of the stratification, which includes weakly-stratified cases, generally corresponding to cloudy conditions. Here, we arbitrarily define weak winds where the 1-m wind speed <1 m s −1 . Almost the entire BPP nocturnal dataset (87 %) satisfies this weak-wind condition while only about 21 % of the FLOSS observations and 57 % of the Rock Springs observations satisfy this condition. Alternatively, one could attempt to define a sub-class of observations where the speeds are less than the transition speeds (Sect. 4.1), but this relies on our ability to define the transition speeds.
For weak winds, σ w /V decreases with increasing stratification for weak stratification at all three sites (Fig. 10), as also found in Sun et al. (2012). For stronger stratification, σ w /V remains small but relatively independent of stratification at the BPP site, increases by a small amount, of unknown significance, at the FLOSS site, and increases more rapidly with increasing stratification at the Rock Springs site. Rσ w /V does not show an obvious increase with increasing stratification, although site differences remain. In any event, these results confirm that the near independence of σ w /V from the stratification does not extend to the weakly stratified regime.

Conclusions
Relationships between variables for the very stable boundary layer become more complex compared to the weakly stable case. The turbulence is simultaneously generated by different non-stationary mechanisms. In addition, site characteristics play a much larger role with stronger stability because even slight slopes and depressions can lead to drainage flows and cold pools. Nearby topography can induce propagating modes that are responsible for some of the turbulence generation.
The intensity of the turbulence, here represented by σ w for 36-s windows, decreases rapidly with increasing bulk Richardson number until the bulk Richardson number reaches a transition value of about unity and then decreases very slowly with further increases of the bulk Richardson number. Although the loss of turbulence variance to vertical path-length averaging is thought to be important in the very stable regime, this loss does not create a detectable decrease of σ w with increasing bulk Richardson number.
In contrast to previous studies, we focused on the variation of the weak turbulence between the intermittent mixing events by selecting periods where the bulk Richardson number >1. In addition to the Richardson number, the turbulence also depends on a companion non-dimensional number proportional to the product of the stratification and the wind speed (V dT , Eq. 2). That is, the same value of the bulk Richardson number can correspond to stronger stratification and winds, or weaker winds and stratification. The turbulence and correlation between w and θ are stronger (but still weak) for larger V dT .
Within the very stable regime, the intensity of the weak turbulence is proportional to the constantly changing submeso wind speed, although with significant scatter. While the scatter is partly related to the short averaging time, the turbulence is also influenced by other mechanisms. Downward diffusion of turbulence, w 3 < 0, augments σ w /V . This diffusion is less dramatic than the downward bursting that occurs with reduced stability, but is an important mechanism for mixing at the surface in very stable conditions. Because of the finite adjustment time of the turbulence, σ w /V is smaller for accelerating flow and larger for decelerating flow.
These investigations were also carried out with Rσ w /V where R is the correlation between the vertical velocity and temperature fluctuations. This combined variable attempts to partially filter out those non-turbulent motions that contribute to σ w but do not contribute to the turbulent heat transport. The correlation-based variables, Rσ w and Rσ w /V , more sharply separate the weakly and strongly stable regimes. However, the relationship of Rσ w to other flow characteristics showed behaviour that was qualitatively similar to σ w .
For very stable conditions, σ w /V is roughly independent of the temperature stratification. That is, for a given wind speed, the impact of suppression of the turbulence by increasing stratification is not observed, as also found in Sun et al. (2012). The reason for this unexpected near independence from the stratification is complex, but is partly associated with development of strong directional shear with strong stratification. Strong directional shear often occurs with cold air drainage over a thin cold pool at the Rock Springs site, transient advection of cold air from a nearby shallow cold pool at the FLOSS site, and a variety of transient situations at the BPP site. These mechanisms require further investigation, preferably with a larger sample size and more sites, before forming definite conclusions on the relationship between the turbulence and mean flow for very stable conditions. The influence of sonic misalignment and coordinate rotation, as well as vertical path-length averaging, need to be examined in more detail.