Hybrid Profile–Gradient Approaches for the Estimation of Surface Fluxes

Abstract

The Monin–Obukhov similarity theory-based wind speed and potential temperature profiles are inherently coupled to each other. We have developed hybrid approaches to disentangle them, and as a direct consequence, the estimation of Obukhov length (and associated turbulent fluxes) from either wind-speed or temperature measurements becomes an effortless task. Additionally, our approaches give rise to two easily measurable indices of atmospheric stability. We compare these approaches with the traditional gradient and profile methods that require both wind-speed and temperature profile data. Using Monte-Carlo-type numerical experiments we demonstrate that, if the input profiles are free of any random errors, the performance of the proposed hybrid approaches is almost equivalent to the profile method and better than the gradient method. However, the proposed hybrid approaches are less competitive in comparison to their traditional counterparts in the presence of random errors.

1 Introduction

More than fifty years ago, in a classic paper, Panofsky (1963) wrote¹:

“In principle, it should be possible to determine the three parameters $z_0$ [aerodynamic roughness length], H [sensible heat flux], and $u_*$ [friction velocity] from three good wind observations close to the ground. But Priestley (1959) has pointed out that a small error in one or more of the winds leads to a huge error in the stress, so that this technique is not practical. Priestley further suggests that temperature data be added to the wind data in order that accurate estimates of stress be made. The present note considers this possibility in some detail.”

After this influential publication, the boundary-layer community at large embraced the idea and decided to focus on the estimation of turbulent fluxes utilizing both wind-speed and temperature data. The so-called gradient and profile methods (Appendix 1) were developed and refined. A few variants, using optimization techniques, were also proposed in parallel (e.g., Nieuwstadt 1978).

In contrast, only a handful of studies did not follow suit. Swinbank (1964), Klug (1967), and Lo (1979) explored the possibility of estimating turbulent fluxes using only wind-speed measurements. Even though they documented reasonably good results, their flux-estimation approaches never received any serious attention in the literature. After all these years, it is difficult to pin-point the exact reasons behind their unpopularity. It is plausible that the inherent complexities of the approaches by Klug (1967) and Lo (1979) utilizing numerical optimization techniques rendered them less desirable in practical applications. Klug’s approach also needed the aerodynamic roughness length ($z_0$) as an input, but accurate prescription of $z_0$ was (and still remains) a challenging task. The algorithm of Lo (1979) did not require $z_0$ as input, but suffered from convergence issues and possible mathematical errors (Zhang 1981). In addition, Lo (1979) did not include any error estimates of the derived variables as pointed out by Nieuwstadt and de Bruin (1981). The flux-estimation approach of Swinbank (1964) was more elegant, but was founded on the strong assumption that the surface-layer wind profile follows an exponential shape (Appendix 2). This assumption departed significantly from the well-accepted logarithmic form (with correction terms) for the wind profile, which likely contributed to its unpopularity.

With the advent of high-resolution, high-accuracy instruments for the measurement of wind speed and temperature (e.g., sodars, lidars, distributed temperature sensors), it is worthwhile to revisit the assertions made by Panofsky (1963). The argument that we need both wind-speed and temperature measurements for flux estimation may no longer be tenable. At the same time, one needs to have a more analytically tractable approach than that advocated by Lo (1979) or Klug (1967). Recently, in a short communication, we proposed such an approach, called the hybrid-wind approach (Basu 2018). With a few mathematical manipulations, we demonstrated that it is actually very straightforward to estimate turbulent fluxes from only wind-speed measurements. Our hybrid approach is similar to Swinbank (1964). In the present study, we first extend this approach to utilize temperature data as input. Next, we compare the proposed hybrid approaches against traditional gradient and profile methods for a wide range of stability conditions. Last and most importantly, through uncertainty propagation experiments, we quantify the errors in estimated fluxes from all the aforementioned approaches.

