Estimation of upper bounds for the applicability of non-linear similarity functions for non-dimensional wind and temperature profiles in the surface layer in very stable conditions

(a) Theoretical background

Using MOS theory, the wind speed (U) and potential temperature (θ) profiles in the horizontally homogeneous and stationary surface layer are given by

2.1

2.2

2.3

2.4

2.5

Another stability parameter Ri_B is defined as

2.6

2.7

The parameter ζ is not known a priori from the observations, and its computation involves an iterative procedure to calculate u_* and θ_* from equations (2.1–2.3). Over the years, efforts have been made to compute ζ, with an indirect approach in terms of Ri_B, which can be readily obtained from routine meteorological measurements. A systematic mathematical analysis was carried out by Sharan et al. (2003a) to find the upper limit of Ri_B for the applicability of the linear functional forms of Φ_M and Φ_H proposed by Businger et al. (1971) and Dyer (1974). However, such an analysis is not possible for the modified functional forms of Φ_M and Φ_H as these lead to non-linear forms of integrated similarity functions, and consequently, equation (2.7) becomes a transcendental equation in variable ζ. To overcome this, an alternative approach based on the behaviour of C_D and C_H (Carson & Richards 1978; Sharan & Kumar 2008) is used here to examine the upper bound of applicability of the non-linear similarity functions.

The coefficients C_D and C_H are defined as

2.8

2.9

2.10

2.11

(b) Upper bound estimation

In stable stratification, C_D and C_H are observed to decrease continuously with increasing values of ζ (Stull 1988; Aditi 2007). Thus, the increasing behaviour of either C_D or C_H with ζ is not physically realistic. Carson & Richards (1978) exploited this fact to infer graphically the critical value of Ri_B up to which the functional forms of Φ_M and Φ_H are applicable for ζ≥1. They assumed Φ_M=Φ_H and pointed out that C_D approaches zero as ζ tends to infinity. However, this approach has a number of limitations; C_D or C_H may tend to a minimum at intermediate values of ζ. In their study, consideration of the functional forms of Φ_M and Φ_H valid only for ζ≥1 in the whole interval (ζ₀, ζ) for the computation of corresponding ψ_M and ψ_H is not correct. In view of this, a systematic analysis is carried out to estimate the lowest values of ζ at which C_D and C_H achieve their minimum values, within the framework of similarity theory under stable conditions. Further, the minimum of these lowest values of ζ for C_D and C_H will yield the upper bound of Ri_B in equation (2.7).

Note that C_D is a positive continuous function of ζ, and its nature completely depends on the behaviour of Y²_M in the denominator of equation (2.8). The lowest value of ζ, where C_D achieves its first minimum value, corresponds to the value of ζ for which the function Y²_M assumes a maximum value. Thus, the problem of minimization of C_D is equivalent to the maximization of Y²_M. Accordingly, we are looking for the lowest value of ζ, which gives maximum Y²_M. In equation (2.10), z₀ and z are parameters, L is considered a variable and ζ₀=z₀/zζ is a function of ζ and thus, Y²_M is actually a function of single variable ζ. The critical points of Y²_M are obtained from

2.12

Here, Y_M cannot be zero as C_D approaches infinity and thus this case is not of physical interest. Accordingly, the critical points of Y²_M are given by

2.13

If equation (2.13) does not have any root, Y²_M is either a continuously decreasing or an increasing function of ζ. For continuously decreasing Y²_M or increasing C_D, the functional form of Φ_M will not be applicable under stable conditions. For continuously increasing Y²_M or decreasing C_D for ζ>0, the functional form of Φ_M can be used for all values of ζ up to infinity, implying that C_D approaches minimum as .

Similarly, the lowest value of ζ at which C_H has its first minimum can be estimated by analysing the behaviour of denominator Y_MY_H in equation (2.9). The critical points of Y_MY_H are given by

2.14

The upper bound of Ri_B is estimated at the minimum of these upper bounds minimizing C_D and C_H. If both C_D and C_H are continuously decreasing functions of ζ approaching minima as , the corresponding upper bound of Ri_B can be inferred by taking the limit as in equation (2.7).

