Flux–Gradient Relationships Below 2 m Over a Flat Site in Complex Terrain

The surface–atmosphere turbulent exchange fluxes are experimentally determined using the eddy-covariance method. Their estimation using profiles of the variables of interest is a less costly alternative, although restricted to certain ranges of stability and assumed to hold for relatively flat and homogeneous terrain. It relays usually on the prescription of the roughness lengths for momentum, heat and matter, the latter two being adjustable parameters with unclear physical significance. The relations are derived with data from screen level to a few tens of metres upward. The application of these expressions using data only at one level in the surface layer implies assuming zero wind speed and the land surface temperature at their respective roughness lengths. The latter is a quantity that experimentally can only be determined radiatively with a substantial uncertainty. In this work the flux-profile relationships for momentum and sensible heat are assessed over a flat site in moderately inhomogeneous complex terrain in the southern pre-Pyrenees, using data between 2 m and the surface. The main findings are that (i) the classical expressions hold in the daytime for most of the dataset, (ii) the iterative estimations using the Obukhov length and the direct ones using the bulk Richardson number provide very similar results, (iii) using a second observation of temperature avoids a radiometric measure of land surface temperature and the prescription of a thermal roughness length value, (iv) the estimations over wet terrain with high irradiance depart largely from observations.


Introduction
The exchange of momentum, energy and matter between the Earth's surface and the atmosphere takes place initially at the interface between both media (Garratt 1994). The amounts exchanged by unit of time and unit of area are the surface fluxes. The adequate knowledge of these transfers are key for a comprehensive description of the climate system, as they are the effective link between its different subsystems. Radiative fluxes can be measured with a good accuracy, while conduction rules essentially in the soil. These two mechanisms are adequately represented in numerical models, as they have well-established ruling equations, and they are mainly sensitive to some quantities, either physiographic parameters like soil texture or vegetation cover, or model variables such as land-surface temperature (LST) or soil moisture.
The sensible and latent heat fluxes in the atmospheric surface layer are mostly due to turbulent transport of thermal or mechanical origin. The only accepted theory of turbulence is of application to homogeneous and stationary conditions (Kolmogorov 1942), which is usually not the case near the ground, where surface heterogeneity is ubiquitous over land and eddies display large anisotropy due to the presence of the surface. Therefore, no wellestablished prognostic equations for the turbulence exist near the wall and empirical equations are prescribed, usually imposing logarithmic-type profiles of the variables in the first metres above the ground, an approach called MOST (Monin-Obukhov similarity theory, Monin and Obukhov 1954).
These formulations can be obtained by relating the observed fluxes directly obtained through eddy-covariance (EC) systems to the vertical gradients of the atmospheric variables, providing the so-called flux-profile relationships (Foken 2008). Expressing the fluxes in terms of vertical gradients allows a straightforward use of the expressions in numerical models and many other applications (Louis 1979;Skamarock et al. 2019). These expressions depend on thermal stratification and have been derived using experimental data during the last third of the twentieth century, following the pioneering work of Businger et al. (1971) and Dyer (1974), preluded by a decade of research (Swinbank 1964;Swinbank and Dyer 1967). Reviews of the first experiments and MOST can be found in Hicks and Baldocchi (2020) and references therein. Due to the diversity of measurements analyzed up to the present, new functions and coefficients are still proposed (Sakagami et al. 2020).
Many of the original experiments were derived in relatively flat and homogeneous terrain using a number of instrumented levels in the first tens of metres above the surface. The obtained expressions shared a similar functional shape and had different coefficients. Högström (1988), using some of those experiments, provided compact expressions, currently widely accepted, that showed a dispersion close to 20% when compared to the different functions proposed in the literature that he employed. The discrepancies can be related to the uncertainty of the EC measurements (Mauder et al. 2013) and to the use of fluxes obtained at different heights, under the basic assumption that the exchange fluxes are independent of height close to the surface, which is often not the case (Haugen et al. 1971;Fazu and Schwerdtfeger 1989).
The dependency on thermal stratification in the flux-profile relationships is made through a stability parameter, either z/L where L is the Obukhov length or the Richardson number Ri. Both L and Ri express a ratio between the thermal and mechanical source terms in the turbulence kinetic energy equation. The resulting expressions are commonly used in parameterizations of the surface fluxes in numerical models (Cuxart et al. 2006). The approach is still left with many challenges (Businger 2005), such as whether there is a constant relationship between friction velocity and wind speed with height in the surface layer (Sun et al. 2020).
Surface fluxes can be estimated experimentally with profile measurements in different ways. The classical approach, more than 50 years old, is to derive flux-profile relationships through the use of measurements at a number of levels, usually installed on a tall tower (Bosveld et al. 2020). These relationships have been used also with only two levels of measurements at varying heights depending on the set-up (Berkowicz and Prahm 1982;Brotzge and Crawford 2000). The simplest configuration is to use screen-level data (around 2 m above ground level, a.g.l.) as the upper level and estimated values close to the surface. The latter implies using a prescribed value for the aerodynamic and thermal roughness lengths (z 0 and z 0T ), the respective heights with zero wind speed and LST as the corresponding temperature. It is worth noting that sensible heat flux estimations using LST imply a larger dispersion because the value of the LST has more uncertainty than temperature measurements obtained with a standard thermometer (Simó et al. 2018). Furthermore, the measurement of LST can be very sensitive to the footprint of the sensor (Burns et al. 2003), specially over heterogeneous surfaces that combine bare soil with vegetation cover (Chehbouni et al. 2001).
There are some other factors that can increase the dispersion of the results. On the one hand, when temperature sensors are housed within naturally ventilated shields, estimated heat fluxes can have errors up to 15% (Brotzge and Crawford 2000). In order to reduce them, careful consideration of the effect of direct insolation and low wind speeds must be taken into account, since lack of ventilation is an important factor (Richardson et al. 1999). On the other hand, the dependency of the roughness length on wind direction is another important factor (Beljaars and Bosveld 1997). Finally, the unclear conceptual definition of z 0T , with values that can vary several orders of magnitude (Sun and Mahrt 1995), oftentimes makes it a tuning parameter.
Flux-profile relationships are derived implicitly assuming homogeneous spatial conditions over flat terrain. Nevertheless, many experiments have been made in locations where significant surface heterogeneity was present, such as land-water discontinuities (Grachev et al. 2018) or the presence of natural and farming land (Huo et al. 2015). The local circulations that these terrain discontinuities induce may add dispersion in the flux estimations when using the classical formulations (Guo and Zuo 2017). Other perturbing factors are the presence of relevant canopies (Novick and Katul 2020), or humidity advections (Haghighi and Or 2015), especially in the case of contrasting semi-arid and well-watered terrains (Li et al. 2019).
Finally, the studies addressing similarity scaling in the surface layer in mountainous terrain are scarce (Stiperski and Rotach 2016) and suggest that MOST does not apply since, at least, one of the required conditions (i.e. surface fluxes must be constant with height) is not fulfilled (Nadeau et al. 2013;Sfyri et al. 2018;Stiperski et al. 2020). However, it must be taken into account that these experiments were made in a variety of conditions. In the case of Nadeau et al. (2013), observations over a steep slope are used, while Sfyri et al. (2018) analyze different sites distributed within an alpine valley at heights above 4 m above ground level (a.g.l.), and Stiperski et al. (2020) focus on mesoscale katabatic flows developed in the stable regime with measurements taken much higher up.
In the current work the applicability of flux-gradient relationships is explored for data obtained very close to the surface over moderately inhomogeneous complex terrain. The site is at the bottom flat area of a pre-Pyrenean valley in the Iberian Peninsula by the city of Iruñea in Navarre (Cantero et al. 2019;Santos et al. 2020). The challenge is to see to what degree expressions derived for flat homogeneous terrain can be used in other conditions. In applications such as numerical models these expressions are used in every grid point over the surface, therefore it is of interest to check their applicability for complex terrain.
This exercise is made using 30-min data for an eight-month period in 2018 for the fluxes of momentum and sensible heat with the aim to discriminate for what range of thermal stratification and surface conditions the standard relations hold. Furthermore relationships for the fluxes, gradients and stability parameters are explored for the whole range of data. Finally the sensitivity of the results to the use of LST or air temperature at a second level is explored, as the second option allows to skip the prescription of a value for the thermal roughness length.

