## Abstract

The computation of surface fluxes by Monin–Obukhov similarity theory with different linear/non-linear similarity functions for non-dimensional wind and temperature profiles becomes limited to specific ranges of *ζ*=*z*/*L* (where *z* is the height above ground and *L* is the Obukhov length) and bulk Richardson number (*Ri*_{B}) under very stable conditions. A systematic mathematical analysis is carried out to estimate the upper bounds of *ζ* and *Ri*_{B} for the extent of applicability of different non-linear similarity functions in the surface layer under these conditions. A generalized methodology is proposed on the basis of momentum drag coefficient (*C*_{D}) and heat exchange coefficient (*C*_{H}) and applied to various non-linear similarity functions available in the literature. A theoretically derived criterion for the applicability of each of the non-linear similarity function is evaluated with observations from Cooperative Atmosphere-Surface Exchange Study-99 and UK Meteorological Office’s Cardington datasets. The evaluation with both datasets for each non-linear similarity function confirms the validity of proposed theoretical results under very stable conditions.

## 1. Introduction

Accurate description of surface fluxes under very stable conditions is required in atmospheric flow and dispersion models as these conditions develop almost every night over land surfaces and assume special significance due to their high pollution potentials. These fluxes determine the mean wind and temperature profiles in the atmospheric boundary layer (ABL).

Monin–Obukhov similarity (MOS) theory is frequently used to compute the surface fluxes under moderately stable conditions, in which the turbulence is continuous and stationary. During clear nights with weak winds, the land surface cools rapidly due to strong long wave radiative cooling, and the overlying boundary layer may become very stable. This very stable boundary layer eludes modelling attempts due to the breakdown of existing formulations of turbulence and due to features found in the atmosphere, which are not normally included in the models (Mahrt 1998). In very stable conditions, the stratified turbulence is weak, sporadic, patchy and no longer continuous and stationary, and hence MOS theory may not be applicable (Nieuwstadt 1984; Sorbjan 2006). In these conditions, the turbulence no longer communicates significantly with the surface, and there is a strong evidence of the so-called *z*-less stratification (where *z* is a vertical height above the surface) in which mean gradients of turbulence become independent of the height *z* and are scaled by the local fluxes (Nieuwstadt 1984). However, based on an analysis of standard deviations covering almost five orders of magnitude in stability parameter *ζ*=*z*/*L* (where *L* is the Obukhov length), Pahlow *et al*. (2001) found that the *z*-less stratification generally does not hold. Recently, Sorbjan (2006) also observed, from an analysis of Cooperative Atmosphere-Surface Exchange Study (CASES)-99 observations, that the ‘flux-based’ local scaling in the stably stratified boundary layer is valid only in cases with strong, continuous turbulence, when the gradient Richardson number is constant and subcritical.

With regard to the characteristics of the turbulent exchange, it is now well recognized that organized turbulent exchange exists even under very stable conditions and that it relates relatively well to larger scale quantities such as vertical gradients. The main reason why most previous studies failed to recognize this fact is that the temporal and spatial scales of such fluxes become largely reduced under very stable stratification, reaching values as low as 5 or 10 s for the temporal scale (Mahrt & Vickers 2006). Beyond that scale, the exchange is performed by mesoscale fluxes, less organized and not directly associated with vertical gradients. Different roles of turbulent and mesoscale exchange are well characterized by Vickers & Mahrt (2003). While applying the similarity relationships in very stable conditions, it is essential to remove all non-turbulent contributions to the fluxes; while balancing the surface energy budgets, heat fluxes might be included at larger time scales, regardless of their origins. Thus, the correct scale of the turbulent exchange must be taken into account, even more so under very stable conditions, when such scale can be severely reduced. However, finding the correct scale that characterizes the turbulent exchange is a problem in employing the similarity relationships. The existing similarity theory fails under very stable conditions, and the reason is that the mesoscale part of the flow, which is highly variable, should be considered for predicting the fluxes (Mahrt 2008). Whether or not the MOS theory is valid under very stable conditions is still an open question (Cheng & Brutsaert 2005). Recently, Cheng & Brutsaert (2005) and Grachev *et al.* (2007) concluded, from an analysis of CASES-99 and Surface Heat Budget of the Arctic Ocean Experiment (SHEBA) observations, respectively, that MOS theory is valid for the wind profile in both weakly stable and very stable conditions. Despite issues involved with MOS theory, it is still widely used in the computation of surface fluxes in stable conditions.