(a) Similarity functions with Φ_M=Φ_H

In the regime of continuous turbulence when Φ_M can be assumed to be same as Φ_H, C_H transforms to C_D as long as the roughness lengths of momentum (z_0m) and heat (z_0h) are equal to z₀. Thus, in this case, it is sufficient to analyse the behaviour of C_D to infer the lowest value of ζ at which C_D achieves its first minimum value.

(i) Similarity function of Webb (1970)

Similarity function proposed by Webb (1970) (hereafter WE-70) is given as

3.1

where β is a positive universal constant estimated from the observations. This is essentially an extension of the linear relationship proposed by Dyer (1974) by a constant function for ζ≥1 such that the continuity is maintained at ζ=1.

The integrated similarity functions ψ_M and ψ_H, corresponding to Φ_M and Φ_H (equation (3.1)), in equations (2.4) and (2.5) are given as:

3.2

For Ψ_M (equation (3.2)) corresponding to ζ≥1, equation (2.13) can be written as

3.3

Note that from equation (3.3), the critical point occurs only at ζ=z/z₀. At ζ=z/z₀, the second derivative Y_M′′=−β/ζ²<0, implying that a local maxima of Y_M occurs at ζ=z/z₀. Y_M′ in equation (3.3) is positive for 1<ζ<z/z₀ and becomes zero at ζ=z/z₀, and it should be noted that at ζ=1, is also positive. This implies that Y_M is a continuously increasing function of ζ and has all positive values in this range of ζ from 1 to z/z₀. Positively increasing behaviour of Y_M in this interval ensures that Y²_M is also an increasing function and accordingly C_D is a decreasing function of ζ in the range from 1 to z/z₀ (figure 1a). This decreasing behaviour of C_D for all values of ζ in the interval (0,z/z₀) is consistent with the physical behaviour of C_D with increasing stability in the stable atmosphere.

• Download figure
• Open in new tab
• Download PowerPoint

For ζ>z/z₀, Y_M′ in equation (3.3) is always negative and thus the negative sign of Y_M′ implies that Y_M is a decreasing function of ζ. Since Y_M is positive at ζ=z/z₀, Y_M decreases as ζ increases and vanishes at ζ=ζ* (say) beyond which Y_M changes its sign and achieves all negative values. Therefore, Y²_M decreases in the interval (z/z₀,ζ*) and approaches zero as ζ → ζ* and then starts increasing continuously for ζ>ζ*. This implies that C_D increases in the interval (z/z₀,ζ*), attaining its maximum value at ζ → ζ* beyond which it decreases as ζ increases further (figure 1a). These increasing and decreasing behaviours of C_D for ζ>z/z₀ are not physically realistic.

For Webb’s function, it is noted that C_D also tends to zero when Y_M approaches infinity as . However, this point is larger than ζ=z/z₀ (where C_D has its first minimum) and is not of physical interest.

After achieving a minimum value at ζ=z/z₀, the increasing behaviour of C_D is not expected in stable stratification, and Webb’s similarity function is therefore not applicable for ζ>z/z₀. Thus, the Webb’s function is applicable as long as ζ≤z/z₀.

The upper value of Ri_B (i.e. Ri_BU) obtained from equation (2.7) at ζ=z/z₀ is given by

3.4

Thus, the similarity function (3.1) is applicable as long as . The upper bound of Ri_B for the applicability of Webb’s function is finite and increases as the ratio z/z₀ increases (Sharan & Kumar 2008).

The Webb function is purely empirical in nature. In general, its applicability and limitations are based on experimental data. A limiting value of ζ about 10 is indicated by observations (Webb 1970). Micrometeorological measurements indicated that the idealized surface layer of the similarity theory may not begin even at the top of the roughness elements of average height h₀, which is an order of magnitude larger than z₀, but at 1.5–2 times h₀. Thus, z/z₀ should be larger than 10 in order to be outside the roughness layer (Arya 1988). From the point of view of observations, the condition ζ≤z/z₀ is expected to be satisfied for all the ratios of z/z₀.

This analysis also shows that the function is not applicable for ζ>z/z₀. This implies that z₀>L. In this regime, stability effects are important in the scale of the surface roughness elements and the MOS theory becomes inapplicable.