Experimental Site
The current work uses data from the Alaitz experiment (ALEX17, Cantero et al. 2019;Santos et al. 2020) devoted to wind energy in complex terrain. The field campaign took place in the Elortz valley, located in the south-west pre-Pyrenees in Navarre, at 450 m above sea level (a.s.l.) following a south-east/north-west orientation with the river flowing westward (Fig. 1). The bottom of the valley is devoted to agriculture with fields of a typical scale of 500 m to 1 km, limited at both sides by mountain slopes at a distance of about 2 km. The southern part of the valley is defined by the Alaitz mountain range, peaking at 1285 m a.s.l., with an average elevation of 550 m above the valley floor that configures a steep slope facing to the north. The northern part is flanked by the low Tajonar ridge that is elevated about 200 m above the valley floor. The mountain range at the south is usually able to block and divert perpendicular flows, while the northern range only modifies the low-level flows.
A description of the main circulations in the valley is given in Santos et al. (2020) and Cantero et al. (2019). Here, the analysis will focus on those periods selected to test the stability functions, consisting in two time frames considered approximately stationary during day (1100-1400 UTC) and night (0000-0300 UTC) and for a defined range of stabilities (see more details in Sect. 2.3). Figure 2 contains the wind roses at 2 m a.g.l. in the centre of the valley (Fig. 1) for the period April-December 2018. The prevalence of directions from the north-west sector is caused by the generation of upvalley winds (west-north-west/northwest) and the presence of northerly flows that here are forced by the Pyrenees to take a north-north-west direction. At night, the most common directions correspond to downvalley flows (east-south-east/south-east), southern synoptic channelled winds (south-east or west- south-west) and downslope flows from the Alaitz mountain range (south-west). The largest wind speeds occurring during the presence of south-east channelled flows are probably due to a funnel effect caused by the valley structure since its cross-section increases to the west (Santos et al. 2020). Figure 3 displays the time evolution for the soil moisture (through the volumetric water content, VWC) and the latent heat flux (LE) between 1100 and 1400 UTC at the same location for the analyzed period, the latter coloured depending on the incoming shortwave radiation (R sw ). The valley floor had enough moisture in the soil to sustain natural green vegetation until the second half of July, the water content recovering together with plant regrowth at the end of October. The evolution of LE along the year is also influenced by soil moisture. Until mid July, LE experienced a wide range of values, reaching maximum quantities (between 300 and 400 W m −2 ) for those days when incoming solar radiation was also large. We mark the summer period of dryness as the one with a soil volumetric water content below 0.14 m 3 m −3 , when most crops have been harvested and the natural terrain is in arid conditions, with mostly bare soil at the bottom of the valley, bushes over the slopes and some coniferous forest areas. During this period, LE was usually below 100 W m −2 due to the lack of water in the upper soil. In November, soil moisture recovered due to the autumn rains and consequently LE increased significantly, although not reaching the high values of spring since incoming solar radiation is smaller this time of the year.

