Surface-Layer Similarity Functions for Dissipation Rate and Structure Parameters of Temperature and Humidity Based on Eleven Field Experiments

In the literature, no consensus can be found on the exact form of the universal funtions of Monin-Obukhov similarity theory (MOST) for the structure parameters of temperature, CT2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_T}^2$$\end{document}, and humidity, Cq2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_q}^2$$\end{document}, and the dissipation rate of turbulent kinetic energy, ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}. By combining 11 datasets and applying data treatment with spectral data filtering and error-weighted curve-fitting we first derived robust MOST functions of CT2,Cq2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_T}^2, {C_q}^2$$\end{document} and ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} that cover a large stability range for both unstable and stable conditions. Second, as all data were gathered with the same instrumentation and were processed in the same way—in contrast to earlier studies—we were able to investigate the similarity of MOST functions across different datasets by defining MOST functions for all datasets individually. For CT2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_T}^2$$\end{document} and ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} we found no substantial differences in MOST functions for datasets over different surface types or moisture regimes. MOST functions of Cq2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_q}^2$$\end{document} differ from that of CT2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_T}^2$$\end{document}, but we could not relate these differences to turbulence parameters often associated with non-local effects. Furthermore, we showed that limited stability ranges and a limited number of data points are plausible reasons for variations of MOST functions in the literature. Last, we investigated the sensitivity of fluxes to the uncertainty of MOST functions. We provide an overview of the uncertainty range for MOST functions of CT2,Cq2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_T}^2, {C_q}^2$$\end{document} and ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}, and suggest their use in determining the uncertainty in surface fluxes.


Introduction
Monin-Obukhov similarity theory (MOST; Monin and Obukhov 1954) is used to describe turbulence characteristics in the atmospheric surface layer (Stull 1988). MOST applies to statistics (e.g. variances, vertical gradients or structure parameters) of scalars such as temperature (T ), specific humidity (q), the three velocity components (u,v,w) and the dissipation rate of turbulent kinetic energy (ε). The application of MOST allows measured turbulence statistics to be related to temperature, humidity and wind scales, which are defined in terms of the turbulent fluxes of heat (H ), moisture (L v E), momentum (τ ), the Obukhov length scale (L) and the height above the surface (z). Through these relations (MOST functions) it is possible to determine turbulent fluxes from measured turbulence statistics once the MOST functions are known. These MOST functions are assumed to be universal, which means that once they are determined they are equal for all similar circumstances.
Here we focus on MOST functions that are relevant to scintillometry, notably that of the structure parameters of temperature (C T 2 ) and humidity (C q 2 ), and ε. Scintillometers measure the path-averaged turbulence intensity of the refractive index of air, which is mainly dependent on temperature and humidity fluctuations (e.g. De Bruin et al. 1995). The sensitivity of a scintillometer for these temperature and/or humidity fluctuations is dependent on the wavelength at which it operates; C T 2 can be obtained from an optical scintillometer, and C q 2 can be acquired from combined optical and microwave scintillometers (e.g. Andreas 1989; Lüdi et al. 2005;Meijninger et al. 2006). Furthermore, ε can be obtained with a smallaperture scintillometer (Hartogensis et al. 2002). The main application that we have in mind is to obtain H, L v E and τ using C T 2 , C q 2 and ε determined from scintillometer data. Several MOST functions for C T 2 , C q 2 and ε have been reported in the literature, of which an overview is given in Online Appendix. 1 This overview illustrates that MOST functions of the same parameter vary largely across different studies. The individual functions have been derived from different measurement instruments, environmental and meteorological conditions and data analysis procedures of separate datasets, as listed in Braam et al. (2014). As all of these factors may have contributed to the variations in the MOST functions, it cannot be concluded from the overview in the online Appendix which factors have caused the variations and if MOST functions can actually be applied across different datasets. Hartogensis et al. (2003) made an inventory of the uncertainty of all parameters that are used in the scintillometer heat-flux calculations and showed that differences in MOST functions cause up to 20 % difference in H at the free convection limit. This difference is substantial and can have a considerable influence on intercomparison studies of different measurement techniques. This clearly demonstrates the need to determine well-defined MOST relations that do not suffer from the issues listed above and that consensus is needed on the MOST functions to be used in flux calculations.
Besides the large variability in the MOST functions, earlier studies found that MOST functions of C T 2 and C q 2 are not always equal (Li et al. 2012;Maronga 2014). This dissimilarity is in contradiction to the assumptions of MOST that heat and moisture are transported by the same mechanism and therefore should have the same MOST function (De Bruin et al. 1993). MOST functions of C T 2 and C q 2 need to be considered separately to investigate if the dissimilarity between MOST functions of C T 2 and C q 2 is consistent for different datasets. The main goal of our study is to define robust MOST functions for C T 2 , C q 2 and ε. These MOST functions are based on 11 datasets that together represent a large stability range, and were gathered with the same instrumentation and processed using the same techniques. With this approach we overcome the shortcomings of previous studies in which MOST functions were obtained from a range of instrumentation and different processing techniques. Besides this main goal our aims include the following: -The large amount of data over a wide stability range allows us to determine a sound uncertainty range of MOST parameters, which can be used in estimating the uncertainty in the calculated fluxes. -We define MOST functions from all individual datasets. Based on the spread in these functions we then discuss plausible reasons for the large range of functions reported in the literature. As part of this we also investigate if the (dis)similarity between MOST functions of C T 2 and C q 2 differs between datasets. -We discuss the validity of the derived MOST relations in two ways. First we show if and to what extent the MOST relations are influenced by cross-correlation of equal terms in the dependant and independent variables of the MOST relations (also known as selfcorrelation). Second, we determine the impact of the uncertainty of the MOST functions on the fluxes. For any issue that we encounter regarding the violation of MOST, we discuss the extent to which the fluxes are affected.
The paper is organized as follows: in Sect. 2 we start with the theoretical background of MOST and an overview of literature functions. In Sect. 3 we describe the data processing and the routines used to fit MOST functions to the data, and in Sect. 4 we first present the overall functions of C T 2 , C q 2 and ε based on data from all combined field experiments. Second, to investigate the variability of MOST functions between datasets we present MOST functions for all datasets separately. Finally, we investigate the influence of self-correlation and the sensitivity of the fluxes H, L v E and τ to variations in MOST functions. Conclusions are given in Sect. 5. Monin and Obukhov (1954) discussed surface-layer similarity, which was set out by Obukhov (1946;English translation in 1971) to express derivatives and statistics of atmospheric surface-layer variables as functions of a dimensionless parameter: z/L, in which L, the Obukhov length, is defined as,