(b) Similarity functions with Φ_M≠Φ_H

In the regime of continuous turbulence in stable conditions, experimental data supported the assumption of Φ_M=Φ_H in the formulations of WE-70, CL-70, H-76 and HD-88. However, in the regime of intermittent turbulence, the exchange of heat is far less efficient than the exchange of momentum, i.e. Φ_H>Φ_M (Beljaars & Holtslag 1991). Consequently, different profiles for the similarity functions of heat and momentum have been proposed by Brutsaert (1982), Beljaars & Holtslag (1991), Cheng & Brutsaert (2005) and Grachev et al. (2007) as follows.

(i) Businger’s extended similarity function (Brutsaert 1982)

 Brutsaert (1982) (hereafter B-82) derived the profiles of Φ_M and Φ_H from an extension of linear relationships proposed by Businger et al. (1971) by constant functions for ζ≥1, such that the continuity is maintained at ζ=1. Thus, in B-82, the formulation of Φ_M is the same as in WE-70 (equation (3.1)), and the formulation of Φ_H is

3.5

where γ is a constant and Pr_t is the turbulent Prandtl number. The values of β and γ are related by β=γPr_t. Corresponding to this Φ_H (equation (3.5)), Ψ_H is given by

3.6

Here, Φ_M is the same as in WE-70 (equation (3.1)), and thus it is applicable as long as ζ≤z/z₀ (§3a(i)).

To estimate the lowest value of ζ for which C_H gets its first minimum value, equation (2.14) can be written as

3.7

Equation (3.7) has its lowest critical point at ζ=z/z₀. For 1<ζ<z/z₀, the derivative Y_MY_H′+Y_M′Y_H of the denominator Y_MY_H in C_H (equation (2.9)) given by equation (3.7) is positive and thus Y_MY_H is an increasing function of ζ. For ζ>z/z₀, Y_MY_H′+Y_M′Y_H is negative up to a point ζ=ζ* (say) and thus Y_MY_H will have the decreasing behaviour. Consequently, C_H decreases in 1<ζ<z/z₀ and gets its first lowest value at ζ=z/z₀ and then starts increasing for ζ>z/z₀ (figure 1b). The decreasing behaviour of C_H in the range of ζ from 1 to z/z₀ is expected in the stable stratification with increasing stability.

In the analysis of the B-82 function, both C_D and C_H decrease in the range of ζ from 1 to z/z₀ and get their minimum values at ζ=z/z₀. The minimum of these two upper bounds of ζ for C_D and C_H will remain at the same point at ζ=z/z₀. The upper bound of Ri_B at this point ζ=z/z₀ is given as

3.8

Thus, B-82 is applicable as long as ζ≤z/z₀ and .

(a) CASES-99 dataset

CASES-99 was an intensive field experiment that was performed during the month of October 1999 near Leon (37^°38′ N; 96^°14′ W) in southeast Kansas in USA (Poulos et al. 2002). The terrain of the main experimental site was relatively flat and homogeneous without any obstacles in the surroundings.

The measurements were taken on a 60 m tower, where the turbulence data were measured at 1.5 (0.5), 5, 10, 20, 30, 40, 50 and 55 m levels by eight sonic anemometers. The turbulence data provided three components of wind (u, v, w) and virtual temperature at a sampling rate of 20 Hz, which were used in the analysis to calculate the fluxes. After the analysis of turbulence data at 10 m level, 254 hourly averaged data corresponding to stable conditions (i.e. L>0) are selected. Out of these, 66 (approx. 26%) hourly data are in stable conditions characterized by z/L>1, and these are used for evaluation purpose in the present study. In general, the fluxes are computed from the hourly observed series from sonic anemometer using the eddy correlation method. Under very stable conditions, the hourly observed turbulent fluxes based on the classical eddy correlation method might have been influenced by the mesoscale part of the flow, which is highly variable (Mahrt 2008). Thus, true observed turbulent fluxes at the 10 m level are determined from sonic observations by using the methodology based on the multi-resolution co-spectral gap proposed by Vickers & Mahrt (2003, 2006). Hourly fluxes are calculated by dividing the 1 h period into a number of subintervals in accordance with the co-spectral gap. In each subinterval, the fluxes are calculated and then their average is taken for the flux corresponding to 1 h.