The structure of the paper is as follows. Sect. 2 introduces the newly proposed hybrid flux-estimation approach, and as by-products of this approach, two atmospheric stability indices are derived. Their characteristics are discussed in Sect. 3. Some caveats of the proposed hybrid approaches are touched upon in Sect. 4, and illustrative examples comparing the proposed approach and traditional flux-estimation approaches are documented in Sect. 5. The uncertainty propagation experiments and the associated results are also elaborated in this section. The concluding remarks including future directions are summarized in Sect. 6. Background information on the traditional flux-estimation approaches, Swinbank’s exponential wind-profile equation, and several relevant stability correction formulations are provided in the Appendices.

2 Methodology

The surface-layer wind speed and potential temperature profile equations based on the Monin–Obukhov similarity theory (MOST; Monin and Obukhov 1954) are written as,

$$\begin{aligned} U\left( z\right)= & {} \frac{u_*}{\kappa } \left[ \ln \left( \frac{z}{z_0} \right) - \psi _m\left( \frac{z}{L}\right) + \psi _m\left( \frac{z_0}{L}\right) \right] , \end{aligned}$$

(1a)

$$\begin{aligned} \Theta \left( z\right) - \Theta _S= & {} \frac{\theta _*}{\kappa } \left[ \ln \left( \frac{z}{z_{0 T}} \right) - \psi _h\left( \frac{z}{L}\right) + \psi _h\left( \frac{z_{0 T}}{L}\right) \right] , \end{aligned}$$

(1b)

where $\psi _m$ and $\psi _h$ are stability correction terms; $u_*$, $\theta _*$, and L denote friction velocity, surface temperature scale, and Obukhov length, respectively. The aerodynamic roughness length and roughness length for temperature are represented by $z_0$ and $z_{0 T}$, respectively. $\Theta _S$ is the surface temperature, while the von Kármán constant is denoted by $\kappa$.

Based on Eq. 1a, the vertical wind-speed difference (aka increment) can be computed as follows,

$$\begin{aligned} \Delta U_{21}= & {} U\left( z_2\right) - U\left( z_1\right) = \frac{u_*}{\kappa } \left[ \ln \left( \frac{z_2}{z_1} \right) - \psi _m\left( \frac{z_2}{L}\right) + \psi _m\left( \frac{z_1}{L}\right) \right] , \end{aligned}$$

(2a)

$$\begin{aligned} \Delta U_{31}= & {} U\left( z_3\right) - U\left( z_1\right) = \frac{u_*}{\kappa } \left[ \ln \left( \frac{z_3}{z_1} \right) - \psi _m\left( \frac{z_3}{L}\right) + \psi _m\left( \frac{z_1}{L}\right) \right] , \end{aligned}$$

(2b)

where $z_1$, $z_2$, and $z_3$, are the heights at which wind speed is measured.

Finally, a ratio of the wind-speed differences can be written as,

$$\begin{aligned} R_W = \frac{\Delta U_{31}}{\Delta U_{21}} = \frac{\ln \left( \frac{z_3}{z_1} \right) - \psi _m\left( \frac{z_3}{L}\right) + \psi _m\left( \frac{z_1}{L}\right) }{\ln \left( \frac{z_2}{z_1} \right) - \psi _m\left( \frac{z_2}{L}\right) + \psi _m\left( \frac{z_1}{L}\right) }. \end{aligned}$$

(3)

In an analogous manner, a ratio of the potential temperature differences can be written as,

$$\begin{aligned} R_T = \frac{\Delta \Theta _{31}}{\Delta \Theta _{21}} = \frac{\ln \left( \frac{z_3}{z_1} \right) - \psi _h\left( \frac{z_3}{L}\right) + \psi _h\left( \frac{z_1}{L}\right) }{\ln \left( \frac{z_2}{z_1} \right) - \psi _h\left( \frac{z_2}{L}\right) + \psi _h\left( \frac{z_1}{L}\right) }. \end{aligned}$$

(4)

We strongly emphasize that the estimation of $R_W$ only requires observed wind-speed data from three levels; similarly, $R_T$ is solely based on temperature measurements at three levels. Due to their explicit functional relationships with L, both these quantities can be considered as independent proxies of atmospheric stability. In other words, both the wind-speed and temperature profile data are not required for the estimations of L and associated fluxes; only one type of variable suffices. Illustrative examples are provided in Sect. 5.