MOST Framework
with the virtual potential temperature θ v , the friction velocity u * = −(u w ) 1/2 , the kinematic heat flux w θ v , the acceleration due to gravity g (9.81 m s −1 ) and the von Kármán constant k = 0.4. The application of surface-layer similarity to structure parameters was first introduced by Obukhov (1960) and applied by . To define the framework of Monin-Obukhov similarity theory, the assumptions of horizontal homogeneity and stationarity were made, as well as the assumption that MOST only applies in the surface layer. Hence, it is assumed that turbulence is solely generated by shear and buoyancy, and is not influenced by non-local conditions such as extra sources of heat or moisture, or no collocation of the sources. The MOST functions that relate structure parameters and ε to the surface fluxes read, where C X 2 is the structure parameter of a conserved scalar X . In this study we focus on scalars θ and q, where θ is potential temperature, and when X is used for scalars θ and q, then X * is θ * = −(w θ )/u * and q * = (w q )/u * respectively. Definitions of C X 2 and ε are given below. Note that z is the height above the zero-plane displacement, calculated as 2/3 of the vegetation height for a closed canopy.
In Sect. 1 we introduced the assumption of MOST that heat and moisture are transported by the same mechanism and therefore in principle should have the same MOST function. However, there are two reasons for dissimilar transport mechanisms between heat and moisture. The first is due to a dissimilar transport efficiency, expressed by different Prandtl and Schmidt numbers, in idealized conditions. The second is due to the inclusion of extra sources of heat or moisture due to horizontal or vertical advection (entrainment) or no collocation of the sources that violate MOST (De Bruin et al. 1999). For example, for horizontal advection, the vegetation-air interaction affects moisture, but not the temperature of the air (Katul et al. 2008). Under these conditions the fluctuations of T and q may not be correlated, and thus C q 2 may be more vulnerable to non-local effects than C T 2 . Moreover, under very stable conditions, turbulence may cease, and the flow is governed by intermittent turbulence, drainage flows along slopes, and pressure perturbations (Mahrt 1999). In Sect. 2.4 past work is discussed that deals with these non-local influences that disturb the heat-moisture similarity. We attempt to filter out the non-local effects as much as possible from our data and attribute differences between MOST relations of C T 2 and C q 2 to differences in the Prandtl and Schmidt numbers.

Definition of C X 2 and ε
With the application of the Taylor frozen turbulence hypothesis (Kaimal and Finnigan 1994) C X 2 and ε are related to the one-dimensional turbulence spectra of X and u (S X ( f ) and S u ( f )) as follows, where the Kolmogorov constant α = 0.55, U is the horizontal wind speed and the natural frequency f is related to the wavenumber k by k = 2π f /U . Note that the relation between C X 2 , ε and the spectra is only valid in the inertial subrange where the Kolmogorov 5/3 power law applies (Kolmogorov 1941).

MOST Functions C X 2 and ε Based on Flux-Profile Relations
Here we describe how MOST functions of C X 2 and ε are related to flux-profile relations through the simplified budget equations of turbulent kinetic energy (TKE) and the variance of a scalar X . These budget equations describe the processes that suppress and enhance

123
Similarity relations for C 2 T , C 2 q and ε based on eleven field 505 turbulence. Under the assumptions of horizontal homogeneity, steady state, and neglection of transport terms, the budget equations read (Kaimal and Finnigan 1994), where ε X is the molecular dissipation of scalar variances. Next, the dimensionless groups of the mean gradient of X and U are defined as, Substituting the flux-profile relations (Eqs. 8 and 9) for the gradients in Eqs. 6 and 7, and the definitions of the dimensionless scales u * , θ * and q * for the fluxes, leads to Next, we introduce relations between ε and structure parameters C X 2 to be able to write a similarity relationship with structure parameters (Monin and Yaglom 1975;Kaimal and Finnigan 1994), in which β is the Obukhov-Corrsin constant. Substituting Eq. 12 with Eqs. 10 and 11 leads to the following MOST function for structure parameters, Furthermore, we can rewrite Eq. 10 to a MOST function for ε, Equations 13 and 14 describe how the MOST functions of C X 2 and ε are related to the MOST flux-profile relations and thus provide a means of checking the consistency of MOST relations between these parameters. In the following we use Eqs. 13 and 14 to predict the shape of f C X 2 and f ε based on well-established flux-profile relations.
Flux-profile relations are here based on the functions of Businger et al. (1971), which were modified by Högström (1988) to account for the choice of k = 0.40 instead of k = 0.35. These are, 123 for unstable conditions (z/L < 0), and, for stable conditions (z/L > 0). Inserting these relations into Eqs. 13 and 14 leads to MOST functions that we call the flux-profile-derived MOST functions, which have a neutral limit equal to 6 for f C X 2 and 1 for f ε . Note that the neutral limit of the flux-profile-derived MOST functions depends on the choice of k. Flux-profile relations have been derived with different values for k, which were summarized and modified by Högström (1988) to k = 0.4. Högström (1996) showed in his review that values of k are mostly in the range between 0.39 and 0.41 and concluded that the average value is indeed k = 0.40 ± 0.01. Almost all functions φ X listed by Högström (1988) lead to a neutral limit of 6 with k = 0.4. For φ U , all flux-profile relations listed by Högström (1988) lead to a neutral limit of 1. The neutral limits of f C X 2 and f ε are also dependent on the value of the Obukhov-Corrsin constant. Hill (1997) concluded that β = 0.86 (based on variations between 0.82 and 1), which is the value used here.