Instrumentation
During ALEX17 a network of nine surface-layer stations (SLSs) was installed to provide both atmospheric and soil measurements throughout the valley. A SLS contains observations for the wind at 2 m a.g.l., two levels of air temperature and humidity at 2 m and 0.36 m and a soil kit with the soil heat flux at a depth of 8 cm, plus the soil moisture and temperature both at 5 cm depth ( Table 1). The hygrothermometers are housed within a self-made radiation shield, the entire system providing an error of 0.3 K and 0.2 g kg −1 for air temperature and humidity, respectively ). The wind is measured by a two-dimensional (2D) sonic anemometer that shows an estimated error of 0.1 m s −1 at 5 m s −1 and an angle error of 3 • (CampbellSci. 2017).
The stations are fully autonomous as they are self-powered with solar panels. To optimize the energy use and data storage, sensors are interrogated once every 5 min. For the wind, the interrogation is performed over 30 s at 1-Hz sampling rate; while for air temperature and humidity sensors, 20 measurements are acquired over 1 min. Soil data is interrogated once every 5 min. This configuration provides 5-min time series with samples of soil measurements (temperature, humidity and heat flux), 1-min averages of air temperature and humidity and 30-s averages of wind speed and direction. Finally, the current analysis uses 30-min averages performed over the raw 5-min time series.
One of the SLS was located at the centre of the valley ( Fig. 1) near the village of Zabalegui together with a surface energy balance (SEB) station that provided turbulent fluxes through the EC technique. The SEB instrumentation ( Table 1) was composed of a three-dimensional (3D) sonic anemometer and an open-path gas analyzer at 1.7 m a.g.l., a four-component net radiometer at 1.2 m, a thermohygrometer at 1.5 m (housed with a Gill-type radiation shield), and four ground-heat-flux plates at different locations (Cantero et al. 2019), allowing to determine the fluxes of momentum, heat and water vapour experimentally. The sonic anemometer was oriented with the positive x-axis forming an azimuth angle of 100 • with respect to the geographical north. Data sampling have a frequency of 10 Hz for all sensors.
After the experiment, an offline data assessment has been undertaken generating new files with averaged intervals of 30 min. The EC data consist in 30-min averaged turbulent fluxes processed with the TK3.11 software package (Mauder and Foken 2015), that allowed to correct and perform a quality control of the dataset. High-frequency data are de-spiked through a threshold of seven standard deviations and only 30-min periods with less than a 10% of missing or bad values are used to compute the turbulence statistics after applying a block averaging. In addition, planar fit (Schotanus et al. 1983;Moore 1986;Wilczak et al. 2001) and density corrections (Webb et al. 1980) were applied, together with a cross-correlation analysis to maximize covariances. Finally, stationarity and integral tests on developed turbulence allow to identify turbulent fluxes with a quality flag, ranging from 0 to 2, the former corresponding to suitable data for fundamental research (Mauder et al. 2013). In the current analysis only EC data flagged with 0 have been used for the comparison against estimated fluxes with the gradient approach.

Strategy of Comparison and Data Filtering
With the disposition of the SLS instrumentation, estimations of the surface fluxes of momentum and sensible heat can be obtained through the flux-gradient approach using MOST. This methodology is validated at the centre of the valley ( Fig. 1), where both SLS and SEB station were operating together during the ALEX17 experiment. In order to get as close as possible to ideal conditions, a number of restrictions have been imposed to select the data that will be treated in the next sections. These requirements can be grouped in three categories:

Restrictions Due to Monin-Obukhov Similarity Theory Applicability
-Stationarity: Half-hourly averages between 1100 and 1400 UTC and between 0000 and 0300 UTC are selected for daytime and night-time computations, respectively. These periods are considered to be the ones closer to stationary conditions and avoid morning and evening transitions at any time of the year. -Limited Range of Stabilities: Since the range for which the similarity functions are supposed to hold is −2 < z/L < 1, we restrain the comparison for values of z/L in this range, computed with the EC fluxes. Cases with z/L > 1 fall well into the z-less regime and those functions become constant (Wyngaard 1973), whereas for z/L < −2 data is entering into the local-free-convection regime (Tennekes 1970). -Limitations in the Definition of Experimental Functions: The original flux-profile relationships were built upon experimental data for which wind measurements could not provide meaningful values under 1 m s −1 since they were performed by cup vanes with mechanical inertia (Businger et al. 1971;Högström 1988). Here we restrict ourselves to wind speeds of at least 1 m s −1 even if the measurements of the sonic anemometers are trustable for weaker values.

Restrictions Due to the Eddy-Covariance Technique
-Mast Perturbation on the Turbulent Flow: Due to sonic anemometer orientation, an airflow from 70 • to 130 • is perturbed by the station structure before reaching the sensor head, providing not suitable turbulent fluxes for the current comparison. A similar range of directions (± 30 • with respect to the orientation of the sonic positive x-axis) was used in Mauder and Zeeman (2018). In addition, the error analysis done for the momentum flux between MOST estimations and EC results confirms that the application of a larger cone of influence does not provide better results. -Limitations of the EC Accuracy: Following the recommendation of Mauder et al. (2020), EC sensible heat fluxes with absolute values smaller than 10 W m −2 are not considered since the uncertainty is usually larger than the value itself. -Restriction to High-Quality Data: As mentioned in a previous section, the current analysis is limited to high quality data provided by the EC system, considering only the 30-min statistics flagged with 0. Thus, cases with large temporal variability are eliminated, as a shower event or intermittent cloud passing, together with those that do not properly follow a flux-variance similarity (Foken 2008).

Restrictions Due to Instrumental Issues
-Limitations in Naturally Ventilated Radiation Shields: When passive shields are used to protect the sensors of temperature and humidity, it has been observed that sometimes ventilation is insufficient for the level closest to the surface, especially for weak wind speeds (Erell et al. 2005). Therefore, for the comparisons using two levels of temperature from a SLS, we choose to eliminate the daylight points that fall outside one standard deviation of a linear relationship between the vertical gradient of air temperature and the incoming solar radiation (Fig. 4), attributing them to an excess of radiation heating in the lowest level.
The original dataset contains 1281 cases during the daytime period (1100-1400 UTC) and 1247 at night (0000-0300 UTC). After the application of the rest of aforementioned restrictions, 58% of data points remain for the momentum flux and 51% for the sensible heat flux during the day. At night, these percentages are reduced to 36% and 14%, respectively. The restriction that has a major impact on the filtering process corresponds to the EC flag, followed during the day by the lack of ventilation of the lower radiation shield and, at night, the omission of small sensible heat fluxes (|H | < 10 W m −2 ) and the wind direction restriction.
The large reduction of cases at night firstly responds to a failure in the stationarity test of the turbulent fluxes calculated over a 30-min window. A smaller averaging period of few minutes would be more appropriate in the nocturnal regime to avoid the influence of small-scale mesoscale motions of diverse origin (Mahrt 2007). However, in the present study we apply a homogeneous average period to the entire dataset as all stability regimes are considered.