We have named our flux-estimation methodology a ‘hybrid’ profile–gradient approach because it borrows ideas from both the traditional profile and gradient methods. Via mathematical manipulations, it disentangles the original MOST equations, which has not been feasible in the traditional methods. Hereafter, we make a further distinction and refer to the proposed approach as ‘hybrid-W’ or ‘hybrid-T’ depending on whether wind-speed or temperature data are being utilized as inputs.

3 Characteristics of $R_W$ and $R_T$

The behaviour of $R_W$ and $R_T$ depends entirely on the stability correction terms ($\psi _m$ and $\psi _h$). For neutral conditions (i.e., $z/L = 0$), $\psi _m = \psi _h = 0$, whence, both ratios simplify to

$$\begin{aligned} R_N = \frac{\ln \left( \frac{z_3}{z_1} \right) }{\ln \left( \frac{z_2}{z_1} \right) }. \end{aligned}$$

(5)

If $z_3> z_2 > z_1$, it is trivial to show that $R_N > 1$. Next, we consider the behaviour of $R_W$ and $R_T$ for non-neutral conditions.

In Fig. 1, the variations of these ratios with respect to 1 / L are shown, where several well-known $\psi _m$ and $\psi _h$ functions are utilized in these plots. More details on these functions can be found in Appendix 3. In these plots, the sensor heights are assumed to be at $z_1 = 5\ \hbox {m}$, $z_2 = 10\ \hbox {m}$, and $z_3 = 20\ \hbox {m}$, respectively. For these specific heights, $R_N = 2$. Clearly, for unstable conditions (left panel), both $R_W$ and $R_T$ monotonically decrease with increasing instability. In contrast, for stable conditions (right panels), these ratios show a monotonically increasing trend with increase in stability. For Businger–Dyer functions (Dyer and Hicks 1970; Businger et al. 1971; Dyer 1974), it can be readily deduced that both $R_W$ and $R_T$ should approach constant values under very stable conditions:

$$\begin{aligned} R_{VS} \rightarrow \frac{\left( z_3 - z_1\right) }{\left( z_2 - z_1\right) }. \end{aligned}$$

(6)

For the chosen sensor heights, $R_{VS} = 3$. This asymptotic behaviour is prominently evident in the right panels of Fig. 1.

Fig. 1

Variations of $R_W$ (top panel) and $R_T$ (bottom panel) with respect to inverse Obukhov length (1 / L). The left and right panels represent unstable and stable conditions, respectively. The legends in these plots correspond to the selected stability correction functions

In summary, for the selected stability correction functions, $R_W$ and $R_T$ are single-valued functions of L. Thus, it should be straightforward to estimate L given measured value of either $R_W$ and $R_T$. In this regard, any suitable root-finding algorithm (e.g., Newton–Raphson approach) can be utilized; we make use of the well-known Levenberg–Marquardt algorithm. Once L is estimated, one can estimate $u_*$ from Eqs. 2a and 2b. Since there are two equations and only one unknown, the conventional linear regression approach with ordinary least squares can be employed. Having determined both L and $u_*$, one can then estimate $\overline{w\theta }$ from the definition of Obukhov length. A similar strategy can be followed in conjunction with $R_T$ as input. Of course, in this case, one solves for $\theta _*$ instead of $u_*$, and from the definition of L, one deduces $u_*$, and subsequently, $\overline{w\theta }$.

4 Limitations of the Proposed Hybrid Approaches

Before delving into the results, we would like to mention a few issues that may limit the applications of the proposed hybrid approaches.

4.1 Validity of MOST

Both the hybrid-W and hybrid-T approaches are deeply rooted in MOST. Hence, they are only applicable when and where MOST is applicable. We would like to remind the readers that MOST is strictly valid in a horizontally homogeneous surface layer. In the surface layer (aka constant flux layer), the turbulent fluxes are assumed to be invariant with height. Thus, all the sensor heights (i.e., $z_1$ , $z_2$, $z_3$) should be within the surface layer to avoid violation of MOST. For strongly stratified conditions, the surface layer may be only a few metres deep; the proposed hybrid approaches should be avoided under that scenario.