MOST Functions of f C X 2 and f ε in the Literature
Different field experiments are described in the literature, which has resulted in different formulations for f C X 2 and f ε . In Online Appendix we give an overview of earlier reported MOST functions. Several formulations of f C X 2 follow the base function proposed by , which consists of two coefficients c 1 and c 2 . In this formulation, c 1 indicates the neutral limit and 1/c 2 the inflection point of the function on a semi-logarithmic scale. For f ε , the most common base function is that proposed by Wyngaard and Coté (1971).
In the online Appendix we specify the neutral limit of the functions and for f C X 2 under unstable conditions we give an additional term c 1 c −2/3 2 when the base function of  was used; c 1 c −2/3 2 is a measure of the shape of the function under free convective conditions (De Bruin et al. 1995). Note that the values for c 1 and c 2 are interconnected; if the neutral range yields a relatively high c 1 , then c 2 should be high as well to conserve the scaling coefficient for free convection. The parameters to compare are therefore c 1 for neutral conditions and c 1 c −2/3 2 for free convection. Note that base functions other than those of  and Wyngaard and Coté (1971) have been used in other studies. Besides that, different MOST function coefficients were found when the same base function was used. For example, the base function of  was applied to other datasets (Hill et al. 1992;De Bruin et al. 1993;Li et al. 2012;Maronga 2014;Braam et al. 2014), which resulted in coefficients that deviate from the coefficients used by . For f C T 2 , c 1 ranges from 4.4 to 8.1, and c 1 c −2/3 2 from 0.93 to 1.58. Overall, the neutral limit of f C T 2 ranges from 4.3 (Kanda et al. 2002) to 8.1 (Hill et al. 1992). In the field of scintillometry, the function of  with a neutral limit of 4.9 is most frequently used to determine H , whereas it differs from the neutral limit of 6 that we derived in the previous section. The neutral limit found by Maronga (2014), which is equal to 6.1, is closest to the neutral limit of 6 that we derived. In the studies from Kohsiek (1982) and Roth et al. (2006) it is shown that the data do not level down to a neutral limit but rather show a further increase towards neutral conditions. We reason that this is because H and θ * go towards zero at the transition from unstable to stable conditions, which makes the dimensionless group go towards infinity. Braam et al. (2014) demonstrated the Similarity relations for C 2 T , C 2 q and ε based on eleven field 507 influence of different stability ranges and regression techniques on MOST functions of f C T 2 under unstable conditions, which explains part of the variation in reported MOST functions. Moreover, they concluded with measurements at various heights between 2.3 m and 58.1 m that MOST functions are height-dependent. For f C q 2 , neutral limits of 3.5 ± 0.1 and 6.3 were found by Li et al. (2012) and Maronga (2014) respectively. As was discussed in Sect. 2.1, MOST functions of C T 2 and C q 2 may differ under the influence of non-local effects or due to different Prandtl and Schmidt numbers. A problem with non-local terms is that they cannot be quantified from typically available turbulence measurements. These data are available through LES models and turbulence can be simulated under ideal and non-ideal circumstances with such a model, which provides a means to study the effect of MOST violation on MOST functions. Although LES modelling results are restricted to near-free convection conditions, Maronga (2014) demonstrated in his LES study that when entrainment was simulated, f C q 2 diverges from f C T 2 . Maronga et al. (2014) quantified the effects of surface heterogeneity on C T 2 and C q 2 and found that an heterogeneous surface generates additional temperature and moisture fluctuations leading to higher C T 2 and C q 2 values. This effect did modify the MOST functions but the effect was very small. Li et al. (2012) and Van de Boer et al. (2014) made an effort to parametrize non-local effects using local turbulence parameters, i.e., the correlation coefficient R T q and the skewness of q. Li et al. (2012) showed that from weakly unstable to neutral conditions, R T q decreases progressively from near one to near zero and derived a dissimilarity ratio that scales with R T q . Van de Boer et al. (2014) observed dissimilarity for MOST functions of the variance of humidity that they were able to relate to entrainment using the skewness of q. They suggested use of an extra term in the flux-variance MOST functions to capture the effect of entrainment, where this term includes the entrainment ratio for humidity and the boundary-layer height. The fact that similarity between f C T 2 and f C q 2 may not always be valid confirms that we should consider functions for f C q 2 and f C T 2 separately. Note, however, that our study does not investigate the dissimilarity between f C T 2 and f C q 2 but rather the differences between MOST functions of different datasets.
For f ε , Frenzen and Vogel (2001), Pahlow et al. (2001) and Hartogensis and De Bruin (2005) reported neutral limits below 1, which implies that the production by shear and the dissipation of TKE are not in balance. Hartogensis and De Bruin (2005) further discuss the neutral limits that were reported before. It is clear from the studies that do not keep the neutral limit fixed at 1 that the data do not fulfil the expectation that the neutral limit it equal to 1. Instead, the deviations from the ideal behaviour are the rule rather than the exception.

Methodology
Essential to our approach is that we use many datasets that are all processed with the same algorithms. In the following we describe the data processing and we discuss how we fit MOST functions to the data.