To compute the surface layer fluxes by MOS theory, hourly averaged temperatures at the surface and 10 m level and wind speed at this level are used. The surface temperature is estimated using the mean temperature measured by five downward-looking narrow-beam Everest infrared radiometers, which were within 300 m of the 60 m tower (Poulos & Burns 2003). The surface roughness length is taken as 0.03 m. With these values, Ri_B is computed from equation (2.6). This computed value of Ri_B with the integrated similarity function in equation (2.7) is used to compute the lowest positive root ζ of equation (2.7) numerically using an iterative method. This value of ζ is used to compute the values of u_* and θ_* from equations (2.1) and (2.2), respectively.

In the analysis of observations in very stable conditions for z/L>1, it is observed that approximately 67 per cent points correspond to Ri_B greater than 0.2. Table 3 shows that out of 66 data hours observed in very stable conditions, approximately 6, 12 and 12 per cent data hours do not satisfy the criteria for applicability of similarity functions of CL-70, H-76 and HD-88, respectively, and thus these functions are unable to compute the surface fluxes for the respective data hours. It is found that the remaining formulations for Φ_M and Φ_H proposed by WE-70, B-82, BH-91, CB-05 and GR-07 fall in the category that compute the surface fluxes for all selected data hours of CASES-99 observations in very stable conditions.

Table 3.Theoretical values of upper bounds of ζ (i.e. ζ_U) and Ri_B (i.e. Ri_BU) for respective similarity functions for the CASES-99 dataset. The last column represents the number of data hours that fall short of satisfying the criteria for applicability of respective similarity function in very stable conditions. All the remaining functions can be used for all Ri_B and ζ as their upper bounds are infinitely large.

Collapse

CASES-99 dataset (z=10 m, z0=0.03 m)
similarity functions	ζU	RiBU	hours (RiB>RiBU)
WE-70	333.33	9.53	—
CL-70	50.92	2.26	4 (approx. 6%)
H-76		1.34	8 (approx. 12%)
HD-88		1.32	8 (approx. 12%)
B-82	333.33	9.58	—

The scatter diagram (figure 2) between computed and observed fluxes for similarity functions of CL-70, H-76 and B-82 shows that most of the points lie along the one-to-one line. Fluxes were also computed for remaining functions, and the results were found to be similar for all formulations. The deviations from the one-to-one line appear to be similar for all Φ_M and Φ_H. However, here our objective is not to analyse the relative performance of the similarity functions, but to estimate the upper bounds for their applicability in very stable conditions.

• Download figure
• Open in new tab
• Download PowerPoint

(b) Cardington dataset

To provide the high-quality data for research on ABL meteorology and surface layer processes, a permanently operating site with surface, subsurface and mast-mounted instrumentation is maintained at Cardington in Bedfordshire, UK (52^°06′ N, 00^°25′ W and 29 m above mean sea level). The site is a large grassy field and is reasonably flat. The dataset contains fast response wind and temperature observations generated from ultrasonic anemometers at 10, 25 and 50 m above the ground surface and slow response temperature measurements at 1.2, 10, 25 and 50 m obtained from platinum resistance thermometers. It also contains the humidity, radiation and subsoil measurements. The dataset has recorded measurements averaged over 1, 10 and 30 min intervals. The turbulence quantities and data from slow response sensors are included in the 10 and 30 min datasets only.

Recently, Luhar et al. (2009) carried out a systematic study on the analysis of turbulence data during months of August, September and November 2005 from the Cardington site in very stable conditions. We have used these 3 months data in this study. For this period, 30 min averaged data at the 10 m level are used to determine the hourly averaged mean wind, temperature and turbulent quantities. The surface temperature θ₀ is taken as the measured mean grass infrared temperature. The roughness length z₀ is taken as 0.05, 0.025 and 0.03 m for August, September and November, respectively (Luhar et al. 2009). In computations of observed L (equation (2.3)) and Ri_B (equation (2.6)) in the bulk layer from z₀ to 10 m, the mean virtual potential temperature () is estimated by averaging the mean temperature observed at 10 and 1.2 m and surface level. After analysis of the data at 10 m level corresponding to stable conditions (i.e. L>0), 972 data hours are selected, from which 228 (approx. 24%) data hours in very stable conditions (i.e. z/L>1) are used in the present study.