Fig. 4
Correlation between incoming shortwave radiation (R sw ) and temperature difference between 2 and 0.36 m a.g.l. Grey dots represent the entire period of study (Apr-Oct 2018) for all 30-min cases with R sw > 10 W m −2 while the subset of data between 1100 and 1400 UTC are represented in black. Black and red lines correspond to the linear fit ± one standard deviation, respectively. Dots falling outside the area limited by both red lines are discarded from the current analysis the turbulent flux of a quantity with its mean vertical gradient through a so-called universal function φ that depends on the stability regime. Such a relationship must be determined experimentally and numerous field campaigns have been used to propose functions with different forms and coefficients (Foken 2008). Wangara 1967 (Australia, Clarke et al. 1971;Hess et al. 1981) and Kansas 1968 (USA, Kaimal and Wyngaard 1990) were among the first experiments that related vertical profiles of wind and temperature with the turbulent fluxes in the surface layer. They used instrumented towers with cup anemometers and thermohygrometers at different levels, from 0.5 to 16 m, over relatively homogeneous terrain with short vegetation.
The flux-profile relationships for momentum (τ ) and sensible heat (H ) fluxes within the surface layer can be expressed as Businger (1988) where the former flux has been expressed in terms of friction velocity u * and a coordinate system aligned with the mean wind direction is taken. We use the Reynolds notation where overbar and primes distinguish respectively an averaged quantity from its fluctuating part. Therefore,ū andθ v correspond respectively to the mean longitudinal wind component and virtual potential temperature of an atmospheric flow with an air density ρ and a specific heat capacity c p , while κ = 0.4 refers to the von Kármán constant. Finally, φ m and φ h represent the universal functions that depend on the dynamic stability, generally expressed in terms of either the Obukhov length L through the non-dimensional height (z/L), or the gradient Richardson number Ri, The experimental determination of the universal functions φ has been an intense research topic in the last 50 years and a vast amount of expressions have been suggested as the best adjusted coefficients depend on the experiment location. Nevertheless, the Businger-Dyer expressions represent a set of functions widely used in the literature (Moene and van Dam 2014), although the exact values of their coefficients depend on the version used (Businger 1988). Högström (1988) reviewed all these expressions together with others and tried to homogenize them by fixing the values of external parameters such as the von Kármán constant or the Prandtl number. Even after this effort, the expressions compared for φ m differ around a ±20% for both stability regimes, while the differences for φ h range between 10 and 25% for the unstable cases, increasing to 35% for the stable regime (Högström 1996).
In the current study different versions of the Businger-Dyer relationships have been tested and the results are shown for those functions that provide the best estimated fluxes with our dataset. In the following sections, the analytical background needed to estimate the turbulent fluxes is specified for those expressions that use z/L or Ri as stability parameters.

Similarity Functions Depending on z/L
To derive the turbulent fluxes from the similarity functions expressed in terms of z/L, Moene and van Dam (2014) suggest to use the integrated form of Eqs. (1) where θ * = −H /(ρc p u * ) and z 0 is the aerodynamic roughness length. The heights z u , z θ 2 and z θ 1 indicate the level where wind is observed and the two levels where temperature are measured, respectively. The integrated version of the universal functions ( m , h ) used in the current work are the Businger expressions (Businger et al. 1971), with the coefficients re-calculated by Högström (1988). Under unstable conditions (z/L < 0), they take the form h (z/L) = 2 ln h (z/L) = −7.8z/L.
Note that the application of the equations above requires an iterative process that eventually converges into a single solution.

Stability Functions Depending on Ri
Alternatively to what is described in the previous section, the turbulent fluxes can be expressed in terms of finite differences by approximating the former Eqs. 1-2 These expressions have been derived after converting the vertical derivatives in finite differences and considering z in Eqs. 1-2 as the logarithmic mean height of the observational levels (Arya 2001), i.e., although other approaches such as the geometric height could be equally made (Garratt and Pearman 2020). The current work uses the Dyer universal functions (Dyer and Hicks 1970;Dyer 1974) but with Ri as the stability parameter instead of z/L. Following Arya (2001), taking z/L = Ri and z/L = Ri/(1 − 5Ri) for unstable (Ri < 0) and stable (Ri > 0) cases respectively, yields the following expressions The current experimental set-up, in which there are two levels of temperature and one level of wind, requires to use the bulk Richardson number: There are several approaches that relate Ri with Ri B , the most simple consisting in assuming both parameters to be equivalent (Ri = Ri B , Oke 2002). However, Byun (1990) develops a complex analytical relationship between both parameters, although it may fail for extremely unstable or free convection. The same work also provides a simplified version (Ri = z z−z 0 ln(z/z 0 )Ri B ) that Lee (1997) modifies with a multiplying factor η. In our case, the relationship used follows the same approach with η = 0.5, This relation is the one that provides the smallest discrepancies against the EC measurements, leading to similar results as for the z/L formulation described in Sect. 3.1.