4.2 Monotonicity of Input Mean Profiles

The hybrid-W approach implicitly assumes that wind speeds monotonically increase with height. Similarly, in the case of the hybrid-T approach, the potential temperature is expected to monotonically increase (decrease) with height for stable (unstable) conditions. If such monotonic conditions are not met, the proposed approaches should not be used.

4.3 Similarity of Footprints

The footprints for scalars and fluxes should be similar in order to estimate fluxes accurately via MOST; under homogeneous surface conditions, this restriction is not that important. However, for heterogeneous cases, the mismatch of footprints could pose a serious limitation. Of course, any application of MOST for these cases will also be questionable.

4.4 Multi-valued Functions

In Sect. 3, we have shown that $R_W$ and $R_T$ variables are single-valued functions of L for a specific set of widely used stability correction functions. However, there are exceptions. In Fig. 2, we compute the same ratios using stability correction formulations proposed by Beljaars and Holtslag (1991) and Cheng and Brutsaert (2005) for stably-stratified conditions (see Appendix 3 for details). Clearly, the resultant functions are multi-valued; in other words, given $R_W$ or $R_T$ , it is not possible to estimate unique values of L. As a consequence, our proposed hybrid approach should not be used in conjunction with these specific stability correction functions.

Fig. 2

Variations of $R_W$ (left panel) and $R_T$ (right panel) with respect to inverse Obukhov length (1 / L). The legends in these plots correspond to the selected stability correction functions

4.5 Turbulent Prandtl Number

In the MOST relation for the potential temperature profile, Eq. 1b, we implicitly assume that the turbulent Prandtl number ($Pr_T$) is equal to one. Since the estimation of L only depends on the ratio $R_T$, this assumption is not relevant. However, its influence on the estimations of $\theta _*$ and $u_*$ via the hybrid-T approach cannot be disregarded. Note that the hybrid-W approach does not involve any information about $Pr_T$.

4.6 Effects of Moisture

Throughout this paper, we have only considered dry atmospheric conditions in the surface layer. It is, however, straightforward to extend the hybrid approaches for moist conditions (e.g., offshore environments). In these cases, one must utilize virtual kinematic heat flux and the virtual potential temperature in the definition of Obukhov length (L) and in Eq. 1b. The stability parameter (z / L) can even be partitioned to account for sensible heat flux and latent heat flux separately. For further details, see Barthelmie et al. (2010) and the references therein.

5 Inter-Comparison of Different Flux-Estimation Approaches

In order to compare the proposed hybrid approaches against the traditional ones, we perform Monte-Carlo-type numerical experiments with the following steps:

• (i) To encompass a wide-range of stability conditions, we assume $u_* \in \left[ 0.1~2\right]$ m s$^{-1}$ and $\theta _* \in \left[ -\,1~0.2\right]$ K. From these sets, we randomly (with uniform probability) select a $\left( u_*, \theta _*\right)$ pair.

• (ii) Furthermore, we assume $z_0 = z_{0 T} = 0.1\ \hbox {m}$ and $\Theta _s = \Theta _0 = 300\ \hbox {K}$.

• (iii) Using these selected inputs, we first estimate L, and then in turn, predict $U\left( z\right)$ and $\Theta \left( z\right)$ via Eqs. 1a and 1b in conjunctions with the Businger–Dyer stability correction functions [i.e., Eqs. 13a–13c].

• (iv) In ‘noise-free input data’ cases, we skip this specific step. Otherwise, we add random noise on the U(z) and $\Theta (z)$ profiles. More details on the characteristics of additive noise are provided later.

• (v) If the estimated $z/|L| < 1$ and mean wind speed $> 1 \ \hbox {m s}^{-1}$, then we proceed to the following step. Otherwise, we discard the selected $\left( u_*, \theta _*\right)$ pair and go back to the first step. In the ‘noisy input data’ cases, we enforce a few more additional exclusion criteria which will be discussed later.