Data Description
We use data from 11 field experiments carried out by the Meteorology and Air Quality group of Wageningen University. An overview of the field experiments is given in Table 1. The datasets are gathered over different surface and vegetation types, and the data cover different stability ranges. Table 1 includes Bowen ratio classifications for the measurement location 123 and surrounding areas. Large differences between these classifications are an indication that a dataset is more sensitive to mesoscale circulations induced by thermal and moisture heterogeneities and may therefore show a larger dissimilarity between f C T 2 and f C q 2 . All measurements were made relatively close to the ground or above the vegetation, mostly at around 2 to 4 m in height. The measurement and vegetation heights of the experiments are given in Table 1, and for a few experiments we used changing vegetation heights as the vegetation grew significantly over the measurement period. At all field experiments, data were obtained with the same eddy-covariance (EC) system that consisted of a sonic anemometer (Campbell Scientific CSAT3) paired with a LI-COR 7500 hygrometer. The EC systems operated at 20 Hz for all field experiments.
We used the flux-software package EddyPro (v5.1.1) from LI-COR Biosciences Inc. (Lincoln, Nebraska, USA) to calculate turbulent fluxes. All standard data treatment and fluxcorrection procedures were included, most notably axis rotation with the planar-fit procedure (Wilczak et al. 2001), raw data screening (Vickers and Mahrt 1997) and low-pass filtering effect (Moncrieff et al. 1997). In addition, error estimates of the fluxes were determined within EddyPro based on Finkelstein and Sims (2001), and quality flags were determined based on Mauder and Foken (2004).
We calculate C X 2 and ε from the power spectra of T, q and u for every 30-min interval using Eqs. 4 and 5. The time series of T, q and u used in this procedure were taken from EddyPro after their processing level 6, i.e., after data pass all screening and correction procedures described above that are applicable to the raw data series. In this way we ensured that the data used for flux calculations and for the turbulence parameters are identical.
We use the MATLAB algorithms as described by Hartogensis (2006), firstly to determine S X ( f ) and S u ( f ) from measured time series using the fast Fourier tranform, and secondly to obtain C X 2 and ε. It is essential in this procedure to distinguish the inertial range of the spectra as Eqs. 4 and 5 are only valid in the inertial range.
We select C X 2 and ε such that, -at least 15 % of the total spectrum has an inertial range behaviour, which is fulfilled when the calculated r.m.s. of C X 2 and ε is within 20 % of the average value for that part of the spectrum, and when the slope of the spectrum is within 10 % of the theoretical −5/3 slope (Hartogensis and De Bruin 2005).
-C T 2 , C q 2 and ε, appropriately scaled with z for neutral conditions, i.e., z 2/3 for C T 2 and C q 2 , and z for ε, have a maximum value of 10 −1.5 , 10 −9 and 10 −4 respectively, whereas the dimensionless groups of C T 2 , C q 2 and ε, are set to a maximum value of 25, 25 and 10 respectively. Furthermore, we only include flux data that meet the following criteria: -the number of samples used for flux calculations should at least be 99 % of the total number of samples expected over a 30-min interval with 20 Hz. -the flux quality control flag based on Mauder and Foken (2004) should equal zero (the highest possible quality class). This flagging system follows the agreement made during the CarboEurope-IP workshop in 2004 (Mauder et al. 2013). -data for unstable conditions can only include daytime data, whereas stable conditions can occur both during the day and night. -heat fluxes H close to zero (−2 < H [W m −2 ] < 5) are discarded.
-data points with excessive large or small error weights were removed from the dataset.
This was done using an initial regression and fitting a Gumbel distribution to the weights of the data points. Next, a power density function of the weights was determined based 123 on 100 bins, and data points in bins that represented less than 0.25 % of all data points were removed.
The choice of data filtering aids in minimizing the influence of non-local effects such that MOST assumptions are met. The strict inertial range test, for example, removes spectra that are distorted due to non-local effects and intermittency that are likely to show up as additional spectral energy at relatively low frequencies. Also important here is that the datasets that we used were gathered relatively close to the ground, this means that the data show a short inertial range that is easily distorted by non-local effects. Furthermore, the daytime stable conditions that are filtered out are indicative of advection. We tried to relate remaining outliers to nonlocal effects through R T q and the skewness of q, but this turned out to be non-conclusive, i.e., there was no clear relation between the data points and these local parameters. We did relate outliers in the dimensionless groups to larger errors in the underlying estimates of C X 2 , ε and fluxes. Therefore we see the need to apply error-weighted curve fitting, which we further discuss below.

Fitting MOST Functions to Data
The base functions for C X 2 that we use to fit a relationship to the data are, for unstable and stable conditions respectively. These functions are the general functions proposed by . Note that other base functions have been proposed as well, however, we use Eqs. 19 and 20 as they are most commonly used and they fulfil the general characteristics of a MOST function that are given by theory (a neutral limit and scaling to (z/L) −2/3 under free convection). For MOST functions of ε we use, for unstable and stable conditions respectively. The function for f ε,unstable was proposed by Högström (1990), and we use this base function because it is consistent with the MOST function that we derived using the flux-profile relations. For f ε,stable we use the base function proposed by Wyngaard and Coté (1971), which assumes a fixed value for the neutral limit, i.e., c 1 = 1. However, in previous studies it was found that f ε in near-neutral conditions can deviate from 1, and so we do not keep c 1 fixed but determine this coefficient from the data. In determining the relation between the MOST dimensionless groups we use an errorweighted, semi-logarithmic, orthogonal distance regression. We do this procedure many times on sub-samples of the dataset that allows us to estimate the uncertainty of the MOST parameters. In the following we explain the regression procedure further. Braam et al. (2014) discuss various regression approaches for MOST functions and conclude that an error-weighted, logarithmic, orthogonal distance regression is the optimal. We follow their approach but take a semi-logarithmic scale as the dimensionless groups do not have a clear log-normal distribution. The conversion of Eqs. 19-22 for application on a semi-logarithmic scale is given in Appendix 1.
Similarity relations for C 2 T , C 2 q and ε based on eleven field 511 The orthogonal distance algorithm, also known as total least-squares regression, is after Petrás and Bednárová (2010). Flux errors determined by EddyPro and estimated errors of turbulent and non-turbulent parameters were combined to obtain uncertainty estimates of the MOST dimensionless groups through standard numerical error propagation. With these uncertainty estimates we determined the so-called internal error (Err int ) of the regression, where Err Y is the error of the dimensionless group of C X 2 or ε, and Err X is the error of z/L. The weight given to each data-pair in the regression is taken as one over the internal error.
The overall MOST functions for all datasets combined are based on 1000 sub-sample regressions consisting of randomly selected sub-sets that contain 20 % of the data points. The functions of individual datasets are based on 200 sub-sample regressions consisting of 30 % of the data points. To ensure that for each sub-sample regression a representative z/L range is included, we divided the z/L range into 10 bins of equal size and sampled an equal amount of samples from each bin. From the distribution of the MOST parameters c 1 and c 2 , thus obtained, we determined a centre-valued c 1 -c 2 pair that represents the overall MOST function. Additionally, we determined a minimum and maximum c 1 and c 2 estimate that represents the uncertainty range of the MOST parameters. The details of the procedure that is followed to determine c 1 and c 2 and their uncertainty estimates are given in Appendix 2.

Results and Discussion
In Sect. 4.1 to 4.3 we present overall functions f C T 2 , f C q 2 and f ε , together with the uncertainty of these functions, based on data from all combined field experiments. The variability of MOST functions from individual datasets is discussed in Sect. 4.4. In Sect. 4.5 we test the influence of self-correlation on MOST functions and in Sect. 4.6 we show the sensitivity of fluxes to varying MOST functions.