Out of these 228 data hours, it is observed that approximately 95 per cent points correspond to Ri_B greater than 0.2. Table 4 shows the number of data hours that fall short of satisfying the upper bounds for applicability of the similarity functions of CL-70, H-76 and HD-88 for the Cardington dataset in each month separately. It shows that out of 51 data hours in August 2005 in very stable conditions, approximately 16, 29 and 29 per cent data hours do not satisfy the criteria for applicability of the similarity functions of CL-70, H-76 and HD-88, respectively. Similarly, out of the 88 data hours in September 2005, approximately 17, 32 and 32 per cent data hours and out of the 89 data hours in November 2005, approximately 6, 23 and 24 per cent data hours, surface fluxes cannot be computed by the similarity functions of CL-70, H-76 and HD-88, respectively. For the 3-month dataset, out of the 228 data hours in very stable conditions, approximately 12, 27 and 28 per cent data hours fail to satisfy the criteria for the applicability of similarity functions of CL-70, H-76 and HD-88, respectively, and thus these are incapable of computing the surface fluxes for the respective hours. The remaining formulations for Φ_M and Φ_H allow computations of the surface fluxes for all selected data hours of Cardington dataset in very stable conditions.

Table 4.The same as table 3 for the 3 months of Cardington dataset.

Collapse

similarity functions	ζU	RiBU	hours (RiB>RiBU)
	August 2005 (z0=0.05 m)
WE-70	200.0	6.26	—
CL-70	41.40	1.91	8 (approx. 16%)
H-76		1.35	15 (approx. 29%)
HD-88		1.32	15 (approx. 29%)
B-82	200.0	6.29	—
	September 2005 (z0=0.025 m)
WE-70	400.0	11.10	—
CL-70	54.45	2.38	15 (approx. 17%)
H-76		1.33	28 (approx. 32%)
HD-88		1.32	28 (approx. 32%)
B-82	400.0	11.15	—
	November 2005 (z0=0.03 m)
WE-70	333.33	9.53	—
CL-70	50.92	2.26	5 (approx. 6%)
H-76		1.34	20 (approx. 23%)
HD-88		1.32	21 (approx. 24%)
B-82	333.33	9.58	—

The observed and computed fluxes for Φ_M and Φ_H of CL-70, H-76 and B-82 are shown in figure 3, in which the computed fluxes are scattered around the one-to-one line and the deviations from it appear almost similar for all Φ_M and Φ_H. Fluxes are also computed for remaining functions, and the results are found to be similar for all formulations. The true observed fluxes based on multi-resolution co-spectral gap are not calculated here as the high frequency data were not available and only 30 min averaged observed turbulent fluxes were available on the website.

• Download figure
• Open in new tab
• Download PowerPoint

(a) Validity of MOS theory in very stable conditions

MOS theory may not be valid in very stable conditions (Holtslag & Nieuwstadt 1986; Mahrt 1998). Several attempts have been made to extend the MOS theory by defining the flux and gradient-based local (z-less) scaling in these conditions (Nieuwstadt 1984; Sorbjan 2006). Our theoretical analysis does not recognize the limitations of the MOS theory/scaling and the superiority of the local similarity of z-less stratification.

Recently, Mahrt (2008) has pointed out that the existing similarity theory fails under very stable conditions, because the mesoscale part of the flow, which is highly variable, is required to be considered for predicting the fluxes. While applying the similarity relationship in these conditions, it is essential to remove all non-turbulent contributions to the fluxes and while balancing the surface energy budgets, heat fluxes might be included at larger time scales, regardless of their origins.

(b) Self-correlation