• (vi) Next, we attempt to do the following inverse computation: given the predicted mean wind-speed and/or temperature profiles, can we accurately estimate the surface fluxes? In hybrid-W (hybrid-T) approach, we estimate the surface fluxes by only using wind-speed (potential temperature) data from $z = 5$, 10, and 20 m.

• (vii) In order to have a direct comparison, we also estimate fluxes using the traditional gradient and profile methods (Appendix 1). In this case, both wind and temperature data from the lowest two levels are utilized.

• (viii) For all the flux-estimation approaches, we quantify the relative errors in the estimations of $u_*$ and $\theta _*$.

• (ix) We repeat all the previous steps until we get $10^5$ admissible samples for all the scenarios.

5.1 Noise-Free Input Data

The relative errors for the estimations of $u_*$ and $\theta _*$ are reported in Table 1. These errors are computed as follows:

$$\begin{aligned} RE = \frac{\chi _{\mathrm{est}}-\chi _{\mathrm{true}}}{\chi _{\mathrm{true}}}\times 100, \end{aligned}$$

(7)

where $\chi$ is either $u_*$ or $\theta _*$. In addition to minimum and maximum values, several percentiles (based on $10^5$ samples for each case) are reported in Table 1.

Table 1 Relative errors (%) in the estimations of $u_*$ and $\theta _*$

Clearly, for both $u_*$ and $\theta _*$, the performance of the traditional profile method is the best among all the approaches as it leads to null errors. In contrast, the traditional gradient method seems to suffer from a systematic error of $O(4\%)$. This error stems from finite-difference approximations, as discussed by Arya (1991).

For both hybrid approaches, the relative errors equal zero for percentiles ranging from 1 to 99. In the case of the hybrid-W approach, negligible errors can occur in the estimation of $u_*$ due to round-off errors during the optimization process. In the case of $\theta _*$, only 17 samples (out of $10^5$) exceeded errors > 1%. Most of these cases had true $\theta _*$ values close to zero and the division by a small number led to very large relative errors. The performance of the hybrid-T approach was perfect for the estimation of $\theta _*$. In the case of $u_*$ estimation, 16 samples (out of $10^5$) exceeded absolute relative error of 1%. In summary, for the noise-free cases, the overall performance of the proposed hybrid approaches is almost at par with the traditional profile method. In the following subsection, we investigate if this conclusion holds in the presence of random errors in input mean profiles.

5.2 Noisy Input Data

We conduct uncertainty propagation experiments to quantify if and how the errors in the input profiles are amplified during various flux estimations. We first add different amounts of noise to the profiles as follows,

$$\begin{aligned} \widetilde{U}= & {} U + \eta _U, \end{aligned}$$

(8a)

$$\begin{aligned} \widetilde{\Theta }= & {} \Theta + \eta _\Theta . \end{aligned}$$

(8b)

The noise terms ($\eta _U$ and $\eta _\Theta$) are generated from a multivariate Gaussian distribution with zero mean and the following covariance matrix,

$$\begin{aligned} \Sigma = \sigma ^2\begin{bmatrix} 1&\rho&\rho \\ \rho&1&\rho \\ \rho&\rho&1 \end{bmatrix}, \end{aligned}$$

(9)

where $\sigma ^2$ is the variance of the noise term. Since we are only concerned with three levels of observations, $\Sigma$ is a $3\times 3$ matrix. The variable $\rho$ captures the correlation of noise between different levels. Such a correlated noise situation is possible when a single instrument (e.g., a lidar) is used to measure wind speeds (or temperature) at different heights.

Table 2 Different scenarios for the noise terms

We consider several noise scenarios which are listed in Table 2. Specifically, we consider two noise levels (with appropriate units): 0.01 (low) and 0.05 (high). In addition, two values of $\rho$ are considered: 0.9 (high) and 0.5 (low). Since the hybrid-W approach only requires wind-speed data, please note that the scenarios 4, 5, and 6 are all the same for this approach.