The Lower Boundary Conditions
The estimation of turbulent fluxes through conventional data normally uses observations at several levels, for instance from a mast, to adjust a function from which the gradients are obtained. For the bulk approach described above only two levels are available and therefore finite differences are taken, instead of obtaining the gradients using an empirical function. As described above, when only one measurement at a single height is present, the second level is taken at the roughness length. The current dataset, with wind and temperature measurements at one and two levels respectively, requires the determination of z 0 but also allows to assess using LST and z 0T as the lower boundary condition. The aerodynamic roughness length z 0 can be determined either by using tabulated values or by adjusting a neutrally stratified profile for the location of interest. Table 2 shows the monthly values obtained for the Zabalegui station, using the turbulent momentum fluxes for events near neutrality (−0.01 < z/L < 0.01) with wind speeds above 2.5 m s −1 , refining the value in a way that it minimizes the error between estimated and observed u * . The latter step is only possible when an EC station is available. However, it has been applied to avoid a dispersion in the results attributable to an inaccurate value of z 0 , as this issue falls outside the scope of the current study.
The yearly evolution of z 0 is consistent with the observed height of the vegetation, which is tallest in spring (0.3 m), when the crops develop, lower in summer (0.10 m), when crops are senescent and harvesting takes place, and small in fall, after the fields are ploughed (0.05 m). According to Arya (2001), z 0 = 0.15 h 0 where h 0 is the vegetation depth, which would provide consistent values for our estimations of z 0 , since a vegetation depth of 30 cm yields z 0 = 0.045 m and, a depth of 5 cm, z 0 = 0.008 m.
When only one measurement of air temperature is available, it is customary to use LST at the thermal roughness length z 0T as the second level (Beljaars and Holtslag 1991), a temperature that in models is an output of the land-surface model, and in observational studies is a quantity derived from the budget of the infrared radiation near the surface. This value can be obtained experimentally at the ground surface, although it is possible to estimate it through remote sensing techniques (García- . Here, in situ LST is obtained using observations from the four-component net radiometer of the SEB station as (Simó et al. 2018) with σ = 5.67 × 10 −8 W m −2 K −4 representing the Stefan-Boltzman constant, L ms being the radiance measured by the longwave sensor pointing to the surface and L dn , the reflected downward longwave radiation coming from the sky. The surface emissivity is set to 0.97 corresponding to senescent sparse shrubs (Snyder et al. 1998), a value applied at similar sites (García-Santos et al. 2019). Measurements of LST have a typical uncertainty of 2-3 K due to radiometer sensitivity (Simó et al. 2018) and surface emissivity (Coll et al. 2005). A priori this associated error would introduce also huge uncertainties in the determination of the estimated fluxes, especially when thermal gradients are small. Furthermore this value has to be applied at the thermal roughness height z 0T , which is a quantity difficult to define conceptually.
In order to reduce the additional uncertainty given by the prescription of z 0T , in the current analysis this parameter is defined as a constant factor of z 0 , z 0T = μ z 0 . For our dataset, μ = 0.4 is found to provide the best estimates for the sensible heat flux without increasing the error in the estimation of u * . This value falls between the range provided by Garratt (1994), where μ = 0.5 is considered for smooth surfaces and μ = 0.14 for randomly distributed elements. A similar discussion could be made concerning the estimation of the latent heat flux and the prescription of surface humidity. However this topic falls beyond the scope of the current work.

Assessment of the Similarity Functions
The remaining set of half-hourly data is inspected in this section to assess if the traditional similarity approach, developed for seemingly flat and nearly homogenous conditions, is of application for the ALEX17 site, which is locally flat, mildly heterogenous and surrounded by significant topography nearby. Since the dataset extends from spring to early winter, a large variety of meteorological and surface conditions occurs, which makes the comparison more general than if the study was focused in a temporally restricted dataset of a few days or weeks.
The analysis is intended to show in which conditions the standard empirical theory holds while identifying the conditions in which data and theory depart significantly. Turbulent fluxes of momentum and sensible heat will be the subject of our study, leaving the analysis of the latent heat flux for a later stage. Let us recall here that we will inspect the behaviour of both fluxes for 3-h intervals at day and at night that are expected to be the closest to stationary in a classical diurnal cycle, which will allow the investigation of the performance of the theory in unstably and stably stratified conditions, respectively. The analysis will show that the theory holds most of the time in the daytime, except for the sensible heat flux when the soil is saturated and the incident solar radiation is high (Fig. 3).
In the direct comparison between measured and estimated values we establish two thresholds of quality. First we will indicate with a line the 20% error, because this uncertainty is of the same order as the one that the similarity equations provide, as indicated in Högström (1996). Therefore if a value has an error lower than 20% we will consider it "good". Then we indicate the 50% error as a "tolerance limit" that may be still of use, whereas values with errors larger than 50% are considered a failure of the theory for our site. Stiperski and Rotach (2016) evaluate the statistical uncertainty of a measured EC flux, considering 50% an optimal compromise between data quality and data loss in mountainous terrain as the ideal 20% is too stringent in such complex environments.
The computation of the momentum and sensible heat fluxes is made through the similarity theory equations introduced in Sect. 3, using both stability parameters (either z/L or Ri B ). The observations of wind and temperature are taken at 2 m a.g.l. The second level for the wind is obtained prescribing a zero value at z 0 for each month (Sect. 4), while the second level of temperature is taken at 0.36 m (indicated from now on as T 36 ). For the stability functions depending on z/L, an additional estimation is done using LST at z 0T instead of T 36 as the second level of temperature. All three estimations (z/L, Ri B and LST) have been performed at day and night to assess which method provides the best results.

Momentum Flux
Some statistics of the direct comparison of the measured and estimated values of the momentum flux are displayed in Table 3. The differences between the daytime (unstably stratified) For the sensible heat flux under unstable conditions, cases are separated between dry and humid cases (see Sect. 5.2), while for stable conditions, statistics are only considered when observed H is between − 10 and − 30 W m −2 . The third column specifies the methodology used for the estimation of turbulent fluxes, where z/L and Ri B refer to the stability parameter used in the similarity theory equations with two levels of air temperature at 2 and 0.36 m a.g.l. and one wind measurement at 2 m a.g.l. LST indicates that turbulent fluxes are estimated with z/L as the stability function and LST at the heat roughness height z 0T instead of air temperature at 0.36 m. The rest of columns refer to: mean error (ME), standard deviation of the difference (SDD), percentage of experimental points with respect to the total with a relative error less or equal to 20% ( p 20% ) and 50% ( p 50% ), and number of points used for the statistics (N) and the night-time (stably stratified) are small for mean error or standard deviation. However, it is noteworthy that the percentage of good estimates is close to 75% in unstable conditions, with very similar values regardless of the form of the stability function (z/L or Ri B ) or the use of LST or T 36 as the lower temperature level. Regarding the 50% tolerance limit, we see that most of the estimates (96%) are below it, which indicates that the theory is of wide applicability for our site in the daytime. The night-time statistics in terms of tolerated error are worse, even if most of the difficult cases have been eliminated. After the removal of EC low quality records, weak wind cases and very stable stratification points, only 36% of the records are left. Of the latter, 50 to 60% of the kept records have an error below 20% and the tolerance limit of 50% is reached by 80 to 92% of the records depending on the method. The best statistics correspond to the use of the functions with z/L and the lower temperature level as T 36 , while the worst results are performed when Ri B is used, since this method underestimates u * for weak winds (not shown). Therefore the performance at night is not as good as in the daytime, but this is in line with the well-known behaviour of the similarity theory at night everywhere. Figure 5a, b depicts the distribution of daytime and night-time data, respectively. The results are shown for z/L, but they are indistinguishable of those obtained using Ri B . The proportion of points with an error below 20% increases with wind speed within the u * range of 0.15 to 0.40 m s −1 . The distribution is very similar at night, with the largest departure (a) (b)

Fig. 5
Error for the friction velocity u * computed through the difference between the estimated u * (z/L) calculated with similarity equations depending on z/L and the observed u * (EC) obtained from the EC system with respect to the observed value u * (EC) for a unstably and b stably stratified conditions. Colour scale indicates wind speed and symbols correspond to different ranges for wind direction (dir). Dashed lines delimit the area with a relative error of 20% (inner cone) and 50% (outer cone) observed for large values of wind and u * , respectively above 7.5 and 0.4 m s −1 , which supposedly correspond to weakly stable conditions since the turbulent mixing is strong. These large wind cases have a direction of 135 • (triangles in Fig. 5), that correspond to the southern general flow funnelled through the valley. Santos et al. (2020) show that the normalized vertical wind profiles within the valley floor for this wind sector vary depending on the stability regime, pointing to differences in the perceived surface roughness due to a more heterogeneous surface at this part of the valley as the most plausible cause. This fact may be responsible for the higher error in the estimation of u * in our case, as seen in Beljaars and Holtslag (1991).
Although cases with wind speeds below 1 m s −1 have been filtered out, the error for points slightly above this threshold (black symbols in Fig. 5b) is larger than for the rest of cases. Weak-wind stable conditions are not well characterized by the standard theory and receive a separate consideration in the literature (Mahrt 2008). However, an upper limit for the value of wind speed that may fall within this category is difficult to determine. Finally, it is interesting to indicate that a plot of u * versus ln(z/z 0 ) for all the stable cases (not shown) yields to an almost linear distribution with a slope close to κ = 0.4, pointing to an observed behaviour close to neutrality.
The experimental points for φ m under unstably and stably stratified conditions are displayed in Fig. 6 against the logarithm of z/L, together with the Businger curve of the form presented by Högström (1988), and specified in Table 4. The points are well distributed following the similarity expressions, departing from them more as the relative error of the estimation becomes larger. As anticipated in the analysis of the statistics, there is a good agreement for the well-behaved cases in stably and unstably stratified regimes. However it is clear in this representation that the points with an error larger than 20% correspond to situations that are not in agreement with MOST basic hypotheses, most likely due to the influence of the surrounding topography.
As a preliminary conclusion, these results show that the surface momentum flux is generally well characterized by the bulk approach of similarity theory for any of the forms (a) (b) Fig. 6 Experimental data for φ m against stability (z/L) for a unstably and b stably stratified conditions. Colour scale indicates the relative error of u * and symbols discriminate between dry (circles) and humid (triangles) cases. Solid black line corresponds to the Businger curve (B) with the coefficients re-calculated by Högström (1988) and specified in Table 4 Table 4  The coefficients α and β are specified for the original formulation of Businger (B) recalculated by Högström (1988). In addition, for φ h , α and β coefficients are displayed for the best-fitting curves calculated with the unstable dry (UD) and unstable humid (UH) experimental datasets in daytime as well as for the stable cases (S) with φ h > 0 in night-time, all of them represented in Figs. 8 and 11 used (z/L, Ri B or LST). A proper characterisation of u * at screen level (i.e. at 1.5 and 5 m a.g.l.) has also been reported for the Cooperative Atmospheric-Surface Exchange Study-1999 (CASES-99) field experiment (Sun et al. 2020). In our analysis, very stable weak-wind conditions are ruled out and the remaining weak wind stable cases provide worse statistics since most probably the wind field is dominated by non-stationary mesoscale flows (Mahrt 2008). Deviations also arise at night for strong winds from the eastern sector that may be linked to strong channelled flow originated by the topography, a regime not typical of flat homogeneous conditions. Fig. 7 a Error for the sensible heat flux H under unstably stratified conditions computed through the difference between the estimated H (z/L) calculated with similarity equations depending on z/L and the observed H (EC) obtained from the EC system with respect to the observed value H (EC), b differences between the estimated values using similarity equations depending on z/L or on Ri B with respect to the observed H through the EC system. Colour scale indicates the observed latent heat flux L E and symbols discriminate between dry (circles) and humid (triangles) cases. Dashed lines delimit the area with a relative error of 20% (inner cone) and 50% (outer cone)

Sensible Heat Flux
Analyzing the data for the sensible heat flux as just made with the momentum heat flux, different conclusions are reached. Since in general larger errors are found, a slight sensitivity to the method used to compute the stability function exists and there are large discrepancies between theory and observations when the soil is moist, especially under large solar irradiation. In Fig. 7a the estimations computed using Eq. 6 with the integrated stability functions in terms of z/L are shown. We observe that there is a large amount of points for values of H under 200 W m −2 that are located outside the 50% limit. The inspection of the data shows that these points correspond mostly to large values of L E (above 250 W m −2 ) for a soil with large available humidity (soil moisture above 0.14 m 3 m −3 ). Most of these records correspond to large values of solar irradiation (Fig. 3) and to low values of the Bowen ratio (not shown).

Daytime Cases with Dry Soil
In order to assess the goodness of the stability functions for the sensible heat flux we will first concentrate on cases with relatively dry soils. These conditions are similar, at least looking at the climatological precipitation amounts, to the ones that were used by Businger and Dyer to derive their expressions, corresponding to data from, respectively, summer in Colorado (USA) and winter in a dry area of Australia.
When the soil VWC value is below 0.14 m 3 m −3 , data are considered to be under dry conditions, indicated by circles in Figs. 7, 8, 9, 10 and 11. Most of the daytime dry cases in Fig. 8 The same as in Fig. 6a but for φ h during the day. Colour scale indicates the Bowen ratio (Bo) and symbols discriminate between dry (circles) and humid (triangles) cases. Solid black line corresponds to the Businger curve (B) with the coefficients re-calculated by Högström (1988). Red and blue curves indicate the best fitting functions for unstable dry (UD) and unstable humid (UH) cases, respectively. Coefficients for all functions are specified in Table 4 (a) (b)  Fig. 7a (using the stability function with z/L) show errors below 20%, with Bowen ratios (Bo) close or larger than 1, meaning that sensible heat flux H is similar or larger than latent heat flux L E. In Fig. 7b the difference of the results using the similarity function with Ri B indicates that in general the latter provides slightly lower estimations than for z/L, although the differences are of the order of a few tens of W m −2 at most, usually below 10% of the EC value and similar to the sum of the mean error of both methods.  Högström (1988). Green curve indicates the best fitting function for night-time cases with φ h ≥ 0 (S). Coefficients for all functions are specified in Table 4 Figure 8 shows experimental φ h values against z/L in logarithmic scale for the unstable regime up to z/L = −0.001, indicating the Bowen ratio with a colour code and separating dry from moist cases with symbols. Dry cases (in circles) are distributed close to the standard similarity function, considered here the Businger formulation re-calculated by Högström (1988). The plot also shows the adjusting functions for dry and humid datasets and Table 4 indicates the values for their respective coefficients. The adjusting function for dry cases (red line) presents minor changes in its α coefficient with respect to Businger's, while β is the same as the value found in Dyer and Bradley (1982) for an additional experiment done in Australia. Nonetheless, the latter falls approximately within the 20% of tolerated discrepancies among the different accepted formulations (Högström 1996).
Therefore we can conclude that the standard similarity equations, both using z/L or Ri B , can be applied to our site when the soil is relatively dry in the daytime. This is also supported by the statistics in Table 3, where it is seen that all methods (including the use of LST instead of T 36 for z/L) display 71-73% of the data below the 20% error and are considered good, while the amount reaches 96-98% if we allow a 50% error. It is clear that the performance of the methods for the sensible heat flux in dry conditions is as good as the one of the momentum flux in all cases.
As an interesting supplementary information, Fig. 9a indicates that there is a relation very close to linearity between the kinematic flux and the temperature gradient between 2 and 0.36 m a.g.l., corresponding to an increase of 0.1 K m s −1 every increase of 0.5 K m −1 in the gradient, regardless of the intensity of stratification. The linearity is less clear when using LST at a height z 0T (Fig. 9b), as there is more dispersion likely linked to the uncertainty of the method for obtaining LST (see Sect. 4 for more details). Still, the fitted line increases approximately 0.1 K m s −1 for every 2.5 K m −1 . The difference between both fitted lines is due to the largest amplitude in the daily cycle of LST with respect to T 36 , generally providing around noon stronger thermal gradients between 2 m and the surface. Nonetheless, this result is very similar qualitatively to the one in Simó et al. (2018), where LST is used as the lowest temperature level.

Night-Time Cases with Either Dry or Moist Soil
The estimation of H at night is done for a smaller dataset since the remaining cases after filtering consist only of 14% of the total, essentially weakly to moderately stable cases. Figure 10a shows that most of the points fall within a range of H between − 10 and − 30 W m −2 , the rest corresponding to strong wind cases in autumn induced by the topography. These selected events have only a 25-26% of cases considered good (error below 20%) and 50-58% below 50% error when using similarity functions depending on z/L or Ri B with T 36 as the lower temperature level (Table 3). The statistics are much worse when LST is used instead, pointing to a failure in either the measurement of this parameter at night probably due to its sensitivity to the particular footprint of the radiometer (Mahrt 2008), or the extension of the temperature profile to a warm surface during most part of the experimental period (Simó et al. 2018). For this regime the wetness of the soil does not appear to be a relevant factor, contrarily to the daytime. The fact that there is condensation or evaporation (whether there is negative or positive L E, respectively) does not seem to imply a distinct behaviour in the results either (Fig. 10a). Figure 10b shows the differences of the estimated H between z/L and Ri-based stability functions, with a dispersion smaller than in Fig. 10a, indicating that the results are independent of the stability parameter. In Fig. 11, the experimental φ h values for stably stratified conditions against a logarithmic z/L > 0.001 is depicted, discriminating the observed L E with colours. The Businger (black) and fitted (green) curves (Table 4) are also included, the latter built upon the experimental points with φ h ≥ 0. The plot shows that for weak stabilities the standard Businger function applies to our data, when counter-gradient cases are excluded (φ h < 0). When stability increases (z/L larger than 0.03), most of the condensation cases appear and the formulation would require a modified equation. Yet the dispersion of φ h for the entire stability range is similar to other analyzed datasets, as in Yagüe et al. (2006).
A plot of the vertical temperature gradient against the EC turbulent sensible heat flux does not display any correlation (not shown). The largest fluxes are due to high winds and are present due to the topographical characteristics of the location. A preliminary conclusion for stably stratified conditions is that the standard similarity formulation is adequate for our site in weakly stable conditions when evaporation prevails, while its applicability in strongly stably stratified situations with dominance of condensation is questionable. However, the few points available in this regime after filtering and excluding the counter-gradient cases suggest to take this conclusion with caution.

Daytime Cases with Moist Soil
The fact that the main discrepancy with the standard theory is linked to the amount of moisture in the soil deserves a more detailed inspection.
As shown in Fig. 3 in spring (between April and June) it is common to have very intense incident solar radiation and a well-watered soil (nearing field capacity as soil V W C is larger than 0.3 m 3 m −3 ), with a heterogeneous distribution of vegetated soil. Many of these cases are characterized by large values of L E and low values of the Bowen ratio. Table 3 indicates that moist cases have a bad behaviour statistically when using the standard stability functions, as the good estimates are just 24-27% compared to the 71-73% of the dry soil cases, as can also be seen in Fig. 7a.
The distinct behaviour of the moist soil cases is evident in Fig. 8, where the corresponding data diverge from the standard curve for values of −z/L smaller than 0.1 or Bowen ratios smaller than 0.5. In such cases most of the energy received from the Sun is used in evaporation and the resulting vertical temperature gradient may not be directly governed by turbulent mixing in the same way as it is when the solar irradiation is smaller or the soil is drier. The fitted curve upon the moist experimental data (blue line) gives a value of φ h at neutrality much higher ((φ h ) z/L=0 = α = 2.1, Table 4) than the Businger expression (α = 0.95). In the same direction, Lee and Buban (2020) provide a value higher than 1.5 for this coefficient with experimental data taken in summer from three towers over soybean fields or grasslands at northern Oklahoma (USA). A rainfall at the beginning of their experiment provided a dataset with high soil water conditions, Bowen ratios close to or lower than one and a high solar irradiation that resemble the environmental conditions of our moist data.
As it is evident in Figs. 7a and 8, for these cases the standard similarity functions estimate fluxes that are about two times larger than the EC observations. Therefore wet soils under intense solar radiation allow larger temperature gradients than would be sustained in dry conditions, as here the turbulence mixing of thermal origin is reduced by the presence of evaporation. The footprint is expected to be short for measurements below 2 m (approximately 200 m around the station), and variations in the soil water content can be assumed to be small for a non-irrigated area. However, one can not rule out that surface heterogeneity could explain partly the obtained results, especially without the presence of additional observations in the very close surroundings. Therefore, this issue is still inconclusive and focus is to be put on the subject in further research efforts by exploring databases for other locations.
The absence of explicit documentation of this effect in the literature might be related to the fact that wet soils with short vegetation and large solar irradiation are not so common [being the aforementioned work of Lee and Buban (2020) an exception]. The pioneering works of Businger or Dyer were made essentially over dry terrain. The reference sites in northern Europe usually lack long episodes of intense solar irradiation, because of their latitude and the more common presence of clouds. In the tropics both conditions can be found, especially in the savannah during and shortly after the rainy season. Instead at the tropical jungles a dense canopy is found above the surface and the physical processes are different. If this preliminary finding is confirmed, adjusted empirical equations for these conditions like the one obtained in Fig. 8 should be applied.

Conclusions
To assess to what degree the standard similarity equations for the turbulent momentum and sensible heat flux apply to complex terrain, they have been computed with screen level values and surface conditions using an 8-month long database obtained at the flat bottom of a mountain valley. Data screening according to MOST conditions and instrumental errors leads to a large reduction of the available cases at night, specially for the sensible heat flux, yielding night-time results that mostly correspond to weakly stratified conditions. A reduced averaging period of few minutes would be more appropriate in a dedicated night-time study due to the usual lack of stationarity for longer time scales, probably increasing also the number of available cases that would pass the quality flag of measured fluxes (see Sect. 2.3).
The estimated momentum fluxes compare well with the EC measurements, especially in the daytime. Night-time comparisons are satisfactory except for weak winds or for strong winds forced by the topography. A key parameter is the aerodynamic roughness length, determined monthly for the site's evolving surface characteristics. Similar results are obtained regardless of the stability parameter used in the stability function (i.e. either the Obukhov length or Richardson number) and of the use of either LST at z 0T or a second level of air temperature close to the surface.
Similarly the estimates of the sensible heat flux are good in relatively dry soil conditions in the daytime, with little sensitivity to the method used in the computation of the stability function. A well defined linear relation between the observed flux and the temperature gradient between 2 and 0.36 m a.g.l. is obtained, and also reproduced for LST although with more scatter due to its larger uncertainty. Therefore, the use of two levels of air temperature is recommended since it also skips the prescription of a thermal roughness length.
At night the estimates for the weakly stable regime are acceptable, independently of the soil moisture content, showing the largest dispersion for weak winds and strong stratification.
The sensible heat flux in the daytime with moist soil and large solar irradiation deserves separated conclusions, as the estimates depart significantly from the measurements. This regime, dominated by large latent heat fluxes, seems to allow vertical temperature gradients larger than those found when the turbulent sensible heat flux exchange is the dominating process.
The results found for this particular site support the use of standard similarity functions over a locally flat surface in moderately inhomogeneous complex terrain for measurements sufficiently close to the ground, sharing similar limitations as in flat and nearly homogeneous surfaces. In contrast, for very complex terrain or over slopes, recent studies have found that MOST fails as fluxes are not constant with height. The large departure of the estimated sensible heat flux for moist soils and large solar irradiation will be further inspected with other databases, together with the assessment of the estimated latent heat flux. Finally, the current study has not considered an important number of nocturnal cases with weak winds where turbulence is not the dominant factor. This subset of data will be the subject of a further work.
projects PCIN-2016-091, CGL2015-65627-C3-1-R (MINECO/FEDER) of the Spanish Government, supplied by the European Regional Development Fund (FEDER). B. Martí is supported by the grant FPI-CAIB (FPI/2165/2018) of the Vicepresidència i Conselleria de Innovació, Recerca i Turisme del Govern de les Illes Balears and the Fons Social Europeu (European Social Fund). The SEB station was rented from the Meteorological Service of Catalonia with funds of the research project PCIN-2014-016-C07-01. We thank Elena Cantero and Fernando Borbón from CENER for their help in the preparation, installation and maintenance of the stations and Javier Sanz for contributing to the direction of the campaign. D. Martínez-Villagrasa acknowledges Ludi Gutiérrez and Ignacio Pérez for their local support during the ALEX17 field campaign. Finally, we wish to thank M. Rotach and another anonymous reviewer for their valuable comments that have contributed to improve the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Data Availability Half-hourly statistics and rawdata from the SEB station are available at Martínez-Villagrasa et al. (2019), while rawdata from the SLS can be found at . Half-hourly processed statistics from SLS are available on reasonable request from the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.