Self-correlation is referred to in the literature as spurious correlation or the shared variable problem that arises when one (dimensionless) group of variables is plotted against another, and the two groups under consideration have one or more common variables (Klipp & Mahrt 2004; Baas et al. 2006). Mahrt (2008) commented that the correlation between C_D and ζ is strongly influenced by self-correlation through the shared variable u_*, and its physical significance cannot be evaluated from the data. In contrast, his analysis showed that the self-correlation between C_H and ζ with shared variables u_* and is of opposite sign to the physical correlation. Thus, the proposed analysis based on the behaviour of C_D and C_H with ζ also suffers from the problem of self-correlation.

(c) Weak wind stable conditions

The weak wind conditions are classified on the basis of surface wind speed (U_s) and the geostropic wind speed (U_g) (Sharan et al. 2003b). The wind is considered as weak when (i) U_s≤2 m s⁻¹ and (ii) U_g≤4 m s⁻¹; otherwise it is considered as strong. During clear nights with weak winds, the land surface cools due to strong long wave radiative cooling and the overlying boundary layer may become very stable. Large values of Ri_B are observed in the observational studies (Sharan et al. 2003b; Aditi & Sharan 2007; Luhar et al. 2009), especially in weak wind stable conditions. For the Cardington dataset, out of the 228 data hours in very stable conditions, approximately 85 per cent hours correspond to weak winds. From these cases of weak wind stable conditions, approximately 98 per cent cases are associated with Ri_B≥0.2. Similarly, in selected CASES-99 data hours, approximately 23 per cent cases correspond to weak winds and in all of these cases, Ri_B is found to be larger than 0.2.

The proposed theoretical analysis is valid in the weak wind stable conditions as long as the underlying assumptions of horizontally homogeneous and stationary flow in MOS theory are satisfied. However, weak wind stably stratified conditions are often non-stationary with very weak turbulence, influenced by the mesoscale motions (Mahrt 2007). To account for the effect of the non-stationary mesoscale motions on the generation of turbulence in very stable conditions with weak winds, Mahrt (2008) generalized the bulk formulation for surface fluxes that leads to a systematic behaviour of drag coefficient. Thus, it is desirable to carry out an analysis, similar to the one proposed here, based on the generalized bulk formulations to estimate the upper bounds for the applicability of the non-linear forms of Φ_M and Φ_H under very stable non-stationary conditions with weak winds.

(d) Measurements uncertainties in very stable conditions

As the fluxes tend to be smaller at night in very stable conditions, it often becomes difficult to have reliable turbulence measurements in these conditions. As a result, several aspects of stably stratified turbulence have remained poorly understood and even controversial until now, and they continue to be the subject of intense research (see Cheng & Brutsaert 2005 and references therein). These cases of very strong stability and weak winds have often been neglected because of the variety of observational difficulties for very weak turbulence (Mahrt 2008). On the one hand, intermittent turbulence in very stable conditions produces strongly non-stationary events during which the validity of turbulent transport and storage measurements is uncertain. Fluxes calculated using a conventional averaging time of several minutes or longer are severely contaminated by poorly sampled mesoscale motions in very weak turbulence (Vickers & Mahrt 2006). One needs to calculate the true observed turbulent fluxes by removing the non-turbulent part of the flow from the turbulent measurements.

(e) Different roughness lengths of momentum (z_0m) and heat (z_0h)

In the present analysis, both roughness lengths of momentum (z_0m) and heat z_0h) are assumed to be equal to the surface roughness length z₀. In general, z₀ is a function of the geometry and the size of the surface elements; however, observational and theoretical studies reveal that z_0m and z_0h differ from each other and can be easily determined from profile measurements or from turbulence data (Blumel 1999). In the case of linear functional forms of Φ_M and Φ_H when the roughness length for momentum differs from that for the heat, an analysis based on the bulk Richardson number Ri_B (Sharan et al. 2003a) is not feasible to estimate the upper bounds for their applicability. However, the analysis presented here is general and can be used to find the upper bounds of these functions with different roughness lengths in very stable conditions.

Here, it should be noted that equations (2.1–2.3) are based on the original MOS hypothesis and are valid when buoyancy effects of water vapour can be ignored. When substantial evaporation occurs and it affects the density stratification in the surface layer, modified MOS relations incorporating the buoyancy effect of water vapour would be more appropriate (Arya 1988). This is not considered in the present analysis.