Illustrative noise values ($\eta _U$) are shown in Fig. 3. Clearly, for $\sigma = 0.05$ m s$^{-1}$, the noise terms can reach up to $\pm \, 0.2$ m s$^{-1}$. A large amount of additive random noise can distort the U(z) and $\Theta (z)$ profiles significantly and can even make them physically unrealistic. To avoid such undesirable situations, we implemented certain exclusion criteria in addition to the ones discussed in the previous sub-section [i.e., step (v)]. If the noisy U(z) and $\Theta (z)$ profiles are not monotonic, we exclude that particular case. If the resultant $R_W$ and $R_T$ values are outside their acceptable ranges (i.e., $1.8< R_W < 3$ and $1.7< R_T < 3$; see Fig. 1), those cases are also excluded.

Fig. 3

Bi-variate probability density functions for scenarios 3 (left panel) and 4 (right panel), respectively. For both the scenarios, $\sigma = 0.05\ \hbox {m s}^{-1}$. However, $\rho$ changes from 0.9 in scenario 3–0.5 in scenario 4. Here $z_1$ and $z_2$ denote two different sensor heights

Fig. 4

Absolute relative errors in the estimation of $u_*$ for four scenarios: 1 (top-left panel), 2 (top-right panel), 3 (bottom-left panel), and 4 (bottom-right panel). Flux-estimation approach: hybrid-W

Fig. 5

Absolute relative errors in the estimation of $\theta _*$ for four scenarios: 1 (top-left panel), 2 (top-right panel), 3 (bottom-left panel), and 4 (bottom-right panel). Flux-estimation approach: hybrid-W

Fig. 6

Absolute relative errors in the estimation of $u_*$ (top panels) and $\theta _*$ (bottom panels) for two scenarios: 5 (left panels), and 6 (right panels). Flux-estimation approach: hybrid-T

Fig. 7

Absolute relative errors in the estimation of $u_*$ (top panels) and $\theta _*$ (bottom panels) for two scenarios: 5 (left panels), and 6 (right panels). Flux-estimation approach: traditional gradient method

Fig. 8

Absolute relative errors in the estimation of $u_*$ (top panels) and $\theta _*$ (bottom panels) for two scenarios: 5 (left panels), and 6 (right panels). Flux-estimation approach: traditional profile method

The results from our uncertainty propagation experiments are shown in Figs. 4, 5, 6, 7, and 8. In these figures, we report various percentiles of absolute relative errors for both $u_*$ and $\theta _*$. The summary of our results is as follows:

• Hybrid-W: for scenarios 1 and 2, the error in $u_*$ estimation is less than 10%. However, the errors increase substantially for scenarios 3 and 4. For low $u_*$ values, the errors can range from 10 to 100%; however, for high $u_*$ values, they are mostly less than 10%. The performance of this approach for $\theta _*$ estimation is somewhat poorer. For stable conditions, the median absolute error values are largely on the order of 10–20%. For unstable conditions, they are higher and seem to be independent of $\theta _*$ values. For near-neutral conditions, large errors can occur due to the division by small numbers.

• Hybrid-T: the estimation of $\theta _*$ is far better than $u_*$ for both scenarios 5 and 6. For unstable conditions, the median error values in $\theta _*$ are largely less than 20%. Marginally higher errors are noticeable for stable conditions.

• Gradient: for scenarios 5 and 6, for low $u_*$ values, the errors could be on the order of 10–100%. Otherwise, for high $u_*$ values, they are much lower than 10%. For all conditions (with the exception of near-neutral), $\theta _*$ errors are less than 10%.

• Profile: similar to the noise-free cases, this approach outperforms others in both the scenarios 5 and 6. Qualitatively, the errors in $u_*$ estimation follow a similar trend as the hybrid-W approach. However, the magnitude of the errors are much smaller. The errors in the estimation of $\theta _*$ also barely exceed 10–20% (other than the near-neutral conditions).

Before closing, we want to stress that our findings from these uncertainty propagation experiments should be used with caution. We selected specific types of additive noise which are correlated across different heights. Other alternatives are also possible. For example, we used a fixed $\sigma$ value for a given scenario; instead, one could use $\sigma$ dependent on the magnitude of U or $\Theta$. In that case, the trends reported in Figs. 4, 5, 6, 7, and 8 would be significantly different. Furthermore, we assumed that the noise in wind-speed and potential temperature profiles are uncorrelated; we do not know if this assumption is realistic or not. In general, high wind speeds lead to lower temperature measurement (radiation) errors; thus, the random errors in wind speeds and temperature might be (anti) correlated.