2
In Fig. 1 we present data from all combined field experiments for which we fitted one overall MOST function for C T 2 under unstable (top) and stable conditions (bottom). The left plots show the data with their weights colour-coded (the higher the uncertainty, the lighter the colour) and the fitted MOST functions. The right plots show the same data but with the functions that were previously reported in the literature. The right plots also include the MOST functions that we derived in Sect. 2.3 based on flux-profile relations. It must be noted that due to the linear y-axis all functions seem to have the same free convection limit. To examine the difference between MOST functions under free convection conditions we discuss the differences in the c 1 c −2/3 2 parameter. The corresponding coefficients c 1 , c 2 and c 1 c −2/3 2 that we find for f C T 2 , f C q 2 and f ε under both unstable and stable conditions are given in Table 2. In the online Appendix we give the abbreviations of the references that we use in Figs. 1, 2 and 3. The MOST functions from the literature are plotted over the z/L range as was specified in their data descriptions. If not specified in the literature, we plot the function over the full z/L range.   Table 2 Coefficients of MOST functions f C T 2 , f C q 2 and f ε including the uncertainty based on the sub-sample regressions as described in Appendix 2. Note that, following Eq. 22, the neutral limit of f ε for stable conditions is equal to c 1 Similarity relations for C 2 T , C 2 q and ε based on eleven field 513 In Fig. 1 we observe that the dimensionless group of C T 2 shows considerable scatter near the neutral limit (−z/L < 10 −1 ). The reason for this is that there are two types of conditions that can occur in this range: 1. True neutral conditions where, due to strong wind, mechanically-generated turbulence is dominant over buoyancy-generated turbulence such that z/L is small but with θ * sufficiently large. These data points tend to have the highest weight in the regression in this range. 2. Day-night transitions where, independent of wind speed, z/L is small since the heat flux changes from positive to negative and θ * has a value near zero. With θ * close to zero, the dimensionless group of C T 2 goes towards infinity. We applied a data filter to avoid the dominance of near-neutral singularity in the dimensionless group of C T 2 . The remaining data that show this singularity tend to have relatively large errors and bare little weight in the overall regression result.
It should be noted that by using the C X 2 base function that was chosen, the heat fluxes close to the day-night transition are slightly overestimated. However, in absolute terms this effect is only small, which we further discuss in Sect. 4.6.
The neutral limit that we find for f C T 2 under unstable conditions (5.6) is within the range of neutral limits that were given in the literature. Neutral limits of f C T 2 from De Bruin et al. (1993) are well below the neutral limit that we find. Furthermore, the neutral limit of f C T 2 and the free convection parameter c 1 c −2/3 2 (1.60) are both mostly consistent with that of the function proposed by Maronga (2014), who determined f C T 2 based on simulated turbulence from a LES model. Note that the function from Maronga (2014) is partially hidden by the flux-profile-derived function. The fact that the MOST functions are very similar to those found by Maronga et al. (2014) based on idealized turbulence in a LES model strengthens our case that the selected data indeed follow MOST.
The neutral limit of f C T 2 for stable conditions (5.5) is consistent with that from unstable conditions (5.6). The function that we find for stable conditions is closest to that of Hartogensis and De Bruin (2005) and Li et al. (2012). Several other functions from the literature lie higher than our function. A potential reason for this is that these functions were determined without an error-weighted regression. Larger scatter is found for very stable conditions and the variation between MOST functions from earlier studies is larger as well. It is well known that MOST breaks down under very stable conditions, where fluxes are near zero and mechanisms other than turbulence dominate the flow (Mahrt 1999;Bou-Zeid et al. 2010). The likely violation of MOST under very stable conditions makes that the correct shape of the MOST function for a particular situation is uncertain. For practical applications, when the MOST functions are used to estimate fluxes, the impact of using a non-exact MOST function for very stable conditions is relatively small in absolute terms as the fluxes are close to zero. Towards neutral conditions turbulence is well developed and fluxes are larger, and this is also where the different functions from the literature converge, giving the MOST function a lower uncertainty. In Sect. 4.6 we discuss how the uncertainty of MOST functions affect the fluxes.
From comparison of f C T 2 that we found with the one that follows from the flux-profile relations, we see that the neutral limit of the flux-profile-derived function (6.0) is within the uncertainty range that we determined for the neutral limit. Moreover, the overall shape of the two functions is also strikingly similar towards more unstable conditions. This indicates that f C T 2 , φ X and φ U are consistent, which is a strong experimental evidence for the validity of the MOST framework.
For stable conditions f C T 2 differs more from the flux-profile-derived function than for unstable conditions, which is in line with Hartogensis and De Bruin (2005). They compared 123 their results with other flux-profile relations, such as that from Beljaars and Holtslag (1991) that corrects for the effect of intermittent turbulence by enhancing the turbulent mixing for very stable conditions ( Van de Wiel et al. 2002). This flux-profile relation showed a better agreement with the corresponding f C T 2 functions. More recently, Katul et al. (2011Katul et al. ( , 2013 proposed flux-profile relations that account for anisotropy, advection and pressure fluctuations, which could in the same way give better correspondence with the MOST relations shown here. Figure 2 has the same set-up as Fig. 1 and shows data and MOST functions of C q 2 ; f C q 2 includes more near-neutral values than f C T 2 because of the stricter near-neutral filtering we adopted for C T 2 to avoid the dominance of the day-night transition singularity of the C T 2 dimensionless group. In Sect. 2.1 we discussed why C q 2 is more vulnerable to non-local effects than C T 2 and that filtering out these effects is more critical for C q 2 . To ensure that

123
Similarity relations for C 2 T , C 2 q and ε based on eleven field 515 non-local effects are brought to a minimum for C q 2 we excluded two datasets from the overall f C q 2 fit, namely Omseson and SudMed. These datasets were gathered over irrigated fields in a desert environment, and are therefore vulnerable to advection induced by thermal and moisture heterogeneities. The Bowen ratio classifications in Table 1 indicate that these two datasets indeed stand out by their heterogeneous surface.
In Fig. 2 we observe that, in contrast to C T 2 , for unstable (near-neutral) conditions the dimensionless group of C q 2 does not go to infinity. We attribute this difference to the fact that L v E mostly remains positive under stable conditions and does not go through zero as does H . Where the term θ * goes to zero during transitions from stable to unstable, the term q * remains large and the dimensionless group does not go to infinity. For this reason f C q 2 has a well-defined neutral limit, which is lower than that of f C T 2 .
We emphasize that by removing the Omseson and SudMed datasets, the amount of scatter has been reduced substantially, thereby decreasing the uncertainty of f C q 2 . To the contrary, removing these datasets for the overall C T 2 function did not change the amount of scatter in the data points of f C T 2 (not shown here). This indicates that C q 2 is indeed more vulnerable to non-local effects than C T 2 . Our result are in line with other reported results; the fact that most of the scatter originated from the datasets with the largest likelihood of being influenced by non-local effects, indicates that non-local influences can indeed affect the shape of f C q 2 . Li et al. (2012) showed that from weakly unstable to neutral conditions, R T q decreases progressively from near one to near zero, reasoning that this was related to non-local effects and derived a dissimilarity ratio that scales with R T q . Maronga (2014) showed with his model that when entrainment was sufficiently small, the MOST functions of f C T 2 and f C q 2 were identical. However, when entrainment was included he found dissimilarity between f C T 2 and f C q 2 . We tried to parametrize non-local effects with R T q , the skewness of q and also the Bowen ratio, but we could not relate any of these parameters to the dissimilarity between C T 2 and C q 2 . We reason that we did not find a relation because we already minimized the influence of non-local effects with our data treatment. Following this argument we cannot ascribe the dissimilarity between f C q 2 and f C T 2 to non-local effects per se. What remains is true dissimilarity, i.e., non-equal turbulent Prandtl and Schmidt numbers. Maronga (2014) and Li et al. (2012) found that the free convection parameter c 1 c −2/3 2 is equal to 1.66 and 1.28 ± 0.14 respectively, where we find 1.2. Just as for f C T 2 we observe that the neutral limit of f C q 2 for stable conditions (4.5) is consistent with that from unstable conditions (4.5). Also for f C q 2 under stable conditions the flux-profile-derived function deviates from the function found herein. This difference may again be related to the fact that intermittent turbulence (Beljaars and Holtslag 1991), or anisotropy, advection and pressure fluctuations (Katul et al. 2011(Katul et al. , 2013 are missing mechanisms in the flux-profile relations as we described in Sect. 4.1.

Overall MOST Functions for ε
In Fig. 3 we present our findings for f ε in a similar way as in Fig. 1 for f C T 2 . We find a neutral limit of f ε smaller than 1 for both stable (0.83) and unstable (0.88) conditions. Note that the neutral limit of f ε under stable conditions is not equal to c 1 , but to c 1 3/2 (Eq. 22). Even taking into account the uncertainty range of the neutral limit f ε does not reach 1, which implies that there is no balance between TKE production by shear and dissipation. These findings are in line with previous studies where neutral limits varied from 0.38 (Roth et al. 2006) to 1.24 (Högström 1990). Hartogensis and De Bruin (2005) discussed the neutral limits that were reported in the previous literature in more detail. The neutral limits that we find are mostly   Frenzen and Vogel (2001), and 0.8 by Hartogensis and De Bruin (2005) and also towards very stable conditions the function found here is closest to the function of Hartogensis and De Bruin (2005). Outside of the neutral range our functions match well with the flux-profile-derived functions for both stable and unstable conditions.

Variability of MOST Functions Between Datasets
MOST functions for C T 2 , C q 2 and ε that were fitted using data from individual datasets are presented in Fig. 4 and the corresponding coefficients are given in Online Appendix. Some datasets only consisted of a limited number of data, and those with less than 90 data points were therefore removed.
The top left plot of Fig. 4 shows that the neutral limit of f C T 2 under unstable conditions ranges between 4.5 (SudMed-less dense) and 7.6 (LITFASS'03-colza). Several datasets have a neutral limit close to 4.9, as found by  and De Bruin et al.
Similarity relations for C 2 T , C 2 q and ε based on eleven field 517  MOST functions for C T 2 (top), C q 2 (middle) and ε (bottom) for every dataset separately for unstable (left) and stable (right) conditions. Functions are only plotted in their line-style over the stability range that was covered by the data. Coefficients c 1 , c 2 , c 1 c −2/3 2 , the number of data points and the z/L range of all fits are given in the online Appendix 123 (1993). Besides variations in the neutral limit we find variations under convective conditions as well; c 1 c −2/3 2 ranges from 1.11 (SudMed-less dense) to 3.65 (Omseson). The neutral limit of f C q 2 under unstable conditions ranges between 2.3 (Omseson) and 4.7 (Transregio'08-sugarbeet). All datasets have a lower neutral limit for f C q 2 compared to f C T 2 , which we reasoned to be an effect of the dimensionless group of C q 2 that does not go to infinity, see Sect. 4.2. Moreover, comparison of the point cloud of f C q 2 and f C T 2 (unstable conditions) in Fig. 4 shows that the scatter is larger for f C q 2 , which implies that C q 2 is more sensitive to non-local effects than C T 2 , as discussed in Sect. 2.1. For stable conditions, for both f C T 2 (upper right plot) and f C q 2 (middle right plot), the largest differences exist for very stable conditions. Note that variability towards very stable conditions causes only small variations in the absolute values of the fluxes, as the fluxes are small under very stable conditions. The neutral limits for MOST functions of ε vary from 0.60 (SudMed-less dense and Transregio'08-sugarbeet) to 1.28 (BLLAST-grass) taking together the unstable and stable functions. Overall, MOST functions from all datasets imply that production by shear and dissipation of TKE are not in balance, as has been observed in previous studies.
Furthermore, we aim to identify plausible causes for variations of MOST functions in the literature. The potential causes that we are able to investigate here are: (1) the accuracy of z measurements, (2) a limited z/L range covered by the data, (3) a limited number of data points, (4) environmental influences such as surface type and humidity, and (5) non-local effects. Other factors that may affect the data are variable instrumentation and processing techniques, where the latter was covered by Braam et al. (2014), but which we are not able to identify in this study. For the analysis below we use data from the online Appendix, which includes c 1 , c 2 and c 1 c −2/3 2 for all individual fits, the number of data points and the z/L range of the data.

Accuracy of z Measurements
We observed in the analysis that when calculations were made with inaccurate values of z, the MOST functions fall outside the range of data points from other datasets. For example, the LITFASS'12-rye dataset was gathered in a period during which the crop was in its growing stage. Calculating the dimensionless groups with the average crop height over that period resulted in more scatter than when the actual crop height was used. Part of the variation that we find in Fig. 4 may therefore be a result of inaccurate z measurements. This may especially be the case for the SudMed datasets, which were gathered over an olive grove where the olive trees were separated by open spaces. The uncertainty in the effective displacement height over these open spaces induces an uncertainty in the further analysis. This height uncertainty would explain why the neutral limit of this dataset deviates from other datasets for all functions f C T 2 , f C q 2 and f ε . The same holds for the free convection parameter c 1 c −2/3 2 of the functions f C T 2 and f C q 2 . This observation of deviating MOST functions with inaccurate z measurements is in line with Hartogensis et al. (2003) who did a sensitivity analysis of all parameters that link C T 2 to H and showed that z is indeed the most sensitive parameter.

Limited z/L Range and Limited Number of Data Points Covered by the Data
When MOST functions are determined over a short z/L range or with a limited number of data points, they have a larger uncertainty and result in an inaccurate or invalid MOST function.
Similarity relations for C 2 T , C 2 q and ε based on eleven field 519 For example, for f C T 2 the LITFASS'12-rye dataset consists of the largest amount of data points and has the smallest uncertainty range for c 1 for both stable and unstable conditions. In contrast, the Omseson dataset only covers a short z/L range and has a limited number of data points, leading to a function with higher uncertainty, which is in this case visible as a higher neutral limit (7.0) of f C T 2 for unstable conditions and a value for c 1 c −2/3 2 (3.65) far from that found for other datasets. The same applies to the BLLAST-grass dataset that has a neutral limit of 1.28 for f ε under unstable conditions, which is substantially higher than the neutral limit of other datasets. Note also that we already excluded a few datasets with less than 90 data points from Fig. 4, because we could not determine a proper function for these datasets. The argumentation about a limited number of data points is therefore not limited to the examples of datasets that we give here.

Surface Variability
We could not find a relation between the MOST functions and environmental influences such as the surface type (e.g. vegetated or bare soil) or humidity of the site, which suggests that the MOST functions of C T 2 , C q 2 and ε can be used across different sites. For f C q 2 we observed differences that are likely related to non-local effects induced by thermal and moisture heterogeneities of the surface. We discuss this below.

Non-local Effects
The Omseson and SudMed datasets were gathered over irrigated areas within very dry environments. This large thermal and moisture heterogeneity makes these datasets highly vulnerable to advection (Hoedjes et al. 2002(Hoedjes et al. , 2007. Figure 4 shows that it is also these two datasets for which the data points of f C q 2 (unstable) are below that of other datasets. In Sect. 4.2 we showed that, when the Omseson and SudMed datasets were excluded, the scatter of the point cloud for f C q 2 was reduced, but not for f C T 2 . This indicates that C T 2 is less vulnerable to non-local effects than C q 2 . The variation between functions from all individual datasets in Fig. 4 is of the same size as the variation between functions from the literature (Figs. 1, 2 and 3). We have now seen that MOST relations that are based on individual datasets of limited length can lead to biased results due to the limited number of data points, limited stability range, or non-local effects specific to that location. If the datasets that we marked as uncertain are removed from Fig. 4 (e.g. SudMed for its uncertainty in the displacement height, Omseson and BLLAST-grass for their short z/L range, Omseson for advection), this reduces the variation between the MOST functions to a smaller range than the variation between literature functions. Based on this analysis we can say that it is not sufficient to determine MOST functions from only one dataset.

Self-Correlation
When applying MOST we have to be aware of self-correlation as θ * and u * are included on both x-and y-axes for f C T 2 and f ε respectively (Hicks 1978;Baas et al. 2006). When self-correlation dominates, we expect to observe a relation between the dimensionless group of C T 2 and z/L when random values of C T  Figure 5 shows the dimensionless group of C T 2 (left) and ε (right) against z/L, both dimensionless groups are calculated with random values of C T 2 and ε. For all other variables involved we take the measured values. For C T 2 we used random values between 10 −1 and 10 0 (exponents were randomly defined between −1 and 0) and similarly we used random values for ε between 10 −3 and 10 −0.5 . Figure 5 shows that there is a relation between the dimensionless groups and z/L, which implies that MOST functions are to some extent sensitive to self-correlation. However, large scatter exists for these relations, which reduces significantly when measured values of C T 2 and ε are used, as evident in Sects. 4.1 and 4.3. The MOST functions that we determined are therefore not a result of self-correlation and can only be defined from dimensionless groups with measured C T 2 and ε.

Sensitivity of Surface Fluxes to Varying MOST Functions
Here we test the sensitivity of fluxes for varying MOST functions. In Fig. 6 we compare fluxes H MOST that are calculated from three MOST functions: (1) the overall MOST function that we found in Sect. 4.1, and (2)-(3) the MOST functions that cover the uncertainty range of the overall MOST function defined through sub-sample regressions. The coefficients of all three MOST functions used here are given in Table 2. The bar plots in Fig. 6 show the comparison between measured (EC) and calculated fluxes (MOST) for H (top), L v E (middle) and τ (bottom) for unstable conditions. We divided the measured fluxes into bins with the same number of data points. For every bin we calculated the median of the relative difference between the measured and calculated fluxes, i.e., (H MO ST − H EC )/H EC for H . For fluxes calculated with the overall MOST function, the median values are shown as bars in the left plots of Fig. 6; the error bars indicate the median when the fluxes were calculated with the uncertainty range of the overall MOST functions. The plots on the right are similar but with the absolute difference between the measured and calculated flux, i.e., H MOST − H EC for H . The bars in Fig. 6 express the adequacy of our MOST functions in capturing the relation between turbulence parameters and fluxes compared to measured fluxes. The fact that the fluxes from the eddy-covariance method and MOST functions are not exactly the same indicates that the base functions for the MOST relations have difficulty covering all the variability over the wide z/L range. From the top left plot of Fig. 6 we can see that the Similarity relations for C 2 T , C 2 q and ε based on eleven field 521 For L v E and τ the relative differences are 4 % and 9 % at an average value of 351 W m −2 and 0.18 N m −2 respectively, and absolute differences of 13 W m −2 for L v E and 0.02 N m −2 for τ . For stable conditions we find relative differences of the same order of magnitude as for unstable conditions, i.e., 4 % for both H and L v E at an average flux of −40 W m −2 for H and 48 W m −2 for L v E. The relative differences that we find here are relatively small because of the narrow uncertainty of the overall MOST functions. However, the relative difference is larger when the MOST function variability from the literature is used. Hartogensis et al. (2003) showed that the uncertainty of literature functions led to a 20 % difference in H in the free convection range. We have now brought this uncertainty down to 6 %, which is based on flux calculations with the uncertainty of MOST functions that we presented in Table 2. Although the errors that we find appear to be small, they are in the same order of magnitude as error estimates of other variables in flux calculations and should therefore not be ignored. It is already generally accepted that when fluxes are calculated, they should include an error estimate of the input parameters such as C T 2 , C q 2 and ε. For future studies we now recommend adding an additional uncertainty in the parameters c 1 and c 2 . We suggest calculating fluxes with the overall functions that we presented in Table 2 and to calculate errors in these fluxes using the uncertainty range of MOST functions that were given in the same table.

Conclusions
Our main goal was to define robust MOST functions for C T 2 , C q 2 and ε, for which we used 11 field experiments carried out over various surface types and under various stability ranges. The data filtering used ensured that MOST assumptions are met. We combined all datasets to derive the MOST functions of C T 2 , C q 2 and ε, with the resulting MOST functions summarized in Table 2. The MOST function for C T 2 that we found for unstable conditions compares well with the MOST function that we derived through the Businger flux-profile relations. The fact that these functions are similar gives confidence in the overall results. A difference between the dimensionless groups of C T 2 and C q 2 is observed for near-neutral conditions. That is, the dimensionless group of C T 2 tends towards infinity for near-neutral conditions, in contrast to that of C q 2 . This difference is primarily driven by the unstable to stable transition. We observe that C q 2 , and thereby also f C q 2 , is more vulnerable to non-local effects than C T 2 . We investigated a parametrization of non-local effects in the combined dataset with R T q , the skewness of q and also the Bowen ratio, but we could not relate any of these parameters to the dissimilarity between C T 2 and C q 2 . We reason that this is because we have already minimized the influence of non-local effects with the data treatment. For ε we found a neutral limit well below 1, which implies that production by shear and dissipation of TKE are not in balance, as found in previous studies.
We also wished to identify if MOST functions are similar for different datasets, which is an indication of whether MOST functions are generally applicable across different sites and environmental conditions. As the effect of variable instrumentation and data treatment Similarity relations for C 2 T , C 2 q and ε based on eleven field 523 is eliminated, we were able to compare MOST functions for different datasets. This analysis showed that functions, determined over a short stability range, with a limited number of data points, or with inaccurate z measurements, have larger uncertainty and therefore can deviate from the overall MOST function. Hence, we conclude that these are likely reasons for variations of MOST functions in the literature in addition to variable instrumentation, and algorithms to estimate fluxes and C T 2 , C q 2 and ε. For MOST functions of C T 2 and ε we found no environmental influence; for example, field experiments over a dry or moist surface did not result in significantly different MOST functions, which implies that MOST functions of C T 2 and ε are indeed generally applicable. MOST functions of C q 2 differed between datasets, which was related to surface heterogeneity, indicating that non-local effects likely have an influence on these datasets. Based on this we can say that it is not sufficient to determine MOST functions based on only one dataset, as has been done in previous studies.
We encourage further research that focuses on indicators for MOST violation and to develop a parametrization to correct for these conditions. A potential powerful approach is that of Maronga (2014) and Maronga et al. (2014) using LES modelling. Such an approach allows the quantification of the effect of non-local conditions on MOST relations. These studies can potentially provide MOST functions that correct for, e.g., the influence of nonlocal effects, such that MOST functions of C q 2 can be applied under all conditions. Lastly, we demonstrated the sensitivity of H, L v E and τ to the uncertainty of the overall MOST functions that we determined using sub-sample regressions. The errors caused by different MOST functions are in the same order of magnitude as error estimates for other variables in the flux calculations. We therefore recommend inclusion of the uncertainty of MOST functions into the uncertainty calculations of surface fluxes. Furthermore, to reach consensus on the MOST functions to be used for flux calculations we recommend using the MOST functions and their uncertainty range that we determined and that are summarized in Table 2.
All functions now apply to a semi-logarithmic scale.
Since the results are presented with functions on a linear scale, we transform the coeffcients c 1 and c 2 from the semi-logarithmic fit  to the linear form of the functions (Eqs. 19-22).

Appendix 2: Sub-sample Regressions
We used sub-sample regressions to determine the distribution of the fit-coefficients c 1 and c 2 . These sub-sample regressions consist of 1000 fits from randomly selected sub-sets that contain 20 % of the full (filtered) dataset. All of these fits together determine a range of the fit-coefficients c 1 and c 2 from which we determine the coefficients of the overall fit. An example of the distribution of c 1 and c 2 is shown in Fig. 7 for fits of f C T 2 under unstable conditions. The value of c 1 for the overall fit was defined as the 0.5 quantile (median) of the distributions of c 1 . For c 2 we define a range of 0.01 quantile around the median of c 1 (shown with a grey background in Fig. 7); c 2 is then defined as the 0.5 quantile of the data points in the area around the median of c 1 . The reason that we define c 2 in the area around c 1 is that in this way the two coefficients are coherent. Furthermore, we defined uncertainties for the coefficients c 1 and c 2 . For c 1 the uncertainties are defined as the 0.1 and 0.9 quantile of the distribution of c 1 and the uncertainty for c 2 was then determined in an area around the 0.1 and 0.9 quantile of c 1 in the same way as c 2 for the overall fit was determined.  Fig. 7 Distribution of c 1 and c 2 for fits of f C T 2 under unstable conditions as determined from sub-sample regressions with 1000 fits from randomly selected sub-sets of the dataset. The middle of the three markers gives the values of c 1 and c 2 that we determined for the overall fit, the two outer markers indicate the uncertainty of the overall fit. Vertical lines indicate the 0.1, 0.5 and 0.9 quantile of the distribution of c 1 , the horizontal line in the magnified plot indicates the 0.5 quantile of the distribution of c 2 in the area around the 0.5 quantile of c 1 (coloured in grey)