Understanding model biases in the diurnal cycle of evapotranspiration: a case study in Luxembourg

The diurnal forcing of solar radiation is the largest signal within the Earth system and dominates the diurnal cycle of the turbulent heat fluxes and evapotranspiration (λE) over land. Incoming solar radiation (Rsd) also shapes temperature, vapor pressure deficit and wind speed known as important controls on λE. Current process-based λE schemes used in remote sensing and land-surface modeling differ in how these controls on λE are represented and which input variables are required. 15 Here, we analyze how well different surface energy balance schemes are able to reproduce the diurnal cycle and how the diurnal signals of observed input variables actually influence the resulting diurnal pattern of λE. As additional constraint for model evaluation we estimate a linear and a non-linear phase shift component of a surface variable (e.g. λE) to incoming solar radiation. We illustrate our analysis with observations from an eddy covariance station at a temperate grassland site in Luxembourg with a focus on clear sky conditions. During the field campaign in 2015 a summer drought led to a dry-down of 20 soil moisture which allows for studying the effect of wet and dry conditions on the diurnal cycle. We found a remarkable, almost linear relationship of λE with Rsd, which exhibits a significant positive phase lag during wet periods. This phase lag in λE was compensated by a preceding phase lag of the sensible heat flux. Vapor pressure deficit (Da, often used as input for Penman-Monteith based approaches) exhibits a strong phase lag, which is driven by air temperature reflecting large diurnal heat storage changes in the lower atmosphere. This large phase lag in Da, which is not seen in λE, 25 explains why actual and potential evapotranspiration approaches can show systematic deviations from observations at the sub-daily time scale and highlight the need for a time-dependent non-linear compensation through the conductance parameterization. The surface to air temperature gradient used as input in energy balance residual approaches corresponds rather well with its linear response and phase lag to the observed sensible heat flux under both, wet and dry conditions. This simplifies the conductance parameterization and explains the better correlation of these models at the sub-daily time scale. 30 We conclude that the analysis of phase lags at the sub-daily time-scale provides valuable information on the drivers of the surface-atmosphere exchange, which can be used to evaluate and improve the representation of land-atmosphere coupling in land-surface schemes. Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Evapotranspiration and the corresponding latent heat flux (λE) couple the surface water and energy budgets and are of high relevance for water resources assessment. λE is generally limited by four physical factors, (i) the availability of energy mostly supplied by solar radiation, (ii) the availability of and the access to water, (iii) the plant physiology, and (iv) the atmospheric transport of moisture away from the surface (Brutsaert 1982). These different limitations have led to different 5 approaches on how to model λE.
Key approaches either focus on the surface energy balance where the surface-air temperature gradient dominates the flux; or approaches which focus on the moisture transfer limitation where vapor pressure gradients dominate the flux. It is critical to recognize that these two limitations are not independent of each other but rather are shaped by land-atmosphere heat and 10 water exchange and thus covary with each other. The diurnal variation of incoming solar radiation (R sd ) causes a strong diurnal imbalance in surface heating leading to the pronounced diurnal cycles of surface states and fluxes (Oke 1987, Kleidon andRenner, 2017). This heat exchange of the surface with the lower atmosphere thus influences the near-surface air temperature (T a ), skin temperature (T s ), vapor pressure (e a ), soil or canopy saturation water pressure (e s ), vapor pressure deficit (D a ), and wind speed (u), that are being regarded as important controls on λE (e.g. Penman 1948). These interactions 15 are particularly dominant at the diurnal time scale (e.g. De Bruin and Holtslag 1982) and depend on meteorological as well as on surface conditions (Jarvis and McNaughton 1986). Ignoring the interdependence of the surface variables may lead to biases in model parameterizations and compensating errors when evaluating the model performance only with respect to a single variable (Matheny et al., 2014, Best et al., 2015, Santanello et al., 2018. 20 To overcome these problems several authors suggested to better evaluate land-atmosphere (L-A) interactions proposing different multivariate metrics (e.g., Betts 1992, Santanello et al. 2009, Matheny et al., 2014, Zhou et al., 2016b, Santanello et al., 2017. These metrics explore internal relationships between state variables to better characterize key processes and to guide a more systematic exploration and understanding of model deficiencies. There is a strong need to investigate and to derive metrics based on comprehensive observation that characterize the whole land surface-atmosphere system (Wulfmeyer 25 et al. 2018).
Here, we propose an alternative metric, which focuses on the diurnal variation of surface states and fluxes with respect to incoming solar radiation (R sd ). It is well established that R sd , which shapes the diurnal cycle of net radiation (R n ) over land, is strongly correlated with λE. Examples are the successful application of equilibrium evapotranspiration (Schmidt 1915, 30 Priestley-Taylor 1972, Miralles et al., 2011, Renner et al., 2016, the semi-empirical Makkink equation (Makkink 1957, De Bruin and Lablans 1997, De Bruin et al., 2015 and using a fixed fraction of available energy for estimating potential evapotranspiration (Milly andDunne 2016, Maes et al., 2018). Observations of a near linear relationship of λE to R sd Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License. (Jackson et al., 1983) and the need for spatial mapping of λE stimulated remote sensing based applications (Crago, 1996).
Although observations justify this simple behavior, the role of other controls on λE remains unclear (De Bruin and Holtslag 1982, Beljaars and Bosveld 1997, Best et al., 2015Haughton et al., 2016;Zhou and Wang 2016a;Milly and Dunne 2016). This is particularly important for remote sensing based λE models since they use observational states which have an imprint of land-atmosphere interactions. 5 The relevance of other atmospheric controls on λE should result in a deviation from linearity between λE and R sd . There could be a consistent lag between the diurnal patterns of λE versus R sd which would reveal a loop when λE is plotted against solar radiation. There is observational evidence of such consistent loops for heat fluxes plotted against net radiation (Camuffo and Bernadi 1982;Mallick et al., 2015), with many studies showing hysteretic loops of the soil heat flux against 10 net radiation (Fuchs and Hadas, 1972;Santanello and Friedl, 2003;Sun et al., 2013). Camuffo and Bernardi (1982) showed that the magnitude and direction of such hysteretic loops can be estimated by a multi-linear regression of the variable of interest against the forcing variables and its first order time-derivative. This simple model allows to estimate storage effects on diurnal (Sun et al. 2013) to seasonal time scales (Duan and Bastiaansen 2017).

15
Here we use R sd as independent forcing of land-atmosphere exchange to quantify the response of surface heat fluxes and states to R sd as well as the magnitude and direction of phase lags. Thereby we obtain a metric to measure the influence of the various drivers of λE. Process-based models represent our knowledge of the various drivers of λE and therefore should be able to capture the magnitude of hysteretic loops under different conditions (Matheny et al., 2014). Hence, we can use the hysteretic behavior as a pattern-based approach and as an additional criterion to evaluate the performance of different 20 models. We focus on two different approaches to estimate λE. The first approach is based on the energy limitation of λE, using the equilibrium evaporation concept (Schmidt 1915) as formulated by Priestley and Taylor (1972) for potential evaporation. For actual evaporation we focus on one source and two source energy balance schemes (OSEB and TSEB, respectively) which derive λE as residual term of the surface energy balance and parameterize the sensible heat flux by a temperature gradient -resistance description (Kustas and Norman 1996) (referred to as 'temperature gradient scheme'). The 25 second approach is based on the Penman-Monteith (PM hereafter) approach (Monteith 1965), which adds water vapor pressure deficit as a driving gradient (referred to as 'vapor gradient scheme'). We use the widely used FAO Penman-Monteith formulation (Allen et al., 1998) for potential or reference evapotranspiration. For actual evapotranspiration we use a modified PM approach which was formulated by Mallick et al. (2014Mallick et al. ( , 2015Mallick et al. ( , 2016Mallick et al. ( , 2018; (see also Bhattarai et al., 2018) and is termed as a surface temperature initiated closure (STIC). STIC is based on finding the analytical solution of the 30 surface and aerodynamic conductances in the PM equation while simultaneously constraining the surface and aerodynamics conductances through both surface temperature and vapor pressure deficit.
Several inter-comparison studies evaluated the performance of these schemes using observations from different landscapes.
OSEB and TSEB which are often used in remote sensing retrievals of λE have been found to perform comparably well in reproducing tower-based energy flux observations (Timmermans et al. 2007;Choi et al. 2009;French et al., 2015). Yang et al. (2015) compared temperature gradient approaches (including TSEB) with the Penman-Monteith approach (based on vapor pressure gradient only) employed by the MODIS evapotranspiration product (MOD16, Mu et al., 2011) and found 5 strongly reduced capability of MOD16 to estimate spatial variability of evapotranspiration. They concluded that the moisture availability information obtained from relative humidity and vapor pressure deficit is not able to capture the surface water limitations as reflected in surface temperature.
In this study, we focus on the ability of these different evapotranspiration models to reproduce the diurnal cycle of λE under 10 wet and dry conditions. In particular, we assess if significant non-linear relationships in form of hysteretic loops exist, if these change under different wetness conditions and if temperature-gradient and vapor-gradient approaches such as PM are able to reproduce this behavior. Further, we evaluate which input variables of the evapotranspiration schemes show a hysteretic pattern and how these patterns influence the flux estimation. To address these questions, we analyze observations and models with respect to internal functional relationships (pattern-based) and use solar radiation as independent driver of 15 land-atmosphere exchange. We focus on wet vs. dry conditions since this is another critical deficiency identified in previous analyses (e.g. Wilson et al., 2003, Matheny et al., 2014, Zhou et al., 2016b. To ensure similar radiative forcing and avoid variability due to cloud cover we focus the evaluation on clear-sky days. We illustrate our approach on a grassland site in a temperate semi-oceanic climate using surface energy balance observations.

20
The analysis will shed light on the capabilities of process-based evapotranspiration schemes to capture the dynamics of diurnal land-atmosphere exchange. We show that the phase lag of surface states and fluxes reveals important imprints of heat storage processes and how this guides the evaluation of the different approaches for modeling λE. This is important for applications in remote sensing with respect to the choice of observational input variables. In doing so, we provide a further, pattern-based metric to assess land-atmosphere interactions, and, thus guide process-based improvements and calibration of 25 land-surface schemes.

Diurnal patterns and hysteresis loop quantification
We first illustrate the pattern-based evaluation of the diurnal cycle using two hypothetical variables Y 1 and Y 2 , as shown in Figure 1. If a variable (Y 1 ) is in phase with R sd , it shows a linear behavior when plotted against R sd (Fig. 1 b). However, if a 30 variable (Y 2 ) has a time lag with respect to R sd , showing a significant difference between morning and afternoon values, it results in a hysteretic loop. The area inside the loop indicates the magnitude of the phase difference, while the direction of Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License. the loop, marked by an arrow at the morning rising limb in Fig. 1b, indicates if a variable is preceding or lagging R sd in time.
If a variable shows consistently larger values during the afternoon as compared to the morning, this will appear as a counterclockwise (CCW) hysteretic loop indicating a positive phase lag with respect to R sd . A negative phase lag appears as a clockwise (CW) loop.

10
To obtain a quantitative measure of the hysteretic pattern, we use the Camuffo-Bernardi equation (Camuffo and Bernardi 1982), which relates the time series of the response variable Y(t) to the forcing variable R sd (t) and its first order time derivative dR sd (t)/dt: (1) 15 Using multi-linear regression, we obtain the coefficients a, b and c assuming a normal distribution of the residuals ε(t). If Y is linear with R sd , the parameter c should be zero. However, if a consistent pattern such as a hysteretic loop exists, then parameter c should be significantly different from zero. Hence, by using regression analysis we can determine if a significant hysteretic relationship between two variables exists, and if the inclusion of such a non-linear term (with c ≠ 0) would 20 improve the model fit.
Both, the parameters b and c contain information on the phase lag between the two variables. Assuming that R sd is a harmonic function with an angular frequency ω, the phase difference can be estimated by: Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.
To derive the first order time derivative of solar radiation, we use a simple difference between time steps. Since the data we use is available in 30 min time steps (see below), we have 48 time steps per day, thus ω = 48/ (2π). To obtain a phase lag between Y and R sd as a time lag t φ [min] we use: (3) 5 Note that the phase lag estimate t φ is somewhat similar to the relative diurnal centroid metric proposed by Wilson et al., (2003) for the analysis of the timing of heat and mass fluxes. The diurnal centroid identifies the timing of the peak of a variable with respect to local time. Since the peak of R sd is at noon local time both metrics are qualitatively comparable.

Field site and observations 10
The study area is a grassland site in Petit-Nobressart, Luxembourg, situated on a gentle east facing slope. The grassland is used as a hay meadow and had short vegetation of about 10-15 cm as the grass was mowed before the start of the experiment. An Eddy Covariance (EC) station (with the setup described in Wizemann et al. 2015)  distance of 0.2 m to the sonic path was used to measure CO2 and H2O fluctuations. The high-frequency signals were recorded at 10 Hz by a CR3000 data logger and the TK3 software was used to compute turbulent fluxes of sensible heat (H), latent heat (λE) and CO2 (Mauder and Foken, 2015). 20 Downwelling and upwelling shortwave and longwave radiation were obtained by a four-component net-radiation sensor (NR01, Hukseflux, Netherlands). The meteorological variables (air temperature, humidity and precipitation) were monitored with a time resolution of 30 minutes. Sensors measuring soil heat flux (two in 8 cm depth, HFP01, Hukseflux, Netherlands), soil temperature (2, 5, 15, 30 cm, model 107, Campbell Scientific Inc., UK), water content (2. 5,15,30 cm,CS616,Campbell 25 Scientific Inc., UK) and matric potential (5,15,30 cm,model 253,Campbell Scientific Inc.,UK) were installed between the turbulence and radiation measurement devices.
Unfortunately, the two upper temperature probes and matric potential sensors showed data gaps and erroneous values from 30 June until excavation on 23 July, 2015. Thus, the ground heat flux was calculated by the heat flux plate method with correction for heat storage (Massman, 1992) only for the period from 11 June to 30 June, 2015. To still obtain soil heat 30 fluxes for the entire measuring period, additionally harmonic wave analysis (Duchon and Hale, 2012)) of the heat flux plate data was applied. Both methods agreed well for the period before 30 June, so that the latter method should provide reliable ground heat flux values for the entire period until 23 July. Table 1 lists the variables obtained from the EC station and used Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License. in this work. For more details on instrumentation and EC data processing see Ingwersen et al. (2011) and Wizemann et al. (2015).

Derived meteorological variables
We derived the saturated water vapor pressure e s (hPa) by the empirical Magnus equation ( Then, the water vapor pressure of the air e a (hPa) was obtained by using air temperature T a and relative humidity (r H ): e a = e s (T a ) r H /100.
To assess the moisture conditions of each date of the site we used the evaporative fraction f E : 5 Since daily averages can be influenced by single large values of the turbulent fluxes and contain missing values, we estimated a daily f E based on the 30 min values of each day using the following linear regression: 10 where f E is the slope of the linear regression, α its intercept and ε R the residuals. Since we use the fluxes of H and λE without energy balance closure correction we obtain the upper range of f E . 15 Since the sonic anemometer measures friction velocity (u*) and the absolute value of wind speed u = sqrt(U 2 +V 2 ), we estimate the aerodynamic conductance for momentum (u* 2 /u) and the aerodynamic conductance (g ah ) for heat including the excess resistance to heat transfer using an empirical formula by Thom (1972): We chose to use this formula for its simplicity and similar performance than more recent, complex parameterizations (Knauer et al., 2018). Also note that effects of atmospheric stability are accounted for in the first term of Eq. (4).

Energy balance closure gap correction 25
Most EC measurements show that the sum of the observed turbulent heat fluxes are smaller than the available energy and thus do not close the energy balance leaving an energy balance closure gap (Q gap ) (Foken et al., 2008) Q gap = R n -(G + H + λE) 30 For our site we observed on average (H + λE) / (R n -G) = 0.81 with an average gap of 67 W m -2 over the whole duration of the field campaign. These values are in the typical range what is commonly found for grassland sites (Stoy et al., 2013).
To correct the turbulent fluxes for the energy balance closure gap, we use a correction based on the Bowen ratio (B R ) (Twine et al., 2000), which is directly related to f E = 1/(B R +1) to obtain corrected fluxes: We use these corrected fluxes in the further analysis.

Clear-sky day classification
In order to achieve comparable conditions with respect to incoming solar radiation, we identified clear-sky conditions. A clear-sky day was defined by its daily sum of incoming solar radiation being larger than 85% of the potential surface radiation (R sd,pot ), which is a function of latitude and day of year (using R package REddyProc, function fCalcPotRadiation): 15 where t corresponds to each time step of measurement and f diff = 0.78 being a constant factor taking account for atmospheric absorption of solar radiation. 20

One and Two Source Energy Balance models
Thermal remote sensing based models estimate evapotranspiration by solving the surface energy balance and rely on land surface temperature (T s ) information as a key boundary condition (Kustas and Norman 1996). A bulk layer formulation of the soil-plus-canopy sensible heat flux is employed and λE is derived by enforcing the surface energy balance. Hence λE is written as: 25 where ρ is density of air, c p is the specific heat of air at constant pressure, and g ah is the effective aerodynamic conductance of heat that characterizes the transport of sensible heat between the surface and the atmosphere. We obtained T s from the 30 observed longwave emission of the surface T s = (R lu /(σ ε s )) 1/4 with σ = 5.67 x 10 -8 W K -4 the Stefan Boltzmann constant and a surface emissivity ε s = 0.98, which is typical for a grassland and agrees with Brenner et al., (2017).
We use two different approaches which are generally classified as one-and two-source models with regard to the implemented treatment of the energy exchange with the surface. While one-source energy balance models (OSEB) treat the surface as a uniform layer, two-source energy balance models (TSEB) partition temperatures, radiative and energy fluxes into a soil and vegetation component. The one-source approach (OSEB) parameterizes the aerodynamic conductance g ah as follows (e.g. Kalma et al., 2008, Tang et al., 2013: where z u and z t are the measurement heights of wind and air temperature, respectively, z 0m and z 0h are roughness lengths for momentum and heat, respectively, k is the von Kármán constant, d is the displacement height, u is the wind speed and Ψ m 10 and Ψ h are the the integrated Monin-Obukhov (MO) similarity functions which correct for atmospheric stability conditions (Brutsaert 2005, Jiménez et al., 2012. For the investigated grassland site, d and z 0m were calculated as fractions of the vegetation height, h c , with d = 0.65 h c and z 0m = 0.125 h c . The roughness length for heat z 0h was set using the dimensionless parameter kB -1 = ln(z 0m /z 0h ), which was set to 2.3 in accordance with Bastiaanssen et al., 1998. Note that this parameterization of aerodynamic conductance does not explicitly distinguish between bare soil and canopy boundary layer 15 conductance, as it is done in two-source approaches.
In addition to OSEB we applied the Two-Source Energy Balance (TSEB) model developed by Norman et al., (1995), Kustas and Norman (1999). For both, the soil and canopy components a separate energy balance (with different component temperatures) and bulk resistance scheme with different aerodynamic conductance are formulated. Then the energy balance 20 equations are solved iteratively. It starts by assuming that a fraction of the canopy (described by vegetation greenness fraction f g ) transpires at a potential rate as described by the Priestley-Taylor equation (Priestley and Taylor 1972): where α PT is the Priestley-Taylor coefficient (1.26), s is the slope of the saturation water vapor pressure curve and γ is the psychrometric constant. However, the canopy latent heat flux λE c = f g λE PT might be too large and the soil component would become negative (condensation at the soil surface) which is unlikely during daytime conditions. To avoid condensation at the soil surface, the α PT coefficient is reduced incrementally until the soil latent heat flux becomes zero or positive. Once this condition is met, all other energy balance components are updated accordingly to satisfy the energy balance equation. For 30 this study we used a constant vegetation fraction of f c = 0.9 and a greenness fraction f g which was derived from close-up pictures taken at the beginning and the end of the field campaign and linearly interpolated in-between.

Penman-Monteith approach
In the Penman-Monteith approach (Monteith 1965) the inclusion of physiological conductance (g s ) imposes a critical control on λE: In equation (8), the transfer of moisture is linked to a supply-demand reaction where the net available energy (R n -G) is the supplies the energy for evaporation and the vapor pressure deficit of the air, D a [= e s (T a )e a ] is the demand for evaporation from the atmosphere. In the PM approach, the two conductances, the aerodynamic conductance g av and the surface conductance g s to water vapor are unknown. A widely used approach to obtain a reference evapotranspiration estimate from 10 meteorological data is the FAO Penman-Monteith reference evapotranspiration (Allen et al., 1998). It defines the two conductances for a well-watered grass surface with a standard height of h = 0.12 m. The aerodynamic conductance is obtained by a bulk approach (eq. 6) with wind speed u measured at 2m above the surface, d = ⅔ h, z 0m = 0.123h, z 0h = 0.1z 0m yielding g av = u/208 (Box 4 in Allen et al., 1998). Surface conductance is fixed at a constant g s = 1/70 m/s. While the FAO estimate is typically intended for estimates of the reference evaporation for well-watered grass on a daily basis, we use it 15 here as a reference for comparison on a subdaily scale. In order to understand the effect of the aerodynamic conductance parameterizations we add another reference evapotranspiration estimate in which the aerodynamic conductance is given by eq. (4) using observations of friction velocity and wind speed, but keeping g s fixed.

Penman-Monteith based Surface Temperature Initiated Closure (STIC) (version STIC1.2)
In order to estimate an actual evapotranspiration rate from meteorological data we employ a method (STIC1.2 hereafter 20 referred to as STIC), which is based on the PM equation, but which in addition integrates surface temperature information.
The STIC methodology is based on finding analytical solutions for the two unknown conductances to directly estimate λE (Mallick et al., 2016. STIC is a one-dimensional physically-based SEB model that treats the vegetation-substrate complex as a single unit (Mallick et al., 2016;Bhattarai et al., 2018). The fundamental assumption in STIC is the first order dependency of g a and g s on soil moisture through T s and on environmental variables through T a , D a , and net radiation. 25 Thereby, surface temperature is assumed to provide information on water-limitation which is linked to the advection-aridity hypothesis (Brutsaert and Stricker 1979). In STIC, no wind speed is required as input data, as opposed to the temperature gradient approaches, but vapor pressure of the air and its saturation value become critical input variables, see Table 2 for an overview. A detailed description of STIC version 1.2 is available in Mallick et al. (2016Mallick et al. ( , 2018 and Bhattarai et al. (2018).

Daily clear-sky and moisture classification
The field campaign was conducted during an exceptionally warm and dry period characterized by clear sky conditions with remarkably high air temperatures with daily maxima above 30°C and little precipitation. Compared to the climatic normal  the precipitation deficit in this region was -44% in June and -41% in July, respectively (source: meteorological 5 station Arsorf, Administration des services techniques de l'agriculture (ASTA)). The air temperature anomaly was higher in July (1.9°C) than in June (0.7°C) (source: meteorological station Clemency, ASTA). The soil water content decreased and parts of the site, especially in the upper part, showed clear signs of vegetation water stress (see Brenner et al., (2017) for an analysis of the spatial heterogeneity of water limitation). However, the dry period was interrupted by a few but strong rainfall events, which significantly changed soil moisture and thus f E with time (Fig. 3a). Based on the observed f E we 10 classified dry days with f E < 0.5 and wet days with f E > 0.6. This separation of dry and wet days is also reflected in the top soil moisture conditions (measured at 5 cm depth) as shown in Fig 3b. precipitation. Panel a) shows the daily time series and panel b) the relationship of f E to soil moisture used to classify "wet" and "dry" days depending on f E > 0.6 or f E < 0.5, respectively. Sunny days are defined using a threshold of 85% of R sd to potential radiation and are marked with solid symbols, with blue circles referring to wet and red squares to dry days. Top soil moisture measured at 5 cm below surface is shown.

20
Based on the classification into wet and dry days under clear-sky conditions we computed composites of the diurnal cycle for each hour. By using only sunny days we aim to achieve similar conditions with respect to downwelling shortwave radiation (R sd ). Figure 4a confirms that R sd and net radiation (R n ) had very similar diurnal cycles and magnitudes for the wet and dry days. However, the downwelling longwave radiation R ld and the soil heat flux were somewhat higher under wet conditions (Fig. 4a). The higher R ld is related to a higher air temperatures and air vapor pressures observed under wet conditions (Fig. 4b), which may explain the greater value of R ld by affecting the atmospheric emissivity for longwave radiative exchange. This has also an impact on the minimum temperatures both for air and skin temperature, which are higher under wet conditions and lower under dry conditions (Fig. 4b). Hence, although we achieve fairly similar conditions 5 for shortwave radiation under wet and dry conditions, we observed a small but significant difference in the longwave radiative exchange.

Diurnal cycle of evapotranspiration under wet and dry conditions
Next, we evaluate how the different evapotranspiration schemes are able to reproduce the fluxes during wet and dry conditions under similar R sd forcing. Figure 5 shows the average diurnal cycle of observations and models for λE. The 15 observations showed a significant difference in λE between dry and wet conditions, with the maximum value of λE of about 200 W m -2 for dry and 350 W m -2 under wet conditions, which amounts to a mean difference of 100 W m -2 for daylight conditions ( Table 3). As reference, we also included two common formulations of potential evapotranspiration, the Priestley-Taylor potential evapotranspiration (PT) and the FAO Penman-Monteith reference evapotranspiration (FAO-PM). Both do not account for water limitation and show a marginal difference of 10 W m -2 between wet and dry conditions. While FAO-20 PM yielded lower mean conditions than PT, it showed lower correlation and RMSE as compared to PT (Table 3). We find that all models for actual λE (rather than PT or FAO-PM) showed differences in λE between wet and dry conditions. Both conditions. In contrast, STIC captured the low λE magnitude under dry conditions (f E <0.5) but underestimated λE under wet conditions (for f E >0.6). Table 3 shows the statistical metrics of the model performances with respect to the Bowen ratio corrected λE. In general, both OSEB and TSEB produced mean λE values within the range of 96 -98% (255 W m -2 and 259 W m -2 ) of the observed 5 λE (264 W m -2 ) in wet conditions, while mean λE from STIC was within 83% (218 W m -2 ) of observed λE for the same conditions. However, for the dry conditions, simulated λE from STIC (180 W m -2 ) was 91% of the observed mean λE (164 W m -2 ), while the simulated mean λE from OSEB and TSEB was 77 -78% of the observed mean λE. Overall, the three models captured 86% (OSEB), 88% (TSEB), and 95% (STIC) of the observed mean λE. Results show that under wet conditions, RMSE of the OSEB /TSEB models is well within the range of the errors when compared with the uncorrected 10 λE, whereas STIC showed relatively higher RMSE. However, under dry conditions the RMSE of OSEB /TSEB models was found to be larger than for STIC. For the entire observation period the three models produced comparable RMSE (41 -46 W m -2 ) but with different correlation. STIC produced relatively low correlation (r 2 = 0.72) as compared to the other two models

Diurnal patterns of evapotranspiration
The evaluation of the diurnal cycle shows that λE was strongly related to the incoming solar radiation, emphasizing that R sd is the dominant driver of λE (Fig. 6). However, under wet conditions we found a marked and consistent difference between morning and afternoon in λE forming a CCW hysteresis loop (Fig. 6b). Using the Camuffo-Bernardi regression we found a significant phase lag for the BR corrected flux (λE BRC ) with a mean t φ = 15 min under wet conditions and no significant lag 10 under dry conditions ( Fig. 7 and Table 4). The uncorrected observations showed a slightly lower wet-dry difference.
The two potential evapotranspiration estimates showed large differences in their phase lag. While the PT estimate showed a small hysteretic loop with a phase lag between t φ = 6-9 min, the FAO-PM estimate showed a substantial loop with a phase lag of t φ = 31 min. This large phase lag of the FAO-PM estimate is very similar to the phase lag when we use a constant g s in 15 the PM equation but with g av obtained from Eq. (4) using friction velocity observations (Table 4). The temperature gradient schemes (OSEB and TSEB) reproduced the observed phase lag relatively well (mean t φ = 9 min for wet and around 0 for dry conditions). However the temperature-vapor gradient scheme (STIC) showed relatively larger phase lags under both dry and wet conditions (t φ = 14-20 min) (Fig. 7, Table 4).

Diurnal patterns of observed fluxes and states
In order to understand the diurnal patterns of λE we also analyzed the hysteresis loops of the observed surface energy balance components [λE = (R n -G) -H] with respect to R sd (Figure 8). Generally, there was only a small hysteresis in the available energy (R n -G) ( Table 4). The turbulent heat fluxes showed significant hysteresis under wet conditions but not 5 under dry conditions. Interestingly, under wet conditions the CCW hysteresis of λE with a phase lag (mean t φ = 15 min) was mostly compensated by a CW hysteresis of H (mean t φ = -22 min) (Figure 8 and Table 4). This compensation is an outcome of net available energy (R n -G) showing little hysteresis for both conditions.

15
We next analyzed the bulk sensible heat flux formulation used in the OSEB and TSEB models to understand how the observations of temperature and the inferred aerodynamic conductances are related to each other. The diurnal patterns of both air and surface temperature revealed a strong CCW hysteresis with R sd (Fig. 9). Air temperature showed a more pronounced hysteretic loop than surface temperature, and with a triangular shape with higher values during the afternoon 20 when solar radiation reduces. Interestingly, the surface-air temperature gradient, being the driving gradient for the sensible Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.
heat flux, showed much less hysteretic behavior. The hysteresis is in a clockwise direction, with a higher gradient in the morning hours compared to the afternoon. It had a similar phase lag as H (see Table 4).

Figure 9: Diurnal patterns of observed anomalies in surface temperature (T s ), air temperature at 2m (T a ), and their gradient (T s -T a ) for (a) dry and (b) wet days. Both T a and T s show a pronounced CCW hysteresis, but the form of the hysteretic loop is
5 significantly different, with air temperature featuring a more pronounced, triangular shape with afternoon values almost independent of incoming solar radiation. The temperature gradient, however, shows a much smaller CW hysteretic loop. Note that the temperature gradient is comparatively higher in the morning than in the afternoon, corresponding to the diurnal course of the sensible heat flux (cf. Fig. 8). 10 We further analyzed different formulations of the aerodynamic conductance (g a ) directly inferred from measurements and from how these are represented in the models evaluated here (FAO-PM, OSEB, TSEB, STIC). We inferred the aerodynamic conductance from observations in three different ways: Firstly, we used the Eddy covariance measurements of friction velocity (u*) and wind speed (u) to estimate the aerodynamic conductance for momentum (g am = u* 2 /u). We then used the 15 empirical formula by Thom (1972) to calculate the aerodynamic conductance for heat including the excess resistance to heat transfer (Eq. 4). Thirdly, we inferred the aerodynamic conductance from the observed sensible heat flux (H BRC ) and temperature gradient (T s -T a ) by inverting H BRC using Eq. (5). The FAO-PM describes the aerodynamic conductance with a simple linear relationship to wind speed. OSEB and TSEB estimates the aerodynamic conductance to heat (g ah ), while STIC estimates the conductance to water vapor (g av ). Thus by comparing these different conducatance eatimates we assume 20 similarity between the fluxes.
The different estimates for the aerodynamic conductance are compared to each other in Fig. 10 for midday conditions.
Although the three observation-based estimates show some variations in the absolute value of the aerodynamic conductance, Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License. they consistently showed a significantly greater conductance for dry days compared to wet days, suggesting a stronger aerodynamic exchange between the surface and the atmosphere under dry conditions. This difference in aerodynamic conductance is partly reproduced by the simple FAO-PM scheme which means that the median wind speed was higher under the drier conditions. The temperature-gradient schemes (OSEB and TSEB) reproduce the wet-dry difference rather well, which also use wind speed but rely on MOST similarity and stability correction. STIC which does not use wind speed did 5 not show any significant differences in g av between wet and dry conditions. Finally we analyze the diurnal patterns of the vapor pressure deficit D a = e s (T a )e a which is a critical driver of the latent heat flux in the PM equation. Since D a is derived from the observations, we analyzed its diurnal patterns in Fig. 11. We found that 15 the vapor pressure in the air remained fairly constant during the day, hence it did not co-vary with R sd , and only showed a small CW hysteresis with higher vapor pressure during the morning than during the afternoon. The saturation vapor pressure, which is a function of air temperature, however, showed a distinct and large CCW hysteresis loop with respect to R sd , which Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.
is consistent with the large hysteresis in air temperature (Figs. 9 and 12). As a consequence, D a also showed a distinct and large CCW hysteresis with a large phase lag of t φ =~ 150 min (see Table 4). This large hysteresis and phase lag is consistent with the respective characteristics of air temperature, but not with those of the temperature gradient. Furthermore, we note that the phase lag in D a did not show any significant influence of wetness, while the phase lag of the temperature gradient became more negative under wet conditions (Fig. 12, Table 4). It would thus seem that the bias in PM-based estimates 5 identified here may relate to a too pronounced role of D a in the evapotranspiration estimate.

Dominant controls of λE at the diurnal cycle
Our analysis of the diurnal cycle showed that λE follows the diurnal course of incoming solar radiation, explaining most of the variance in λE. However, a significant non-linearity in form of a phase lag between λE and R sd was detected which showed larger λE for the same R sd in the afternoon as in the morning. We found that the lag in λE is accompanied by a preceding phase lag of the sensible heat flux, while the other surface energy balance components (e.g., net radiation and soil 10 heat flux) revealed very small phase lags with R sd . Hence, there is compensation between the phase shifts of sensible and latent heat fluxes, which becomes more apparent under the wet conditions. Our results are consistent with comprehensive Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.
FLUXNET studies of Wilson et al. (2003) and Nelson et al. (2018) which used a different metric (median centroid) for assessing diurnal phase shifts. Wilson et al. (2003) found that H precedes λE at most sites, with the exception of sites in a Mediterranean climate. Using the FLUXNET2015 dataset, Nelson et al. (2018) found that the median centroid of λE occurs predominantly in the afternoon across all plant functional types when f E > 0.35, while for very dry conditions (f E < 0.2) a shift of the λE centroid towards the morning was found. This indicates that our results are not just applicable to Luxembourg, 5 but are a general phenomenon which justifies a wider interpretation within temperate climates.
The obtained phase lags of the surface fluxes and variables allow for a process-based insight into the diurnal heat exchange of the surface with the atmosphere. Since there is only limited heat storage in the surface layer itself, which explains the small phase lags of the heat fluxes, the heating imbalance caused by solar radiation must be effectively redistributed. While the soil heat flux is limited by the relatively slow heat conduction, the lower atmosphere acts as efficient heat storage. Thus, 10 the lower atmosphere is effectively heated by surface longwave emission and the sensible heat flux, which in combination with the diurnal cycle of vertically transported turbulent kinetic energy (TKE) leads to the development of the convective planetary boundary layer (CBL) (e.g., Oke 1987). The changes in heat storage in the CBL are reflected by the very large phase lags for air temperature and longwave downwelling radiation, which both have a phase lag of about 2.5h. This large phase lag of air temperature then shapes (i) the vertical surface-air temperature gradient, which drives the sensible heat flux, 15 and (ii) the vapor pressure deficit of the air. Despite the complexity of processes within the convective boundary layer, including the morning transition and entrainment at its top, we find that all surface energy components correlate strongly with solar radiation (Table 4). What this suggests is that the state of the surface-atmosphere system is predominantly shaped by fluxes, particularly by solar radiation as its primary driver, with the state in terms of temperatures and humidity gradients adjusting to these fluxes, rather than the reverse, where the state (in terms of temperature and humidity) drives the fluxes. 20 These interactions are also affected by soil water availability, as reflected in the phase lags. This is most clearly seen for the surface-air temperature gradient and the sensible heat flux, whose phase lag is two times larger for wet than for dry days.
This means that for the same solar radiation forcing we find higher values of the sensible heat flux in the morning than in the afternoon. The effect of water availability is also seen for the phase lag of the latent heat flux and to a lesser extent for the 25 soil heat flux. These findings agree well with studies which use the diurnal centroid method, showing that moisture limitation decreases the lag in timing of λE (Wilson et al. 2003, Xiang et al., 2017, Nelson et al., 2018. The phase shift of D a might enhance evaporation at the cost of the sensible heat flux during the afternoon under sufficient moisture availability. However, under drier conditions, our findings suggest that the surface heats more strongly and generates more buoyancy, which is reflected by higher aerodynamic conductances as compared to the wet conditions (Fig. 10). The larger aerodynamic 30 conductance would then enable a more effective sensible heat exchange, and would thus lower the phase difference between the sensible and latent heat fluxes.
Note, that our interpretation disregards the effects of horizontal advection of moisture and temperature. Events of strong advection, e.g. of temperature can add heat to the surface energy balance and thus alter the diurnal cycle. Similarly, events of dry air advection may enhance local λE at the cost of the sensible heat flux. Since we used composite averages and statistics over a set of days we aim to reduce the impact of such advective events. We expect that it is unlikely that such events occurred throughout all wet / dry days in a consistent manner. 5

Using phase lags to identify model biases
Our comparison of different modeling approaches shows that by using phase lags, one can identify biases in evapotranspiration parameterizations and relate these towards processes for a better understanding of surface-atmosphere interactions under different conditions of water availability. One of our main findings is that the surface energy balance 10 fluxes and the temperature gradient have a comparatively small phase lag to the incoming solar radiation, while air temperature and vapor pressure deficit have substantial phase lags. This difference in phase lags can then be used to infer biases in estimates of evapotranspiration. In our application of this approach to observations of one site in a temperate climate we found that evapotranspiration exhibits a comparatively small phase lag, indicating that it was dominantly driven by solar radiation and temperature gradients, and not by air temperature and the water vapor pressure deficit. This 15 interpretation is consistent with studies of non-water-stressed evapotranspiration that is best represented by potential evapotranspiration schemes which are primarily driven by net radiation, as demonstrated for FLUXNET observations by Maes et al., (2018) and for climate model simulations by Milly and Dunne, (2016). Milly and Dunne (2016) interpreted these findings in terms of strong feedbacks between the surface and the atmosphere, which couple the surface variables and which result in a top-down energy constraint that is well captured by energy-only formulations. 20 Our analysis allows to better understand the relevance of the feedbacks which occur at a sub-daily time scale. These feedbacks are driven by the redistribution of heat gained by absorption of solar radiation at the surface, which causes a significant co-variation of the input variables to incoming solar radiation (Table 4). This is especially important for the vapor pressure deficit of the air which acts as a driver of λE and is also known to affect the stomatal conductance (Jarvis 1976, 25 Jarvis andMcNaughton 1986). De Bruin and Holtslag (1982) showed that a positive correlation between D a and R n allows simplifying the complex PM equation to a form similar to equilibrium evaporation (Eq. 7) with net radiation as the dominant driver. Therefore, simpler, energy based formulations for λE show similar performance as PM based approaches, but with less input parameters Holtslag 1982, Beljaars andBosveld 1997). The challenge of the PM equation is then a parameterization of the conductances, which must capture the feedbacks included in the input data. Since the co-variation 30 originates from the diurnal redistribution of heat, a mismatch would then clearly be seen at the diurnal timescale. Hence by focusing on the internal relationship of the modeled λE to R sd at the diurnal time scale we found that (i) the Penman-Monteith based approaches showed a consistently larger phase lag than what was actually observed and (ii) these approaches did not Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.
show a reduction of the phase lag under dry conditions. The PM approaches use the vapor pressure deficit as input variable which showed a substantial hysteresis loop in the order of 2.5 h lagging R sd . This is due to the temperature dependency of the saturation vapor pressure, while the actual vapor pressure shows no relationship with R sd . Besides D a , all other input variables to the Penman-Monteith approaches used here (both FAO and STIC) showed minor phase lags with respect to R sd .
Since the surface conductance in FAO Penman-Monteith is fixed with time, the resulting prediction of potential λE showed a 5 significant and large phase lag in the order of 0.5 h. Even when we use the observed aerodynamic conductance as input, the effect remains the same, which emphasizes that a constant surface conductance results in biases in the diurnal cycle of λE. In contrast to assuming a constant g s , STIC computes λE through analytical estimation of g s and g av from the information of both, the surface-air temperature gradient and the vapor pressure deficit. Given that the lag of D a to R sd was similar for both dry and wet conditions (Fig. 12), λE from STIC also revealed a similar pattern and there was no substantial difference in 10 phase lag of λE from STIC between the wet and dry conditions. Hence, our analysis indicates that the PM-based approaches used here overestimated the effect of water vapor deficit on actual evapotranspiration, which, in the end, reflects in the estimation of the surface and the aerodynamic conductance to water vapor.
The temperature-gradient approaches used here (OSEB and TSEB) are structurally different from the PM approaches, since 15 they infer λE from the residual of the surface energy balance and thus do not explicitly deal with the aerodynamic and surface conductance of water vapor. The phase lag analysis of the environmental variables used to drive the predictive models of λE helped to identify an important benefit of the temperature-gradient approaches over the Penman-Monteith based approach.
The temperature-gradient approaches employ the vertical temperature gradient (T s -T a ) which showed a significant counter-20 clockwise, i.e. a leading hysteretic loop, which is in the order of the phase shift detected for the sensible heat flux (Fig. 12).
In addition, there is a distinct and significant increase of the phase shift in both the temperature gradient and the sensible heat flux under the wet conditions. Hence, the temperature gradient as input contains valuable information on water limitation in terms of the magnitude (i.e. the slope of (T s -T a ) to R sd ) and the diurnal phase lag (see Table 4).

25
While the PM approaches must identify two conductances simultaneously, the temperature-gradient approaches only need to parameterize the aerodynamic conductance to heat (g ah ) using wind-speed as input. Thereby we found that these approaches agreed well with the approximated g ah from the EC tower, which shows an enhanced conductance under dry conditions. Contrarily, the diagnosed g av from STIC did not show substantial differences between dry and wet conditions, pointing to the difficulty of the analytical approach and its associated assumptions to identify two bulk conductances parameters from the 30 available radiometric and meteorological data  for the climatic conditions in which these were evaluated here.
Note, that we evaluated a temperate grassland site which experienced an exceptional summer drought. Thereby, the evaporative fraction did not decline below 0.3. In semi-arid ecosystems the evaporative fraction may decrease substantially below 0.3 and Nelson et al., 2018 showed that there is a morning shift of λE (analogous to a negative phase lag) under very dry conditions (f E < 0.2). This points towards a different stomatal regulation changing the diurnal course of surface conductance. While it was shown by Bhattarai et al., (2018) that STIC performs well also under semi-arid conditions, 5 temperature-gradient approaches can show larger biases under semi-arid conditions (Morillas et al., 2013). The difficulty of temperature-gradient approaches are predominantly in the parameterization of aerodynamic conductance of heat which becomes more challenging under these very dry conditions (Kustas et al., 2016).
The relevance of the diurnal time scale for the problem of surface conductance parameterizations was already highlighted by 10 Matheny et al., 2014. However, they and others evaluated the diurnal patterns of the hysteretic loops between λE and D a (see also Zhang et al., 2014, Zheng et al., 2014. Given that solar radiation is the cause of the strong L-A feedbacks at the diurnal time scale we believe that solar radiation is better suited as a reference variable than D a . Our results show that the new metric of the phase lag of heat fluxes and surface states to incoming solar radiation reveals important biases of evapotranspiration schemes often used in remote sensing. These biases may well be compensated for at a longer time scales (Matheny et al.,15 2014) but would lead to biased sensitivities with respect to climate change (Milly and Dunne 2016). Here, we applied the phase lag metric to observationally driven evapotranspiration schemes. In the future, we plan to apply these new metrics based on hysteretic loops to model outputs of land-surface models (such as NOAH-MP, Niu et al., 2011) as well as of fully coupled surface-atmosphere simulations in order to detect and to identify errors in the parameterization of state-of-the-art LSMs. 20

Conclusions
We quantified the non-linear relationship of evapotranspiration to incoming solar radiation in terms of a phase lag at the diurnal time scale. Our findings have practical implications for remote sensing based λE retrievals and for land-surface model evaluation and calibration. We evaluated three structurally different schemes which are used in remote sensing based applications for a temperate grass site which experienced a summer drought. Our results showed that temperature gradient 25 schemes show a higher correlation and thus better represent the diurnal cycle than a combined temperature-vapor gradient scheme. What this means is that for our site, the temperature gradient between the surface and the near-surface atmosphere contained more important information of water limitation than the water vapor pressure deficit (which is also difficult to retrieve from remote sensing (Kalma et al., 2008, Yang et al., 2015). Furthermore, schemes which use vapor pressure deficit as additional input (such as the Penman-Monteith formulation), require a dynamic, i.e. time dependent characterization of 30 surface conductance to account for the strong phase lag in vapor pressure deficit. Hence, our results suggest that simpler λE approaches based on the surface energy balance and surface temperature may be more suitable to estimate evapotranspiration Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License. from observational data (e.g. remote sensing data) in climates without substantial water stress. It would seem that these surface observations already contain the imprint of land-atmosphere interactions, whereas in the case of coupled land surface-atmosphere models these interactions are explicitly resolved. Hence, detailed models of aerodynamic conductance and its interaction with the environment are of crucial importance for skillful climate predictions including the carbon cycle (Prentice et al., 2014, Wolf et al., 2016Konings et al., 2017). 5 We suggest that an evaluation of these schemes should be based on the diurnal time scale, because a land-atmosphere exchange scheme must accomplish a balance between the surface energy balance with small imprints of heat storage changes and the lower atmosphere with strong imprints of heat storage changes (Kleidon and Renner 2017). Although a mismatch of the diurnal patterns may not be seen at the aggregated time scales of days and months, it may lead to biased model 10 sensitivities (Matheny et al., 2014) which would be crucial for assessing climate change impacts (Milly and Dunne 2016).
Here, we analyzed observationally driven evapotranspiration schemes and their inputs, which revealed an apparent energy constraint. This constraint, which appears as strong correlation of surface fluxes and gradients to incoming solar radiation should be correctly represented by any land-surface model which resolves the land-atmosphere interaction. While this may sound trivial, recent benchmarking studies showed that current state-of-the-art land-surface models have difficulties to 15 represent the strong link of turbulent heat fluxes to solar radiation (Best et al., 2015;Haughton et al., 2016). Our findings provide an explanation of why this is so difficult to achieve and we suggest that further information is gained by evaluating land surface schemes in terms of phase lags in surface fluxes and states such as the sensible and soil heat flux including the diurnal dynamics of surface and air temperature. Correctly representing these metrics will lead towards a more accurate representation of the diurnal heat and mass exchange of the land with the atmosphere. 20 MR merged the data and performed data analysis. MR and LC developed the phase lag computation. JW provided ancillary simulation data. IT provided climate information. MR prepared the manuscript with contributions from all co-authors. 5

Competing interests
The authors declare that they have no conflict of interest.

Special issue statement
This draft shall contribute to the upcoming Joint special issue in HESS (lead journal) and ESSD: Linking landscape organisation and hydrological functioning: from hypotheses and observations to concepts, models and 10 understanding

Acknowledgements
This study was supported by the German Research Foundation (DFG) through funding of the research unit "From Catchments as Organised Systems to Models based on Dynamic Functional Units -CAOS" (FOR 1598) within the subproject "Understanding and characterizing land surface-atmosphere exchange and feedbacks" (Project number 182331427  Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2018-310 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 9 July 2018 c Author(s) 2018. CC BY 4.0 License.