6 Concluding Remarks

We have developed new approaches to estimate surface fluxes utilizing either wind-speed or temperature profile data. We have compared our approaches against traditional gradient and profile methods that require both wind-speed and temperature profile data. For noise-free input data, the hybrid approaches perform as well as the traditional profile method. However, in the presence of random errors in input data, the proposed approaches lead to somewhat more flux-estimation errors than the traditional ones.

Given the unique one-to-one relationships between the ratio of wind-speed differences (or the ratio of potential temperature differences) with the Obukhov length, we propose that either of these ratios could be utilized as a proxy for atmospheric stability. In Basu (2018), we demonstrated that the ratio of wind-speed differences was able to categorize observational data in a physically meaningful way. However, further direct verifications are needed.

We believe that the hybrid-W approach is ideally suited for sodar and lidar-based wind-speed measurements owing to their high vertical resolution in the surface layer. Similarly, the distributed temperature sensing-based high-resolution temperature profiles can be utilized as inputs for the hybrid-T approach. In our future work, observational datasets from various field campaigns will be utilized to make an in-depth assessment of the proposed hybrid-W and hybrid-T approaches. Of course, we will pay close attention to the issues of non-stationarity and heterogeneity, as under such circumstances, the usage of the proposed hybrid approaches (and MOST in general) is not appropriate.

7 Appendix 1: Traditional Gradient and Profile Methods

In the traditional gradient method, the following normalized gradient equations are solved in a coupled and iterative manner (Arya 2001),

$$\begin{aligned} \left( \frac{\kappa z}{u_*}\right) \left( \frac{\partial U}{\partial z}\right)= & {} \phi _m\left( \frac{z}{L} \right) , \end{aligned}$$

(10a)

$$\begin{aligned} \left( \frac{\kappa z}{\theta _*}\right) \left( \frac{\partial \Theta }{\partial z}\right)= & {} \phi _h\left( \frac{z}{L} \right) . \end{aligned}$$

(10b)

The vertical gradients are approximated by the finite-difference formulation as follows: $\frac{\partial U}{\partial z} \approx \frac{\Delta U}{\Delta z} = \frac{U(z_2) - U(z_1)}{\left( z_2 - z_1\right) }$, and $\frac{\partial \Theta }{\partial z} \approx \frac{\Delta \Theta }{\Delta z} = \frac{\Theta (z_2) - \Theta (z_1)}{\left( z_2 - z_1\right) }$. The estimated gradients are applicable at the mid-point height $z_m = \frac{z_1+z_2}{2}$. Even though this approach (based on linear approximation) is the most popular, an alternative approach utilizing logarithmic approximation was proposed by Arya (1991). For unstable (stable) conditions, the logarithmic (linear) approximation-based approach was found to outperform its counterpart.

Application of the profile method typically requires the following variables as input: wind speed at one level, temperature at two levels, and aerodynamic roughness length (Berkowicz and Prahm 1982). In a slightly modified version, one uses wind-speed from an additional level instead of the roughness length. One then utilizes the MOST-based profile equations and solves for the unknown fluxes. Brotzge and Crawford (2000) utilized this modified profile approach to estimate fluxes from the Oklahoma mesonet.

8 Appendix 2: Swinbank’s Exponential Wind Profile

Swinbank (1964) proposed the following equation for surface-layer wind profile,

$$\begin{aligned} U\left( z_2\right) - U\left( z_1\right) = \frac{u_*}{\kappa }\ln \left[ \frac{\exp \left( \frac{z_2}{L}\right) - 1}{\exp \left( \frac{z_1}{L}\right) - 1}\right] , \end{aligned}$$

(11)

and further derived

$$\begin{aligned} \frac{U\left( z_3\right) - U\left( z_1\right) }{U\left( z_2\right) - U\left( z_1\right) } = \frac{\ln \left[ \frac{\exp \left( \frac{z_3}{L}\right) - 1}{\exp \left( \frac{z_1}{L}\right) - 1}\right] }{\ln \left[ \frac{\exp \left( \frac{z_2}{L}\right) - 1}{\exp \left( \frac{z_1}{L}\right) - 1}\right] }, \end{aligned}$$