MOS theory is based on the specification of universal similarity functions of momentum (*Φ*_{M}) and heat (*Φ*_{H}) for non-dimensional wind and temperature profiles in terms of the stability parameter *ζ*. The linear functional forms of *Φ*_{M} and *Φ*_{H} in terms of *ζ* are commonly used to compute surface fluxes in stable conditions (Businger *et al.* 1971; Dyer 1974). The applicability of these linear formulations is limited (Sharan *et al.* 2003*a*) to the bulk Richardson number (*Ri*_{B}), smaller than its critical value *Ri*_{BC}. In moderately stable conditions, *Ri*_{BC} is found to be approximately equal to 0.2. However, large values of *Ri*_{B} are observed in field studies in very stable conditions (Poulos & Burns 2003; Aditi & Sharan 2007; Luhar *et al.* 2009), especially in weak wind stable conditions.

To overcome this, the linear functional forms of *Φ*_{M} and *Φ*_{H} have been modified over the years by various investigators (Clarke 1970; Webb 1970; Hicks 1976; Holtslag & De Bruin 1988; Beljaars & Holtslag 1991; Cheng & Brutsaert 2005; Grachev *et al.* 2007). A comparison of the prediction of several sets of *Φ*_{M} and *Φ*_{H} functions in the limit of very stable stratification has been given by Andreas (2002). Sharan *et al.* (2003*a*) have proposed a systematic mathematical analysis for finding the extent of applicability of linear functional forms. However, such an analysis is not feasible for these modified similarity functions as each of them leads to a non-linear form of integrated similarity function. Blumel (2000) used a function similar to Webb’s (1970) formulation to find the upper value of *Ri*_{B} for the linear part. However, an approximate analytical approach was proposed by exploiting the fact that *Ri*_{B} is known at *ζ*=1 in order to compute *ζ* and the fluxes for *Ri*_{B} larger than 0.2. The extent of the applicability of the similarity function has not been analysed. By assuming *Φ*_{M}=*Φ*_{H}, an alternative approach based on the fact that the drag coefficient (*C*_{D}) approaches zero as to infer the critical value of *Ri*_{B} was proposed by Carson & Richards (1978). However, this approach has a number of limitations. *C*_{D} may tend to a minimum at other intermediate values of *ζ*. In their study, the calculation of an integrated similarity function, assuming that the functional forms of *Φ*_{M} and *Φ*_{H} are valid for *ζ*≥1 for the whole interval (*ζ*_{0}, *ζ*), is questionable when *ζ*_{0}=*z*_{0}/*L*.

In this study, we propose a systematic mathematical analysis to estimate the upper bound for the applicability of these non-linear functional forms of similarity functions in very stable conditions. The analysis is performed by examining the behaviour of *C*_{D} and heat exchange coefficient (*C*_{H}) as functions of *ζ*. The methodology is applied to infer the upper bounds of *Ri*_{B} and *ζ* for which these similarity functions are applicable. The methodology is given in the general framework in the next section and is applied to different well-known similarity functions in §3. The theoretically derived upper bounds for the applicability of each non-linear similarity function are evaluated with the observations obtained from two experimental programmes in §4; §5 discusses the theoretical issues and the limitations of the analysis and §6 summarizes the conclusions drawn from this study.

## 2. Methodology

### (a) Theoretical background

Using MOS theory, the wind speed (*U*) and potential temperature (*θ*) profiles in the horizontally homogeneous and stationary surface layer are given by2.1and2.2where *u*_{*} is the friction velocity, *θ*_{*} the friction temperature, *z*_{0} the surface roughness length, *θ*_{0} the potential temperature at *z*_{0} and *k* the von Karman constant, and the stability parameter *ζ*=*z*/*L* (where *L* is the Obukhov length) is defined as2.3in which *g* is the acceleration due to gravity and is the mean potential temperature. Functions *ψ*_{M} and *ψ*_{H} are, respectively, integrated similarity functions of momentum and heat given by2.4and2.5where *ζ*′=*z*′/*L* and *ζ*_{0}=*z*_{0}/*L*=*z*_{0}/*zζ*. Note that the value of lower limit *ζ*_{0}=*z*_{0}/*L* (physically *z*_{0} is much smaller than *L* in the surface layer) in the derivation of *ψ*_{M} and *ψ*_{H} from the integration (equations (2.4) and (2.5)) is less than one.

Another stability parameter *Ri*_{B} is defined as2.6Using equations (2.1–2.3) in equation (2.6), one gets: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.* (2003*a*) 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 as2.8and2.9in which2.10and2.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 (*ζ*_{0}, *ζ*) 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*^{2}_{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*^{2}_{M} assumes a maximum value. Thus, the problem of minimization of *C*_{D} is equivalent to the maximization of *Y*^{2}_{M}. Accordingly, we are looking for the lowest value of *ζ*, which gives maximum *Y*^{2}_{M}. In equation (2.10), *z*_{0} and *z* are parameters, *L* is considered a variable and *ζ*_{0}=*z*_{0}/*zζ* is a function of *ζ* and thus, *Y*^{2}_{M} is actually a function of single variable *ζ*. The critical points of *Y*^{2}_{M} are obtained from2.12where the prime denotes the derivative of *Y*_{M} with respect to *ζ*.

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*^{2}_{M} are given by2.13As *ζ*>0 in stable conditions, we are primarily interested in the lowest positive root of equation (2.13). This root (say *ζ*_{crt}) is further investigated for maxima of *Y*^{2}_{M} (or minima of *C*_{D}) by ensuring that the second derivative of *Y*^{2}_{M} (i.e. ) with respect to *ζ* at *ζ*_{crt} is negative and *Y*^{2}_{M} increases in the range 0<*ζ*<*ζ*_{crt}. If at this critical point *ζ*=*ζ*_{crt}, it implies that this lowest root minimizes *Y*^{2}_{M} (or maximizes *C*_{D}), which is physically not desirable, and the functional form used for *Φ*_{M} in this case is not appropriate under stable conditions. The lowest root *ζ*=*ζ*_{crt} of equation (2.13), maximizing *Y*^{2}_{M} or minimizing *C*_{D}, is used in equation (2.7) to estimate the upper bound of *Ri*_{B} for applicability of the similarity function under stable conditions.

If equation (2.13) does not have any root, *Y*^{2}_{M} is either a continuously decreasing or an increasing function of *ζ*. For continuously decreasing *Y*^{2}_{M} or increasing *C*_{D}, the functional form of *Φ*_{M} will not be applicable under stable conditions. For continuously increasing *Y*^{2}_{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*_{M}*Y*_{H} in equation (2.9). The critical points of *Y*_{M}*Y*_{H} are given by2.14The maximum or minimum values of the function *Y*_{M} *Y*_{H} on these critical points can be investigated by the sign of second derivative of *Y*_{M} *Y*_{H} with respect to *ζ*. If equation (2.14) does not have any root, the function *Y*_{M} *Y*_{H} is either a continuously decreasing or an increasing function of *ζ*. For continuously decreasing *Y*_{M}*Y*_{H} or increasing behaviour of *C*_{H}, the functional forms of *Φ*_{M} or *Φ*_{H} will not be applicable under stable conditions. For continuously increasing *Y*_{M} *Y*_{H} or decreasing *C*_{H} for *ζ*>0, the functional forms of *Φ*_{M} and *Φ*_{H} can be used for all values of *ζ* up to infinity, implying that *C*_{H} approaches minimum as .

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

## 3. Applications

Now, we apply the theory discussed in §2 to infer the upper bounds of *ζ* and *Ri*_{B} for the applicability of each of the available non-linear similarity functions in the literature. As we are primarily interested in very stable conditions, the analysis has been carried out for *ζ*>1 in most of the cases. The non-linear functional forms of *Φ*_{M} and *Φ*_{H} proposed by various investigators are listed in tables 1 and 2, and their corresponding integrated similarity functions *ψ*_{M} and *ψ*_{H} obtained from equations (2.4) and (2.5) are given at the place where they are analysed.

### (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*_{0}. 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 as3.1where *β* 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.2For *Ψ*_{M} (equation (3.2)) corresponding to *ζ*≥1, equation (2.13) can be written as3.3Note that from equation (3.3), the critical point occurs only at *ζ*=*z*/*z*_{0}. At *ζ*=*z*/*z*_{0}, the second derivative *Y*_{M}′′=−*β*/*ζ*^{2}<0, implying that a local maxima of *Y*_{M} occurs at *ζ*=*z*/*z*_{0}. *Y*_{M}′ in equation (3.3) is positive for 1<*ζ*<*z*/*z*_{0} and becomes zero at *ζ*=*z*/*z*_{0}, 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*_{0}. Positively increasing behaviour of *Y*_{M} in this interval ensures that *Y*^{2}_{M} is also an increasing function and accordingly *C*_{D} is a decreasing function of *ζ* in the range from 1 to *z*/*z*_{0} (figure 1*a*). This decreasing behaviour of *C*_{D} for all values of *ζ* in the interval (0,*z*/*z*_{0}) is consistent with the physical behaviour of *C*_{D} with increasing stability in the stable atmosphere.

For *ζ*>*z*/*z*_{0}, *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*_{0}, *Y*_{M} decreases as *ζ* increases and vanishes at *ζ*=*ζ** (say) beyond which *Y*_{M} changes its sign and achieves all negative values. Therefore, *Y*^{2}_{M} decreases in the interval (*z*/*z*_{0},*ζ**) and approaches zero as *ζ* → *ζ** and then starts increasing continuously for *ζ*>*ζ**. This implies that *C*_{D} increases in the interval (*z*/*z*_{0},*ζ**), attaining its maximum value at *ζ* → *ζ** beyond which it decreases as *ζ* increases further (figure 1*a*). These increasing and decreasing behaviours of *C*_{D} for *ζ*>*z*/*z*_{0} 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*_{0} (where *C*_{D} has its first minimum) and is not of physical interest.

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

The upper value of *Ri*_{B} (i.e. *Ri*_{BU}) obtained from equation (2.7) at *ζ*=*z*/*z*_{0} is given by3.4Thus, 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*_{0} 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*_{0}, which is an order of magnitude larger than *z*_{0}, but at 1.5–2 times *h*_{0}. Thus, *z*/*z*_{0} 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*_{0} is expected to be satisfied for all the ratios of *z*/*z*_{0}.

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

#### (ii) Similarity function of Clarke (1970)

Clarke (1970) (hereafter CL-70) proposed a rational form of similarity function from Webb’s results to cover the whole boundary layer under stable conditions. The forms of *Φ*_{M} and *Φ*_{H} for CL-70 are given in table 1. Based on a mathematical analysis given in the electronic supplementary material, ESM-1.1, it is concluded that the function CL-70 is applicable only for *ζ*≤*ζ*_{c} and (table 1). It is observed that the upper bound of *Ri*_{B} for the applicability of CL-70 is finite and increases as the ratio *z*/*z*_{0} increases.

Note that the upper bound of *Ri*_{B} depends on the dimensionless ratio *z*/*z*_{0}. For illustration, we have taken *z*=10 m and *z*_{0}=0.1 m. Table 1 for CL-70 yields *ζ*_{c}=29.92 and *Ri*_{BU}=1.48. For these values of *z* and *z*_{0}, Carson & Richards (1978) pointed out that Clark’s formulation leads to and *Ri* as , whereas the present analysis reveals that *C*_{D} → 0 as . Further, they were not able to figure out from their analysis that *C*_{D} can take a minimum value for a finite *ζ* and *Ri*_{B}. However, they found graphically that *C*_{D} takes a minimum value at *Ri*_{B}=1.5.

These discrepancies in their analysis are attributed to incorrectly assuming the form of similarity function valid for *ζ*≥1 for the whole interval (*ζ*_{0},*ζ*) in the computation of integrated similarity functions *Ψ*_{M} and *Ψ*_{H}. Because of this, Carson & Richards (1978) concluded that the Clarke’s function does not give a physically realistic picture and cannot be used in the computation of surface fluxes in very stable conditions.

#### (iii) Similarity function of Hicks (1976)

In strong stability, Hicks (1976) (hereafter H-76) suggested a linear profile of similarity function above the height *ζ*>10 and recommended the forms of *Φ*_{M} and *Φ*_{H} given in table 1 for H-76. Based on a mathematical analysis of the behaviour of C_{D} in the electronic supplementary material, ESM-1.2, it is concluded that the Hicks function is applicable as long as (table 1):

*ζ*≤*ζ*_{H}and for 1<*z*/*z*_{0}≤*z*_{c}andall and

*Ri*_{B}<(1−*z*_{0}/*z*)(*c*−*βz*_{0}/*z*)^{−1}for*z*/*z*_{0}>*z*_{c},

where *z*_{c}≈7.593 and *ζ*_{H} is given in table 1 for H-76.

The Hicks function is purely empirical in nature, and its applicability is based on experimental observations. Micrometeorological measurements indicate that the ratio *z*/*z*_{0} should be larger than 10 in the surface layer (Sharan & Kumar 2008). Thus, the value of *z*/*z*_{0} lying in the range 1<*z*/*z*_{0}≤*z*_{c} is not of physical interest. Accordingly, Hicks formulation is applicable in very stable conditions as long as *Ri*_{B}<(1−*z*_{0}/*z*)(*c*−*βz*_{0}/*z*)^{−1}.

In contrast, Carson & Richards (1978) found in their analysis that the upper bound of *Ri*_{B} for Hicks function is [*c*(1−*z*_{0}/*z*)]^{−1}. The difference between the upper bounds of *Ri*_{B} is due to the discrepancy in their analysis attributed to wrongly applying the form of *Φ*_{M} valid only for *ζ*≥1 to the whole interval (*ζ*_{0}, *ζ*) in the computation of *Ψ*_{M}.

#### (iv) Similarity function of Holtslag & De Bruin (1988)

*Φ*_{M} proposed by Holtslag & De Bruin (1988) (hereafter HD-88) is given in table 1. Based on a mathematical analysis in the electronic supplementary material, ESM-1.3, it is shown that *C*_{D} continuously decreases with increasing values of *ζ* in stable conditions. This clearly states that the HD-88 function is applicable for all values of *ζ* in stable conditions (electronic supplementary material, ESM-1.3). The upper bound of *Ri*_{B} estimated from equation (2.7) by taking the limit as is 1/*a* (table 1). Thus, the similarity function of HD-88 is applicable for all positive values of *ζ* and *Ri*_{B}<1/*a* under stable conditions (table 1).

### (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} is3.5where *γ* 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 by3.6

Here, *Φ*_{M} is the same as in WE-70 (equation (3.1)), and thus it is applicable as long as *ζ*≤*z*/*z*_{0} (§3*a*(i)).

To estimate the lowest value of *ζ* for which *C*_{H} gets its first minimum value, equation (2.14) can be written as3.7Equation (3.7) has its lowest critical point at *ζ*=*z*/*z*_{0}. For 1<*ζ*<*z*/*z*_{0}, the derivative *Y*_{M}*Y*_{H}′+*Y*_{M}′*Y*_{H} of the denominator *Y*_{M}*Y*_{H} in *C*_{H} (equation (2.9)) given by equation (3.7) is positive and thus *Y*_{M}*Y*_{H} is an increasing function of *ζ*. For *ζ*>*z*/*z*_{0}, *Y*_{M}*Y*_{H}′+*Y*_{M}′*Y*_{H} is negative up to a point *ζ*=*ζ** (say) and thus *Y*_{M}*Y*_{H} will have the decreasing behaviour. Consequently, *C*_{H} decreases in 1<*ζ*<*z*/*z*_{0} and gets its first lowest value at *ζ*=*z*/*z*_{0} and then starts increasing for *ζ*>*z*/*z*_{0} (figure 1*b*). The decreasing behaviour of *C*_{H} in the range of *ζ* from 1 to *z*/*z*_{0} 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*_{0} and get their minimum values at *ζ*=*z*/*z*_{0}. The minimum of these two upper bounds of *ζ* for *C*_{D} and *C*_{H} will remain at the same point at *ζ*=*z*/*z*_{0}. The upper bound of *Ri*_{B} at this point *ζ*=*z*/*z*_{0} is given as3.8

Thus, B-82 is applicable as long as *ζ*≤*z*/*z*_{0} and .

#### (ii) Similarity function of Beljaars & Holtslag (1991)

Realizing the different transfer efficiencies between momentum and heat, Beljaars & Holtslag (1991) (hereafter BH-91) modified the HD-88 similarity function *Φ*_{M} (HD-88 in table 1) by simply changing the values of parameters as *a*=1 and *b*=0.667 and defined a different similarity function for *Φ*_{H} in the regime of intermittent turbulence. The *Φ*_{H} proposed by BH-91 is given in table 2. Based on a mathematical analysis in the electronic supplementary material, ESM-1.4 for BH-91, it is shown that the functional forms of *Φ*_{M} and *Φ*_{H} for BH-91 are applicable for all values of *ζ* and *Ri*_{B} in very stable conditions (table 2).

#### (iii) Similarity function of Cheng & Brutsaert (2005)

Cheng & Brutsaert (2005) (hereafter CB-05) analysed the wind and temperature profiles in the stable boundary layer by using the CASES-99 observations and proposed the non-linear forms for *Φ*_{M} and *Φ*_{H}, covering the entire range of stability from neutral to very stable conditions. *Φ*_{M} and *Φ*_{H} of CB-05 are given in table 2. The mathematical analysis in the electronic supplementary material, ESM-1.5 shows that the similarity function of CB-05 is theoretically valid for all values of *ζ* in stable conditions (table 2). The upper bound of *Ri*_{B} tends to infinity as in equation (2.7) and thus the function of CB-05 can be used for any positive value of *ζ* and *Ri*_{B}. This is consistent with the analysis of Guo & Zhang (2007).

#### (iv) Similarity function of Grachev *et al.* (2007)

A theoretical analysis is also carried out for the functional forms of *Φ*_{M} and *Φ*_{H} proposed by Grachev *et al.* (2007) (hereafter GR-07). These *Φ*_{M} and *Φ*_{H} (table 2) are based on the measurements of atmospheric turbulence made during the SHEBA. Although they have used the data corresponding to *Ri*_{B}≤0.2, these relations cover the entire range of stability from neutral to very stable conditions. In contrast, the present analysis reveals that the similarity function proposed by GR-07 is applicable for all values of *ζ* and *Ri*_{B} (table 2) in very stable conditions (electronic supplementary material, ESM-1.6).

## 4. Evaluation with field datasets

All the similarity functions used in the above analysis are purely empirical in nature, and, in general, their applicability and limitations are based on experimental observations. The theoretically derived upper bound for the applicability of each non-linear similarity function is evaluated with the observations obtained from two experimental programmes: CASES-99 (Poulos *et al.* 2002) and UK Meteorological Office’s Cardington monitoring facility (http://badc.nerc.ac.uk/data/cardington/).

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

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.

### (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 *θ*_{0} is taken as the measured mean grass infrared temperature. The roughness length *z*_{0} 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*_{0} 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.

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.

## 5. Theoretical issues and limitations

In this section, some underlying theoretical issues and limitations associated with the theoretical analysis described in the previous section are discussed.

### (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.* 2003*b*). The wind is considered as weak when (i) *U*_{s}≤2 m s^{−1} and (ii) *U*_{g}≤4 m s^{−1}; 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.* 2003*b*; 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*_{0}. In general, *z*_{0} 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.* 2003*a*) 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.

## 6. Conclusions

In very stable conditions, in which the stratified turbulence is weak and intermittent, the computation of surface fluxes by MOS theory with different linear/non-linear similarity functions becomes limited to specific ranges of *ζ* and *Ri*_{B}. A systematic mathematical analysis is presented to find the upper bounds of *ζ* and *Ri*_{B} for the extent of applicability of different non-linear similarity functions for non-dimensional wind and temperature profiles in the surface layer in these conditions. A generalized methodology is proposed on the basis of momentum drag coefficient *C*_{D} and heat exchange coefficient *C*_{H} and applied to various non-linear similarity functions available in the literature. The analysis is carried out for both kinds of similarity functions with *Φ*_{M}=*Φ*_{H} and *Φ*_{M}≠*Φ*_{H}.

The theoretical analysis for the criteria of each non-linear similarity function is evaluated with the observations obtained from CASES-99 and UK Meteorological Office’s Cardington datasets. Hence, these formulations for similarity functions can be used for the computation of surface fluxes in very stable conditions within their respective ranges of validity. The computed fluxes are scattered around the one-to-one line with the observations. The deviations are similar for all non-linear functions, implying that all these formulations for *Φ*_{M} and *Φ*_{H} lead to comparable fluxes. The proposed analysis clearly demonstrates the extent of applicability of some of the similarity functions such as those by Clarke (1970) and Hicks (1976) that are no longer in use in very stable conditions, and the fluxes computed from them are comparable to those obtained by other functions. The analysis proposed here is valid, in general, for finding the extent of applicability of any non-linear functional form of *Φ*_{M} and *Φ*_{H} in very stable conditions even with different roughness lengths of momentum and heat. This information is useful for accurately computing the surface fluxes in atmospheric models.

## Acknowledgements

The authors gratefully acknowledge the UK Meteorological Office for use of Cardington measurements and National Centre for Atmospheric Research (NCAR) for CASES-99 observations. The authors would like to thank the reviewers for their valuable comments and suggestions.

- Received April 27, 2010.
- Accepted July 9, 2010.

- This journal is © 2010 The Royal Society