(12)

commenting that Eq. 12 permits the determination of L from observed wind-speed data at three levels using numerical or graphical interpolation. Once L is determined, $u_*$ can be estimated from Eq. 11. Our proposed hybrid-W approach is almost identical, albeit it makes use of Eq. 3.

9 Appendix 3: Stability Correction Functions

Over the years, numerous stability correction functions have been proposed in the literature. A few of them are listed below:

Dyer and Hicks (1970), Businger et al. (1971), Dyer (1974):

$$\begin{aligned} \psi _m= & {} 2\ln \left( \frac{1+x}{2} \right) + \ln \left( \frac{1+x^2}{2} \right) - 2\tan ^{-1}x + \frac{\pi }{2}; \quad \text{ for } \frac{z}{L} \le 0 \end{aligned}$$

(13a)

$$\begin{aligned} \psi _h= & {} 2\ln \left( \frac{1+x^2}{2} \right) ; \quad \text{ for } \frac{z}{L} \le 0 \end{aligned}$$

(13b)

$$\begin{aligned} \psi _m= & {} \psi _h = -5\frac{z}{L}; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(13c)

where $x = \left( 1 - 16\frac{z}{L}\right) ^{1/4}$.

Beljaars and Holtslag (1991):

$$\begin{aligned} \psi _m= & {} -a\frac{z}{L} -b\left( \frac{z}{L} - \frac{c}{d} \right) \exp \left( -d\frac{z}{L} \right) -\frac{bc}{d}; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(14a)

$$\begin{aligned} \psi _h= & {} -\left( 1 + \frac{2a}{3}\frac{z}{L} \right) ^{3/2} -b\left( \frac{z}{L}-\frac{c}{d}\right) \exp \left( -d\frac{z}{L} \right) -\frac{bc}{d}+1; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(14b)

where $a = 1$, $b = \frac{2}{3}$, $c = 5$, and $d = 0.35$.

Duynkerke (1991):

$$\begin{aligned} \psi _m= & {} -\left( 1+ \frac{\beta _m}{\alpha _m}\frac{z}{L}\right) ^{\alpha _m}; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(15a)

$$\begin{aligned} \psi _h= & {} -\left( 1+ \frac{\beta _h}{\alpha _h} \frac{z}{L}\right) ^{\alpha _h}; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(15b)

where $\alpha _m = \alpha _h = 0.8$, $\beta _m = 5$, and $\beta _h = 7.5$.

Wilson (2001):

$$\begin{aligned} \psi _m= & {} 3\ln \left( \frac{1+\sqrt{1+\gamma _m {|}z/L|^{2/3}}}{1+\sqrt{1+\gamma _m |z_0/L|^{2/3}}} \right) ; \quad \text{ for } \frac{z}{L} \le 0 \end{aligned}$$

(16a)

$$\begin{aligned} \psi _h= & {} 3\ln \left( \frac{1+\sqrt{1+\gamma _h |z/L|^{2/3}}}{1+\sqrt{1+\gamma _h |z_{0 T}/L|^{2/3}}} \right) ; \quad \text{ for } \frac{z}{L} \le 0 \end{aligned}$$

(16b)

where $\gamma _m = 3.6$ and $\gamma _h = 7.9$.

Cheng and Brutsaert (2005):

$$\begin{aligned} \psi _m= & {} -a\ln \left( \frac{z}{L} + \left( 1 + \left( \frac{z}{L}\right) ^b\right) ^{1/b}\right) ; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(17a)

$$\begin{aligned} \psi _h= & {} -c\ln \left( \frac{z}{L} + \left( 1 + \left( \frac{z}{L}\right) ^d\right) ^{1/d}\right) ; \quad \text{ for } \frac{z}{L} \ge 0 \end{aligned}$$

(17b)

where a, b, c,  and d equal to 6.1, 2.5, 5.3, and 1.1, respectively.