Articles | Volume 23, issue 1
Research article
28 Jan 2019
Research article |  | 28 Jan 2019

Using phase lags to evaluate model biases in simulating the diurnal cycle of evapotranspiration: a case study in Luxembourg

Maik Renner, Claire Brenner, Kaniska Mallick, Hans-Dieter Wizemann, Luigi Conte, Ivonne Trebs, Jianhui Wei, Volker Wulfmeyer, Karsten Schulz, and Axel Kleidon

While modeling approaches of evapotranspiration (λE) perform reasonably well when evaluated at daily or monthly timescales, they can show systematic deviations at the sub-daily timescale, which results in potential biases in modeled λE to global climate change. Here we decompose the diurnal variation of heat fluxes and meteorological variables into their direct response to incoming solar radiation (Rsd) and a phase shift to Rsd. We analyze data from an eddy-covariance (EC) station at a temperate grassland site, which experienced a pronounced summer drought. We employ three structurally different modeling approaches of λE, which are used in remote sensing retrievals, and quantify how well these models represent the observed diurnal cycle under clear-sky conditions. We find that energy balance residual approaches, which use the surface-to-air temperature gradient as input, are able to reproduce the reduction of the phase lag from wet to dry conditions. However, approaches which use the vapor pressure deficit (Da) as the driving gradient (Penman–Monteith) show significant deviations from the observed phase lags, which is found to depend on the parameterization of surface conductance to water vapor. This is due to the typically strong phase lag of 2–3 h of Da, while the observed phase lag of λE is only on the order of 15 min. In contrast, the temperature gradient shows phase differences in agreement with the sensible heat flux and represents the wet–dry difference rather well. We conclude that phase lags contain important information on the different mechanisms of diurnal heat storage and exchange and, thus, allow a process-based insight to improve the representation of land–atmosphere (L–A) interactions in models.

1 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 approaches on how to model λE.

Key approaches either focus on the surface energy balance where the surface-to-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 water exchange and thus covary with each other. The diurnal variation of incoming solar radiation (Rsd) causes a strong diurnal imbalance in surface heating leading to the pronounced diurnal cycles of surface states and fluxes (Oke, 1987; Kleidon and Renner, 2017). This heat exchange of the surface with the lower atmosphere thus influences the near-surface air temperature (Ta), skin temperature (Ts), vapor pressure (ea), soil or canopy saturation water pressure (es), vapor pressure deficit (Da), and wind speed (u), which are being regarded as important controls on λE (e.g., Penman, 1948). These interactions are particularly dominant at the diurnal timescale (e.g., De Bruin and Holtslag, 1982) and depend on meteorological as well as on surface conditions (Jarvis and McNaughton, 1986; van Heerwaarden et al., 2010). 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).

There is a strong need to investigate and to derive metrics based on comprehensive observations that characterize the whole land-surface–atmosphere system (Wulfmeyer et al., 2018). Several authors proposed different multivariate metrics to better evaluate land–atmosphere (L–A) interactions in observations and models. Generally, 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. A number of metrics focus on the diurnal evolution of the heat and moisture budgets in the planetary boundary layer (e.g., Betts, 1992; Santanello et al., 2009, 2018). Also statistical metrics exploring the strength of linear relationships between surface heat fluxes and states to surface radiation components have been employed to evaluate the performance of reanalysis with observations (Zhou and Wang, 2016; Zhou et al., 2017, 2018).

Furthermore, there are pattern-based metrics which focus on nonlinear interactions at the diurnal timescale. Wilson et al. (2003) proposed the method of a diurnal centroid to measure the timing of the surface heat fluxes and their timing difference, which was more recently used by Nelson et al. (2018) to quantify the timing of evapotranspiration under different dryness condition for the FLUXNET dataset. In contrast, Matheny et al. (2014) and Zhang et al. (2014) explored the diurnal relationship of the latent heat flux to vapor pressure deficit showing a pronounced hysteresis loop. Zheng et al. (2014) also included air temperature and net radiation as reference variables and showed that the hysteresis loops of λE to Da or Ta are large, while there are only small hysteresis effects when Rn was used. Hysteresis loops have also been found when heat fluxes are plotted against net radiation (Camuffo and Bernardi, 1982; Mallick et al., 2015), with many studies showing hysteretic loops of the soil heat flux against net radiation (Fuchs and Hadas, 1972; Santanello and Friedl, 2003; Sun et al., 2013). The presence of a hysteresis loop indicates that there is a time-dependent nonlinear control on the variable of interest, typically induced by heat storage processes. Camuffo and Bernardi (1982) showed that the magnitude and direction of such hysteretic loops can be estimated by a multilinear regression of the variable of interest against the forcing variables and its first-order time derivative. This simple model allows estimating storage effects on diurnal (Sun et al., 2013) to seasonal timescales (Duan and Bastiaansen, 2017).

Here, we choose the Camuffo and Bernardi (1982) model because it provides an objective measure of the magnitude of hysteresis loops and it allows for an assessment of statistical significance. We extend the Camuffo and Bernardi (1982) model in two ways.

First, we use incoming solar radiation (Rsd) as a reference variable instead of net radiation to estimate the phase lag of surface heat flux observations and models. And secondly, we use a harmonic transformation of the Camuffo and Bernardi (1982) regression model to estimate the phase lag in time units. This extension allows us to compare the diurnal phase lag signatures of the different model inputs and how these influence the resulting diurnal course of the latent heat flux estimate.

We specifically choose incoming solar radiation Rsd as the reference for the phase-shift analysis, since Rsd can be regarded as an independent forcing of the surface energy balance (e.g., Ohmura, 2014):

(1) R sd ( 1 - α ) + R ld - H - λ E - G = σ T 4 + m ,

with surface albedo α, incoming longwave radiation Rld, sensible heat flux H, latent heat flux λE, the conductive soil heat flux G, the outgoing longwave radiation σT4, and storage terms of the surface layer summarized in m. This formulation of the surface energy balance provides the direction of the energy exchange processes at the surface, illustrating that the terms on the right-hand side depend on heat fluxes on the left-hand side of Eq. (1) (Ohmura, 2014). As a consequence, the term net radiation Rn, which resembles the radiation budget of the shortwave and longwave components, Rn=Rsd(1-α)+Rld-σT4, cannot be regarded as an independent surface forcing. Consequently, we choose Rsd instead of Rn or RnG as the reference variable for the phase-shift analysis of the latent heat flux and the main input variables of evapotranspiration model approaches.

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 the residual term of the surface energy balance and parameterize the sensible heat flux by a resistance description of the surface-to-air temperature gradient (Kustas and Norman, 1999). The 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 the “vapor-gradient scheme”). We use the widely used Food and Agriculture Organization of the United Nations (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. (2014, 2015, 2016, 2018) (see also Bhattarai et al., 2018) and is termed as a the Surface Temperature Initiated Closure (STIC). STIC is based on finding the analytical solution of the surface and aerodynamic conductances in the PM equation while simultaneously constraining the surface and aerodynamic 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 strongly reduced capability of MOD16 to estimate spatial variability of evapotranspiration. They concluded that the moisture availability information obtained from the relative humidity and vapor pressure deficit of the air 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 wet and dry conditions. In particular, we assess if significant nonlinear relationships in the 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 an independent driver of 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 and Wang, 2016). 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.

Figure 1Illustration of a pattern-based evaluation of the diurnal cycle. (a) shows the diurnal cycle of Rsd under clear-sky conditions and the diurnal cycle of two variables Y1 and Y2, one in phase with Rsd and another lagging Rsd. (b) illustrates the relationship of these variables when plotted against Rsd. The bold arrow indicates the direction of the loop and the area inside the hysteresis describes the magnitude of the phase shift.


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 land-surface schemes.

2 Methods and data

2.1 Diurnal patterns and hysteresis loop quantification

We first illustrate the pattern-based evaluation of the diurnal cycle using two hypothetical variables Y1 and Y2, as shown in Fig. 1. If a variable (Y1) is in phase with Rsd, it shows a linear behavior when plotted against Rsd (Fig. 1b). However, if a variable (Y2) has a time lag with respect to Rsd, 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 the loop, marked by an arrow at the morning rising limb in Fig. 1b, indicates if a variable is preceding or lagging Rsd 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 Rsd. A negative phase lag appears as a clockwise (CW) loop.

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 Rsd(t) and its first-order time derivative dRsd(t)∕dt:

(2) Y ( t ) = a + b R sd ( t ) + c d R sd ( t ) / d t + ε ( t ) .

Using multilinear regression, we obtain the coefficients a, b, and c assuming a normal distribution of the residuals ε(t). If Y is linear with Rsd, 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 nonlinear term (with c≠0) would improve the model fit.

Figure 2Location of the EC site at Petit-Nobressart, Luxembourg. (a) shows the location within Western Europe. (b) shows an orthophoto of the surroundings of the system (ESRI® World Imagery). (c) shows a picture of the mast with micrometeorological sensors. The soil sensors are located on the right of the solar panel. Photo: Elisabeth Thiem.


Although significance testing of the coefficient c is an advantage, it is clear from Eq. (2) that the magnitude of c depends on the units and magnitude of the response variable Y. In order to estimate a comparable estimate of the phase lag we employ a harmonic transformation of the regression model. Assuming that Rsd is a harmonic function with an angular frequency ω, the phase difference φ can be estimated from the two regression coefficients b and c:

(3) φ = tan - 1 ( c ω / b ) .

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 Rsd as a time lag tφ (min) we use

(4) t φ = tan - 1 ( 48 / ( 2 π ) c / b ) ( 60 × 24 / ( 2 π ) ) .

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 Rsd is at noon local time, both metrics are qualitatively comparable.

2.2 Field site and observations

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) was installed at the grassland close to the village of Petit-Nobressart (Fig. 2; exact coordinates: 4946.77 N, 0548.22 E). The EC station was operated from 11 June until 23 July 2015. The three-dimensional wind and temperature fluctuations were measured at 2.41 m above ground by a sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, USA) facing to the mean wind direction of 290. A fast-response open-path CO2∕H2O infrared gas analyzer (IRGA LI-7500, LI-COR, USA) installed at a lateral 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).

Downwelling and upwelling shortwave and longwave radiation were obtained by a four-component net-radiation sensor (NR01, Hukseflux, the Netherlands). The meteorological variables (air temperature, humidity, and precipitation) were monitored with a time resolution of 30 min. Soil heat flux was measured by heat flux plates (two in 8 cm depth; HFP01, Hukseflux, the Netherlands), soil temperature was measured at 2, 5, 15, 30 cm depth (model 107, Campbell Scientific Inc., UK), water content at 2.5, 15, 30 cm depth (CS616, Campbell Scientific Inc., UK), and matric potential at 5, 15, 30 cm depth (model 253, Campbell Scientific Inc., UK). All soil sensors were installed between the turbulence and radiation measurement devices.

Unfortunately, the two upper-temperature probes and soil-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 to 30 June 2015. To still obtain soil heat fluxes for the entire measuring period, additionally harmonic wave analysis (Duchon and Hale, 2012) of the heat flux plate data was applied. The harmonic wave analysis calculates the wave spectrum at the soil surface from the Fourier transform of the soil heat flux measured by the heat flux plates in a few centimeter depth (here: 8 cm) by correcting for wave amplitude damping and phase shift. The surface ground heat flux is then obtained by an inverse Fourier transformation of the corrected wave spectrum. The method has a dependence on soil moisture affecting the damping depth. The dependence is, however, weak for clayey soils with soil water contents >10 % (Jury and Horton, 2004) as observed at the site. The damping depth was obtained by the exponential decay of the soil temperature amplitude measured at the various depths. Differences in the damping depth between wet and drier soil moisture conditions only yielded differences in G smaller than 10 W m−2. Therefore, we used a constant damping depth for the whole period.

Both methods for deriving the total soil heat flux agreed well for the period before 30 June, so 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 in this work. For more details on instrumentation and EC data processing, see Ingwersen et al. (2011) and Wizemann et al. (2015).

2.2.1 Derived meteorological variables

We derived the saturated water vapor pressure es (hPa) using the empirical Magnus equation (Magnus, 1844) as a function of air temperature T (C) with empirical coefficients from Alduchov and Eskridge (1996):


Then, the water vapor pressure of the air ea (hPa) was obtained by using air temperature Ta and relative humidity (RH):


To assess the moisture conditions of each date of the site we used the evaporative fraction fE:


Since daily averages can be influenced by single large values of the turbulent fluxes and contain missing values, we estimated a daily fE based on the 30 min values of each day using the following linear regression:


where fE 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 fE.

Table 1Variables provided by the surface energy balance station and used for this work.

Download Print Version | Download XLSX

Since the sonic anemometer measures friction velocity (u*) and the absolute value of wind speed u=(U2+V2), we estimate the aerodynamic conductance for momentum (u*2/u) and the aerodynamic conductance (gah) for heat including the excess resistance to heat transfer using an empirical formula by Thom (1972):

(5) g ah , Thom = u u * 2 + 6.2 u * 2 3 - 1 .

We chose to use this formula for its simplicity and similar performance than more recent, complex parameterizations (Knauer et al., 2018; Mallick et al., 2016). Also note that effects of atmospheric stability are accounted for in the first term of Eq. (5).

2.2.2 Energy balance closure gap correction

Most EC measurements show that the sum of the observed turbulent heat fluxes is smaller than the available energy and thus does not close the energy balance, leaving an energy balance closure gap (Qgap) (Foken, 2008):


For our site we observed on average a slope of (H+λE)(Rn-G)=0.81 (by linear regression) with an average gap of 37 W m−2 over the whole duration of the field campaign. These values are in the typical range of what is commonly found for grassland sites (Stoy et al., 2013).

To correct the turbulent fluxes for the energy balance closure gap (evaluated at the 30 min time steps), we use a correction based on the Bowen ratio (BR) (Twine et al., 2000), which is directly related to the evaporative fraction fE=1/(BR+1) to obtain corrected fluxes:




The correction is applied at 30 min time steps using the daily fE estimates. We use these corrected fluxes in the further analysis.

2.2.3 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 (Rsd,pot), which is a function of latitude and day of year (using R package REddyProc, function fCalcPotRadiation):


where t corresponds to each time step of measurement and with fdiff=0.78 being a constant factor taking into account atmospheric extinction of solar radiation.

2.3 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 (Ts) information as a key boundary condition (Kustas and Norman, 1999). 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

(6) λ E = R n - G - H = R n - G - ρ c p T s - T a g ah ,

where ρ is the density of air, cp is the specific heat of air at constant pressure, and gah is the effective aerodynamic conductance of heat that characterizes the transport of sensible heat between the surface and the atmosphere. We obtained Ts from the observed longwave emission of the surface Ts=(Rlu/(σεs))1/4 with σ=5.67×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 treat the surface as a uniform layer, two-source energy balance models partition temperatures as well as radiative and energy fluxes into a soil and vegetation component. The one-source approach (OSEB) parameterizes the aerodynamic conductance gah as follows (e.g., Kalma et al., 2008; Tang et al., 2013):


where zu and zt are the measurement heights of wind and air temperature, respectively; z0m and z0h 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 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 z0m were calculated as fractions of the vegetation height, hc, with d=0.65hc and z0m=0.125hc. The roughness length for heat z0h was set using the dimensionless parameter kB-1=ln(z0m/z0h), 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 conductance, as it is done in two-source approaches.

In addition to OSEB we applied the two-source energy balance model developed by Norman et al. (1995) and 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 equations are solved iteratively. It starts by assuming that a fraction of the canopy (described by vegetation greenness fraction fg) transpires at a potential rate as described by the Priestley–Taylor equation (Priestley and Taylor, 1972):

(8) λ E PT = α PT s s + γ R n - G ,

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 λEc=fgλEPT 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 this study we used a constant vegetation fraction of fc=0.9 and a greenness fraction fg which was derived from close-up pictures taken at the beginning and the end of the field campaign and linearly interpolated in-between.

Table 2Input variables used in the different evapotranspiration schemes.

Download Print Version | Download XLSX

2.4 Penman–Monteith approach

In the Penman–Monteith approach (Monteith, 1965) the inclusion of physiological conductance (gs) imposes a critical control on λE:

(9) λ E = s R n - G + ρ c p g av e s T a - e a s + γ 1 + g av g s .

In Eq. (9), the transfer of moisture is linked to a supply–demand reaction where the net available energy (RnG) is the supply energy for evaporation and the vapor pressure deficit of the air Da [=es(Ta)-ea] is the demand for evaporation from the atmosphere. In the PM approach, the two conductances, the aerodynamic conductance gav and the surface conductance gs, to water vapor are unknown. A widely used approach to obtain a reference evapotranspiration estimate from 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. 7) with wind speed u measured at 2 m above the surface, d=2/3h, z0m=0.123h, z0h=0.1z0m, yielding gav=u/208 (Box 4 in Allen et al., 1998). Surface conductance is fixed at a constant gs=1/70 m s−1. Here, we use the latter definitions of the conductances and use direct measurements for the other input variables of Eq. (9) to obtain the FAO Penman–Monteith estimate. While the FAO estimate is typically intended for estimates of the reference evaporation for well-watered grass on a daily basis, we use it here as a reference for comparison on a sub-daily 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. (5) using observations of friction velocity and wind speed, but keeping gs 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 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, 2018). STIC is a one-dimensional physically based surface energy balance 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 ga and gs on soil moisture through Ts and on environmental variables through Ta, Da, and net radiation. Therefore, 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. (2016, 2018) and Bhattarai et al. (2018).

Figure 3Daily observations of soil moisture, evaporative fraction, ratio of observed to potential solar radiation, and mean precipitation. (a) shows the daily time series and (b) the relationship of fE to soil moisture used to classify wet and dry days depending on fE>0.6 or fE<0.5, respectively. Sunny days are defined using a threshold of 85 % of Rsd to potential radiation and are marked with solid symbols, with blue circles referring to wet days and red squares to dry days. Top soil moisture measured at 5 cm below surface is shown.


3 Results

3.1 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 (1981–2010) the precipitation deficit in this region was −44 % in June and −41 % in July, respectively (source: meteorological 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 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 fE with time (Fig. 3a). Based on the observed fE we classified dry days with fE<0.5 and wet days with fE>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.

Figure 4Observations of average diurnal cycles of energy fluxes (a), with Rsd representing the shortwave downwelling flux, Rld the longwave downwelling flux, Rn the net radiation, G the ground heat flux, Ts and Ta the surface and air temperatures, and ea the air vapor pressure, comparing wet and dry days (b).


Based on the classification of 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 (Rsd). Figure 4a confirms that Rsd and net radiation (Rn) had very similar diurnal cycles and magnitudes for the wet and dry days. However, the downwelling longwave radiation Rld and the soil heat flux were somewhat higher under wet conditions (Fig. 4a). The higher Rld is related to higher air temperatures and air vapor pressures observed under wet conditions (Fig. 4b), which may explain the greater value of Rld 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 for shortwave radiation under wet and dry conditions, we observed a small but significant difference in the longwave radiative exchange.

Table 3Statistics for all days and sunny wet or dry days based on 30 min values during daytime hours 06:00–18:00 LT. Performance statistics, root mean square error (RMSE), and explained variance r2 are computed with respect to the observed latent heat flux corrected for the closure gap by the Bowen ratio (λEBRC). As a reference we also provide statistics for the uncorrected, observed latent heat flux (λEuncor). Potential evapotranspiration estimates are Priestley–Taylor (PT) and FAO Penman–Monteith (FAO-PM) reference evapotranspiration. Actual λE estimates are provided by the three schemes. Statistics are computed for all days and for clear-sky days classified as wet and dry.

Download Print Version | Download XLSX

Figure 5Average diurnal cycle of λE estimates for (a) dry and (b) wet days. Error bars denote the standard deviation obtained for each hour. The bold black line with squares shows the observed latent heat flux corrected for the surface energy balance closure (λEBRC). The grey-shaded area depicts the range induced by the energy balance closure gap.


3.2 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 Rsd forcing. Figure 5 shows the average diurnal cycle of observations and models for λE. The observations showed a significant difference in λE between dry and wet conditions, with a 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-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 OSEB and TSEB showed a tendency to overestimate λE under dry conditions but captured the high λE values under wet conditions. In contrast, STIC captured the low λE magnitude under dry conditions (fE<0.5) but underestimated λE under wet conditions (for fE>0.6).

Figure 6Diurnal hysteresis of λE to Rsd for (a) dry and (b) wet conditions of observations and different models. Bold arrows indicate the rising limb in the morning hours (07:00 to 08:00 LT) showing a counterclockwise hysteresis of λE under wet conditions. Vertical arrows depict the standard deviation of λEBRC for each hour.


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 and 259 W m−2) of the observed λ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 λ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 (r2=0.72) as compared to the other two models (r2=0.84–0.85). Therefore, we find that the correlation of the schemes is distinctly larger under wet conditions as compared to dry conditions. The correlations under wet conditions of OSEB and TSEB are in the range of the correlation of the uncorrected λE (r2=0.91), whereas STIC and FAO-PM showed lower correlation. Under dry conditions the correlation was significantly lower than the correlation of the uncorrected λE (r2=0.93). While OSEB–TSEB explained 62 % of the observed uncorrected λE variability in dry conditions (STIC explained 44 %), both models produced higher RMSE (57–58 W m−2) as compared to STIC (45 W m−2) under these conditions.

3.3 Diurnal patterns of evapotranspiration

The evaluation of the diurnal cycle shows that λE was strongly related to the incoming solar radiation, emphasizing that Rsd 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 Bowen-ratio-corrected flux (λEBRC) with a mean tφ=15 min under wet conditions and no significant lag under dry conditions (Fig. 7 and Table 4). The uncorrected observations showed only a slightly lower wet–dry difference, highlighting that the method to close the energy balance closure gap does not significantly influence the estimated phase lag.

Figure 7Boxplot of the daily phase lag of λE to Rsd for observed (without and with Bowen ratio correction) and modeled latent heat flux using sunny dry (red) and wet (blue) days. A positive phase lag means that λE follows Rsd and thus forms a CCW hysteresis as shown in Fig. 6. Dots show the actual data for each day with filled symbols indicating significant phase lags (P<0.05, t test of coefficient significantly different from 0).


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 gs in the PM equation but with gav obtained from Eq. (5) 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- and vapor-gradient scheme (STIC) showed relatively larger phase lags under both dry and wet conditions (tφ=14–20 min) (Fig. 7, Table 4).

Table 4Results of the Camuffo–Bernardi regression model with mean (standard deviation) for wet and dry days. The slope of the variable against Rsd is represented by b (note that the unit of b depends on the unit of the variable) and the phase lag to incoming solar radiation is converted to minutes. The adjusted explained variance by the multilinear regression model is given in column r2. The phase lag to RnG is reported in the last column for comparison.

Download Print Version | Download XLSX

Figure 8Diurnal patterns of observed surface energy balance components for (a) dry and (b) wet days. The lines show the composite average and vertical bars the standard deviation for available energy (black), latent heat (blue), and sensible heat (red) of each hour. There is a nearly linear response of all surface heat fluxes to Rsd under dry conditions and a systematic hysteresis loop under wet conditions. Under wet conditions the CCW hysteresis of λE is mostly compensated for by a CW hysteresis of H.


Since all evapotranspiration schemes use RnG as forcing, we also computed the phase lags with RnG as a reference variable (see Table 4). The differences to Rsd as reference are, however, rather small with slightly lower phase lags and in the range of the standard deviation of the daily estimates. This small difference can be attributed to a negligible phase lag between Rsd and Rn as well as the rather small magnitude and the phase lag of the soil heat flux.

Figure 9Diurnal patterns of observed anomalies in surface temperature (Ts), air temperature at 2 m (Ta), and their gradient (TsTa) for (a) dry and (b) wet days. Both Ta and Ts show a pronounced CCW hysteresis, but the form of the hysteretic loop is 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 (see Fig. 8).


3.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=(Rn-G)-H] with respect to Rsd (Fig. 8). Generally, there was only a small hysteresis in the available energy (RnG) (Table 4). The turbulent heat fluxes showed significant hysteresis under wet conditions but not under dry conditions. Interestingly, under wet conditions the CCW hysteresis of λE with a phase lag (mean tφ=15 min) was mostly compensated for by a CW hysteresis of H (mean tφ=-22 min) (Fig. 8 and Table 4). This compensation is an outcome of net available energy (RnG) showing little hysteresis for both conditions.

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 Rsd (Fig. 9). Air temperature showed a more pronounced hysteretic loop than surface temperature, and with a triangular shape with higher values during the afternoon when solar radiation reduces. Interestingly, the surface-to-air temperature gradient, being the driving gradient for the sensible heat flux, showed much lower 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 to H (see Table 4).

We further analyzed different formulations of the aerodynamic conductance (ga) 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 (gam=u*2/u). We then used the empirical formula by Thom (1972) to calculate the aerodynamic conductance for heat including the excess resistance to heat transfer (Eq. 5). Thirdly, we inferred the aerodynamic conductance from the observed sensible heat flux (HBRC) and temperature gradient (TsTa) by inverting HBRC using Eq. (6). The FAO-PM describes the aerodynamic conductance with a simple linear relationship to wind speed. OSEB and TSEB estimates the aerodynamic conductance to heat (gah), while STIC estimates the conductance to water vapor (gav). Thus by comparing these different conductance estimates we assume 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, 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, and they also use wind speed but rely on Monin–Obukhov similarity theory and stability correction. STIC, which does not use wind speed, did not show any significant differences in gav between wet and dry conditions.

Figure 10Boxplot of the different estimates of aerodynamic conductance under dry (red) and wet (blue) conditions. Only sunny days are sampled and midday values (10:00–15:00 LT) are used in the comparison. The top three estimates are directly inferred from observations, as described in the text.


Finally we analyze the diurnal patterns of the vapor pressure deficit Da=es(Ta)-ea, which is a critical driver of the latent heat flux in the PM equation. Since Da is derived from the observations, we analyzed its diurnal patterns in Fig. 11. We found that the vapor pressure in the air remained fairly constant during the day; hence it did not co-vary with Rsd and only showed a small CW hysteresis with higher vapor pressure during the morning compared to 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 Rsd, which is consistent with the large hysteresis in air temperature (Figs. 9 and 12). As a consequence, Da 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 (see Fig. 9). Furthermore, we note that the phase lag in Da 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 identified here may relate to a too-pronounced role of Da in the evapotranspiration estimate.

Figure 11Diurnal patterns of vapor pressure in air (black), the saturated vapor pressure evaluated at observed air temperature (red), and the vapor pressure deficit (blue) for (a) dry and (b) wet days.


4 Discussion

4.1 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 nonlinearity in the form of a phase lag between λE and Rsd was detected, which showed larger λE for the same Rsd 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 heat flux) revealed very small phase lags with Rsd. 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 the comprehensive 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 fE>0.35, while for very dry conditions (fE<0.2) a shift of the λE centroid towards the morning was found. This indicates that our results are not just applicable to Luxembourg, but are a general phenomenon which justifies a wider interpretation within temperate climates.

Figure 12Phase lag to solar radiation of surface energy fluxes and surface state variables used as input for the evapotranspiration models for dry (red) and wet (blue) days. Boxplot and daily estimates with filled symbols showing significant phase lag estimates.


It is important to emphasize here that the observed phase lags are not dominated by diurnal heat storage changes below the surface, since both the diurnal magnitude and the phase lag of the soil heat flux were relatively small compared to the turbulent heat fluxes. The models employed here use available energy (RnG) as input to estimate λE. However, the phase lag of the latent heat flux would only reduce by about 3 min when choosing RnG instead of Rsd as the reference variable to calculate the phase lags. Hence, the observed phase lags of λE and H to Rsd are not an artifact of the analysis, but can be considered as an imprint of L–A interaction.

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. Over land it is the lower atmosphere which acts as efficient heat storage to buffer most of the diurnal imbalance caused by solar radiation, because the heat storage of the subsurface is limited by the relatively slow heat conduction into the soil. Thus, 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.5 h. This large phase lag of air temperature then shapes (i) the vertical surface-to-air temperature gradient, which drives the sensible heat flux; 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.

We also found that the phase lag of the turbulent heat fluxes is affected by soil water availability. This is most clearly seen for the surface-to-air temperature gradient and the sensible heat flux, whose phase lag is 2 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 soil heat flux.

Our 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 Da 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 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.

4.2 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 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 the water vapor pressure deficit. Our findings are in line with observations of a near-linear relationship of λE to Rsd by Jackson et al. (1983) which stimulated remote-sensing-based spatial mapping of λE (Crago, 1996). Also the semi-empirical Makkink equation to estimate potential evapotranspiration (Makkink, 1957; de Bruin and Lablans, 1998; de Bruin et al., 2015) uses Rsd as the main driver.

Further support of the argument is given by the successful application of equilibrium evapotranspiration (Schmidt, 1915; Priestley and Taylor, 1972; Miralles et al., 2011; Renner et al., 2016) which uses Rn and air temperature as key inputs.

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

Our analysis allows for the better understanding of the relevance of the feedbacks which occur at a sub-daily timescale. 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; Jarvis and McNaughton, 1986). De Bruin and Holtslag (1982) showed that a positive correlation between Da and Rn allows simplifying the complex PM equation to a form similar to equilibrium evaporation (Eq. 8) with net radiation as the dominant driver. Therefore, simpler, energy-based formulations for λE show similar performance to PM-based approaches, but with less input parameters (De Bruin and Holtslag, 1982; Beljaars and Bosveld, 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 originates from the diurnal redistribution of heat, a mismatch would then clearly be seen at the sub-daily timescale. Hence by focusing on the internal relationship of the modeled λE to Rsd at the sub-daily timescale 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 show a reduction of the phase lag under dry conditions. The PM approaches use the vapor pressure deficit as an input which showed a substantial hysteresis loop on the order of 2.5 h lagging Rsd. This is due to the temperature dependency of the saturation vapor pressure, while the actual vapor pressure shows no relationship with Rsd. Besides Da and s, all other input variables to the Penman–Monteith approaches used here (both FAO and STIC) showed minor phase lags with respect to Rsd. Since the surface conductance in the FAO Penman–Monteith formulation is fixed with time, the resulting prediction of potential λE showed a significant and large phase lag on 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 gs, STIC computes λE through analytical estimation of gs and gav from the information of both the surface-to-air temperature gradient and the vapor pressure deficit. This dynamic treatment of gs reduced the phase lag to values similar to the observations under the wet conditions. However, under dry conditions STIC still showed significant phase lags, which may be related to the lag of Da to Rsd which was similar for both dry and wet conditions (Fig. 12). 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 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 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 (TsTa) which showed a significant counterclockwise, i.e., a leading, hysteretic loop, which is on the order of the phase shift detected for the sensible heat flux (Fig. 12). In addition, there is a distinct and significant increase in 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 (TsTa) to Rsd) and the diurnal phase lag (see Table 4).

While the PM approaches must identify two conductances simultaneously, the temperature-gradient approaches only need to parameterize the aerodynamic conductance to heat (gah) using wind speed as input. Therefore, we found that these approaches agreed well with the approximated gah from the EC tower, which shows an enhanced conductance under dry conditions. In contrast, the diagnosed gav 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 conductance parameters from the available radiometric and meteorological data (Mallick et al., 2018) for the climatic conditions in which these were evaluated here.

Note that we evaluated a temperate grassland site which experienced an exceptional summer drought. Therefore, 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 (fE<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, temperature-gradient approaches can show larger biases under semi-arid conditions (Morillas et al., 2013). The difficulty of temperature-gradient approaches is 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 timescale for the problem of surface conductance parameterizations was already highlighted by Matheny et al. (2014). However, they and others evaluated the diurnal patterns of the hysteretic loops between λE and Da (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 timescale we believe that solar radiation is better suited as a reference variable than Da. 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. These biases may well be compensated for at a longer timescales (Matheny et al., 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.

5 Conclusions

We analyzed the relationship of surface heat fluxes and states to incoming solar radiation at the sub-daily timescale for a temperate grassland site which experienced a summer drought. Most variables showed significant hysteresis loops which we objectively quantified by a linear component and a nonlinear phase lag component using multiple linear regression and harmonic analysis. We then compared these diurnal signatures obtained from observations of an eddy-covariance station with commonly used but structurally different approaches to model actual and potential evapotranspiration. The models have been forced by the observational data such that the differences to observations can be attributed to model formulation and signals contained in the input data. Our analysis guides model selection with a preference for the temperature-gradient approaches, because the vertical temperature gradient contains relevant signals of soil moisture limitation as opposed to the vapor pressure deficit of the air. 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 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 from observational data (e.g., remote sensing data) in climates without substantial water stress. Apparently, the 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 and surface 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., 2016; Konings et al., 2017).

We suggest that an evaluation of these schemes should be based on the sub-daily timescale, 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 detected at the aggregated timescales of days and months, it may lead to biased model sensitivities (Matheny et al., 2014). For example, an overly sensitive formulation of evapotranspiration to vapor pressure deficit and thus to temperature would predict larger changes in potential evapotranspiration under global warming (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 a 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 in representing the strong link of turbulent heat fluxes to solar radiation (Best et al., 2015; Haughton et al., 2016). Our findings provide an explanation of this model deficiency 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 temperatures. Correctly representing these metrics will lead towards a more accurate representation of the diurnal heat and mass exchange of the land with the atmosphere.

Code availability

The data analysis was performed with the open-source environment R (R Core Team, 2015). Functions to calculate the phase lag are provided as R package “phaselag” (Renner, 2019a), which is available from GitHub at, 2019. The script to perform data analysis and figures is published on Zenodo (Renner, 2019b) and can also be obtained from, 2019. Code to perform OSEB and TSEB simulations is published as Brenner (2019) and can also be found at, 2019. Code to simulate STIC1.2 simulations is available upon request from Kaniska Mallick (LIST,

Data availability

Data of observations (Wizemann et al., 2018, and model output (Renner et al., 2018, used in this study can be obtained from the research data repository GFZ Data Services (

Author contributions

MR and AK conceived the analysis of the diurnal cycle. VW, KS, and IT designed the field campaign. HW carried out the EC measurements and EC data processing. CB performed OSEB–TSEB simulations. KM performed STIC1.2 simulations. MR merged the data and performed the 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.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Linking landscape organisation and hydrological functioning: from hypotheses and observations to concepts, models and understanding (HESS/ESSD inter-journal SI)”. It is not associated with a conference.


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 sub-project “Understanding and characterizing land-surface–atmosphere exchange and feedbacks” (project number 182331427). Maik Renner and Axel Kleidon were funded by DFG grant number KL 2168/2-1. Claire Brenner was supported by the Austrian Science Fund (FWF) through funding of the CAOS Research Unit (I 2142-N29). Kaniska Mallick was supported by the Luxembourg Institute of Science and Technology (LIST) through the project BIOTRANS (grant number 00001145), CAOS-2 project grant INTER/DFG/14/02 funded by Fonds National de la Recherche (FNR)–DFG, and HiWET project funded by the Belgian Science Policy (BELSPO)–FNR under the program STEREO III (INTER/STEREOIII/13/03/HiWET; contract no. SR/00/301). IT was supported by the Luxembourg Institute of Science and Technology (LIST) through the project BIOTRANS (grant number 00001145).

The article processing charges for this open-access
publication were covered by the Max Planck Society.

Edited by: Shraddhanand Shukla
Reviewed by: Chunlüe Zhou and one anonymous referee


Alduchov, O. A. and Eskridge, R. E.: Improved Magnus Form Approximation of Saturation Vapor Pressure, J. Appl. Meteor., 35, 601–609,<0601:IMFAOS>2.0.CO;2, 1996. 

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspiration – Guidelines for computing crop water requirements, FAO Irrigation and drainage paper 56, FAO, Rome, 300, 6541, 1998. 

Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., and Holtslag, A. A. M.: A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation, J. Hydrol., 212–213, 198–212,, 1998. 

Beljaars, A. C. M. and Bosveld, F. C.: Cabauw Data for the Validation of Land Surface Parameterization Schemes, J. Climate, 10, 1172–1193,<1172:CDFTVO>2.0.CO;2, 1997. 

Best, M. J., Abramowitz, G., Johnson, H. R., Pitman, A. J., Balsamo, G., Boone, A., Cuntz, M., Decharme, B., Dirmeyer, P. A., Dong, J., Ek, M., Guo, Z., Haverd, V., van den Hurk, B. J. J., Nearing, G. S., Pak, B., Peters-Lidard, C., Santanello, J. A., Stevens, L., and Vuichard, N.: The Plumbing of Land Surface Models: Benchmarking Model Performance, J. Hydrometeorol., 16, 1425–1442,, 2015. 

Betts, A. K.: FIFE atmospheric boundary layer budget methods, J. Geophys. Res.-Atmos., 97, 18523–18531, 1992. 

Bhattarai, N., Mallick, K., Brunsell, N. A., Sun, G., and Jain, M.: Regional evapotranspiration from an image-based implementation of the Surface Temperature Initiated Closure (STIC1.2) model and its validation across an aridity gradient in the conterminous US, Hydrol. Earth Syst. Sci., 22, 2311–2341,, 2018. 

Brenner, C.: ClaireBrenner/pyTSEB_Renner_et_al2018: OSEB/TSEB code for Renner et al. 2019, Zenodo,, 2019. 

Brenner, C., Thiem, C. E., Wizemann, H.-D., Bernhardt, M., and Schulz, K.: Estimating spatially distributed turbulent heat fluxes from high-resolution thermal imagery acquired with a UAV system, Int. J. Remote Sens., 38, 3003–3026,, 2017. 

Brutsaert, W.: Evaporation into the Atmosphere, Springer Netherlands, Dordrecht, 1982. 

Brutsaert, W.: Hydrology: an introduction, Cambridge University Press, Cambridge, 2005. 

Brutsaert, W. and Stricker, H.: An Advection-Aridity Approach to Estimate Actual Regional Evapotranspiration, Water Resour. Res., 15, 443–450,, 1979. 

Camuffo, D. and Bernardi, A.: An observational study of heat fluxes and their relationships with net radiation, Bound.-Lay. Meteorol., 23, 359–368,, 1982. 

Choi, M., Kustas, W. P., Anderson, M. C., Allen, R. G., Li, F., and Kjaersgaard, J. H.: An intercomparison of three remote sensing-based surface energy balance algorithms over a corn and soybean production region (Iowa, US) during SMACEX, Agr. Forest Meteorol., 149, 2082–2097, 2009. 

Crago, R. D.: Conservation and variability of the evaporative fraction during the daytime, J. Hydrol., 180, 173–194,, 1996. 

De Bruin, H. A. R. and Holtslag, A. A. M.: A Simple Parameterization of the Surface Fluxes of Sensible and Latent Heat During Daytime Compared with the Penman–Monteith Concept, J. Appl. Meteorol., 21, 1610–1621,<1610:ASPOTS>2.0.CO;2, 1982. 

de Bruin, H. A. R. and Lablans, W. N.: Reference crop evapotranspiration determined with a modified Makkink equation, Hydrol. Process., 12, 1053–1062,<1053::AID-HYP639>3.0.CO;2-E, 1998. 

de Bruin, H. a. R., Trigo, I. F., Bosveld, F. C., and Meirink, J. F.: A Thermodynamically Based Model for Actual Evapotranspiration of an Extensive Grass Field Close to FAO Reference, Suitable for Remote Sensing Application, J. Hydrometeorol., 17, 1373–1382,, 2015. 

Duan, Z. and Bastiaanssen, W. G. M.: Evaluation of three energy balance-based evaporation models for estimating monthly evaporation for five lakes using derived heat storage changes from a hysteresis model, Environ. Res. Lett., 12, 024005,, 2017. 

Duchon, C. and Hale, R.: Time series analysis in meteorology and climatology: an introduction, in: Vol. 7, John Wiley & Sons, 2012. 

Foken, T.: The energy balance closure problem: An overview, Ecol. Appl., 18, 1351–1367,, 2008. 

French, A. N., Hunsaker, D. J., and Thorp, K. R.: Remote sensing of evapotranspiration over cotton using the TSEB and METRIC energy balance models, Remote Sens. Environ., 158, 281–294, 2015. 

Fuchs, M. and Hadas, A.: The heat flux density in a non-homogeneous bare loessial soil, Bound.-Lay. Meteorol., 3, 191–200,, 1972. 

Haughton, N., Abramowitz, G., Pitman, A. J., Or, D., Best, M. J., Johnson, H. R., Balsamo, G., Boone, A., Cuntz, M., Decharme, B., Dirmeyer, P. A., Dong, J., Ek, M., Guo, Z., Haverd, V., van den Hurk, B. J. J., Nearing, G. S., Pak, B., Santanello, J. A., Stevens, L. E., and Vuichard, N.: The Plumbing of Land Surface Models: Is Poor Performance a Result of Methodology or Data Quality?, J. Hydrometeorol., 17, 1705–1723,, 2016. 

Ingwersen, J., Steffens, K., Högy, P., Warrach-Sagi, K., Zhunusbayeva, D., Poltoradnev, M., Gäbler, R., Wizemann, H.-D., Fangmeier, A., Wulfmeyer, V., and Streck, T.: Comparison of Noah simulations with eddy covariance and soil water measurements at a winter wheat stand, Agr. Forest Meteorol., 151, 345–355,, 2011. 

Jackson, R. D., Hatfield, J. L., Reginato, R. J., Idso, S. B., and Pinter, P. J.: Estimation of daily evapotranspiration from one time-of-day measurements, Agr. Water Manage., 7, 351–362, 1983. 

Jarvis, P. G.: The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field, Philos. T. Roy. Soc. Lond. B, 273, 593–610,, 1976. 

Jarvis, P. G. and McNaughton, K. G.: Stomatal Control of Transpiration: Scaling Up from Leaf to Region, in Advances in Ecological Research, in: Vol. 15, edited by: MacFadyen, A. and Ford, E. D., Academic Press, London, 1–49, 1986. 

Jiménez, P. A., Dudhia, J., González-Rouco, J. F., Navarro, J., Montávez, J. P., and García-Bustamante, E.: A revised scheme for the WRF surface layer formulation, Mon. Weather Rev., 140, 898–918, 2012. 

Jury, W. A. and Horton, R.: Soil Physics, 6th Edn., John Wiley & Sons, Inc, Hoboken, NJ, 2004. 

Kalma, J. D., McVicar, T. R., and McCabe, M. F.: Estimating Land Surface Evaporation: A Review of Methods Using Remotely Sensed Surface Temperature Data, Surv. Geophys., 29, 421–469,, 2008. 

Kleidon, A. and Renner, M.: An explanation for the different climate sensitivities of land and ocean surfaces based on the diurnal cycle, Earth Syst. Dynam., 8, 849–864,, 2017. 

Knauer, J., Zaehle, S., Medlyn, B. E., Reichstein, M., Williams, C. A., Migliavacca, M., De Kauwe, M. G., Werner, C., Keitel, C., and Kolari, P.: Towards physiologically meaningful water-use efficiency estimates from eddy covariance data, Global Change Biol., 24, 694–710, 2018. 

Konings, A. G., Williams, A. P., and Gentine, P.: Sensitivity of grassland productivity to aridity controlled by stomatal and xylem regulation, Nat. Geosci., 10, 284–288,, 2017. 

Kustas, W. P. and Norman, J. M.: Evaluation of soil and vegetation heat flux predictions using a simple two-source model with radiometric temperatures for partial canopy cover, Agr. Forest Meteorol., 94, 13–29,, 1999. 

Kustas, W. P., Nieto, H., Morillas, L., Anderson, M. C., Alfieri, J. G., Hipps, L. E., Villagarcía, L., Domingo, F., and Garcia, M.: Revisiting the paper “Using radiometric surface temperature for surface energy flux estimation in Mediterranean drylands from a two-source perspective”, Remote Sens. Environ., 184, 645–653,, 2016. 

Maes, W. H., Gentine, P., Verhoest, N. E. C., and Miralles, D. G.: Potential evaporation at eddy-covariance sites across the globe, Hydrol. Earth Syst. Sci. Discuss.,, in review, 2018. 

Magnus, G.: Versuche über die Spannkräfte des Wasserdampfs, Ann. Phys., 137, 225–247,, 1844. 

Makkink, G.: Testing the Penman formula by means of lysimeters, J. Inst. Water Eng., 11, 277–288, 1957. 

Mallick, K., Jarvis, A. J., Boegh, E., Fisher, J. B., Drewry, D. T., Tu, K. P., Hook, S. J., Hulley, G., Ardö, J., Beringer, J., Arain, A., and Niyogi, D.: A Surface Temperature Initiated Closure (STIC) for surface energy balance fluxes, Remote Sens. Environ., 141, 243–261,, 2014. 

Mallick, K., Boegh, E., Trebs, I., Alfieri, J. G., Kustas, W. P., Prueger, J. H., Niyogi, D., Das, N., Drewry, D. T., Hoffmann, L., and Jarvis, A. J.: Reintroducing radiometric surface temperature into the Penman–Monteith formulation, Water Resour. Res., 51, 6214–6243,, 2015. 

Mallick, K., Trebs, I., Boegh, E., Giustarini, L., Schlerf, M., Drewry, D. T., Hoffmann, L., von Randow, C., Kruijt, B., Araùjo, A., Saleska, S., Ehleringer, J. R., Domingues, T. F., Ometto, J. P. H. B., Nobre, A. D., de Moraes, O. L. L., Hayek, M., Munger, J. W., and Wofsy, S. C.: Canopy-scale biophysical controls of transpiration and evaporation in the Amazon Basin, Hydrol. Earth Syst. Sci., 20, 4237–4264,, 2016. 

Mallick, K., Toivonen, E., Trebs, I., Boegh, E., Cleverly, J., Eamus, D., Koivusalo, H., Drewry, D., Arndt, S. K., Griebel, A., Beringer, J., and Garcia, M.: Bridging Thermal Infrared Sensing and Physically-Based Evapotranspiration Modeling: From Theoretical Implementation to Validation Across an Aridity Gradient in Australian Ecosystems, Water Resour. Res., 54, 3409–3435,, 2018. 

Massman, W. J.: Correcting errors associated with soil heat flux measurements and estimating soil thermal properties from soil temperature and heat flux plate data, Agr. Forest Meteorol., 59, 249–266,, 1992. 

Matheny, A. M., Bohrer, G., Stoy, P. C., Baker, I. T., Black, A. T., Desai, A. R., Dietze, M. C., Gough, C. M., Ivanov, V. Y., Jassal, R. S., Novick, K. A., Schäfer, K. V. R., and Verbeeck, H.: Characterizing the diurnal patterns of errors in the prediction of evapotranspiration by several land-surface models: An NACP analysis, J. Geophys. Res.-Biogeo., 119, 1458–1473,, 2014. 

Mauder, M. and Foken, T.: Eddy-Covariance Software TK3, University of Bayreuth, Bayreuth, 2015. 

Milly, P. C. D. and Dunne, K. A.: Potential evapotranspiration and continental drying, Nat. Clim. Change, 6, 946–949,, 2016. 

Miralles, D. G., Holmes, T. R. H., De Jeu, R. A. M., Gash, J. H., Meesters, A. G. C. A., and Dolman, A. J.: Global land-surface evaporation estimated from satellite-based observations, Hydrol. Earth Syst. Sci., 15, 453–469,, 2011. 

Monteith, J. L.: Evaporation and environment, Symp. Soc. Exp. Biol., 19, 205–234, 1965. 

Morillas, L., García, M., Nieto, H., Villagarcia, L., Sandholt, I., Gonzalez-Dugo, M. P., Zarco-Tejada, P. J., and Domingo, F.: Using radiometric surface temperature for surface energy flux estimation in Mediterranean drylands from a two-source perspective, Remote Sens. Environ., 136, 234–246, 2013. 

Mu, Q. Zhao, M., and Running, S. W.: Improvements to a MODIS global terrestrial evapotranspiration algorithm, Remote Sens. Environ., 115, 1781–1800,, 2011. 

Nelson, J. A., Carvalhais, N., Migliavacca, M., Reichstein, M., and Jung, M.: Water-stress-induced breakdown of carbon–water relations: indicators from diurnal FLUXNET patterns, Biogeosciences, 15, 2433–2447,, 2018. 

Niu, G.-Y., Yang, Z.-L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res., 116, D12109,, 2011. 

Norman, J. M., Kustas, W. P., and Humes, K. S.: Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature, Agr. Forest Meteorol., 77, 263–293,, 1995. 

Ohmura, A.: The development and the present status of energy balance climatology, J. Meteorol. Soc. Jpn. Ser. II, 92, 245–285, 2014. 

Oke, T. R.: Boundary layer climates, Routledge, Londan, New York, 1987. 

Penman, H. L.: Natural evaporation from open water, bare soil and grass, P. Roy. Soc. Lond. A, 193, 120–145, 1948. 

Prentice, I. C., Dong, N., Gleason, S. M., Maire, V., and Wright, I. J.: Balancing the costs of carbon gain and water transport: testing a new theoretical framework for plant functional ecology, Ecol. Lett., 17, 82–91,, 2014. 

Priestley, C. and Taylor, R.: On the assessment of surface heat flux and evaporation using large-scale parameters, Mon. Weather Rev., 100, 81–92, 1972. 

R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, available at: (last access: January 2019), 2015. 

Renner, M.: Phaselag: R package to compute the phase lag between two variables, Zenodo,, 2019a. 

Renner, M.: R script to reproduce phase lag analysis of the diurnal cycle of evapotranspiration, Zenodo,, 2019b. 

Renner, M., Hassler, S. K., Blume, T., Weiler, M., Hildebrandt, A., Guderle, M., Schymanski, S. J., and Kleidon, A.: Dominant controls of transpiration along a hillslope transect inferred from ecohydrological measurements and thermodynamic limits, Hydrol. Earth Syst. Sci., 20, 2063–2083,, 2016. 

Renner, M., Wizemann, H.-D., Brenner, C., Mallick, K., Trebs, I., Wulfmeyer, V., Schulz, K., and Kleidon, A.: Surface energy balance observations and model output at a grassland site in Luxembourg, GFZ Data Services,, 2018. 

Santanello, J. A. and Friedl, M. A.: Diurnal Covariation in Soil Heat Flux and Net Radiation, J. Appl. Meteorol., 42, 851–862,<0851:DCISHF>2.0.CO;2, 2003. 

Santanello, J. A., Peters-Lidard, C. D., Kumar, S. V., Alonge, C., and Tao, W.-K.: A Modeling and Observational Framework for Diagnosing Local Land–Atmosphere Coupling on Diurnal Time Scales, J. Hydrometeorol., 10, 577–599,, 2009. 

Santanello, J. A., Dirmeyer, P. A., Ferguson, C. R., Findell, K. L., Tawfik, A. B., Berg, A., Ek, M., Gentine, P., Guillod, B. P., van Heerwaarden, C., Roundy, J., and Wulfmeyer, V.: Land–Atmosphere Interactions: The LoCo Perspective, B. Am. Meteorol. Soc., 99, 1253–1272,, 2018. 

Schmidt, W.: Strahlung und Verdunstung an freien Wasserflächen; ein Beitrag zum Wärmehaushalt des Weltmeers und zum Wasserhaushalt der Erde, Ann. Calender Hydrographie und Maritimen Meteorologie, 43, 111–124, 1915. 

Stoy, P. C., Mauder, M., Foken, T., Marcolla, B., Boegh, E., Ibrom, A., Arain, M. A., Arneth, A., Aurela, M., Bernhofer, C., Cescatti, A., Dellwik, E., Duce, P., Gianelle, D., van Gorsel, E., Kiely, G., Knohl, A., Margolis, H., McCaughey, H., Merbold, L., Montagnani, L., Papale, D., Reichstein, M., Saunders, M., Serrano-Ortiz, P., Sottocornola, M., Spano, D., Vaccari, F., and Varlagin, A.: A data-driven analysis of energy balance closure across FLUXNET research sites: The role of landscape scale heterogeneity, Agr. Forest Meteorol., 171–172, 137–152,, 2013. 

Sun, T., Wang, Z.-H., and Ni, G.-H.: Revisiting the hysteresis effect in surface energy budgets, Geophys. Res. Lett., 40, 1741–1747,, 2013. 

Tang, R., Li, Z.-L., Jia, Y., Li, C., Chen, K.-S., Sun, X., and Lou, J.: Evaluating one- and two-source energy balance models in estimating surface evapotranspiration from Landsat-derived surface temperature and field measurements, Int. J. Remote Sens., 34, 3299–3313,, 2013. 

Thom, A. S.: Momentum, mass and heat exchange of vegetation, Q. J. Roy. Meteorol. Soc., 98, 124–134,, 1972. 

Timmermans, W. J., Kustas, W. P., Anderson, M. C., and French, A. N.: An intercomparison of the Surface Energy Balance Algorithm for Land (SEBAL) and the Two-Source Energy Balance (TSEB) modeling schemes, Remote Sens. Environ., 108, 369–384,, 2007. 

Twine, T. E., Kustas, W. P., Norman, J. M., Cook, D. R., Houser, P., Meyers, T. P., Prueger, J. H., Starks, P. J., and Wesely, M. L.: Correcting eddy-covariance flux underestimates over a grassland, Agr. Forest Meteorol., 103, 279–300, 2000. 

van Heerwaarden, C. C., Vilà-Guerau de Arellano, J., Gounou, A., Guichard, F., and Couvreux, F.: Understanding the Daily Cycle of Evapotranspiration: A Method to Quantify the Influence of Forcings and Feedbacks, J. Hydrometeorol., 11, 1405–1422,, 2010. 

Wilson, K. B., Baldocchi, D., Falge, E., Aubinet, M., Berbigier, P., Bernhofer, C., Dolman, H., Field, C., Goldstein, A., and Granier, A.: Diurnal centroid of ecosystem energy and carbon fluxes at FLUXNET sites, J. Geophys. Res.-Atmos., 108, 4664,, 2003. 

Wizemann, H.-D., Ingwersen, J., Högy, P., Warrach-Sagi, K., Streck, T., and Wulfmeyer, V.: Three year observations of water vapor and energy fluxes over agricultural crops in two regional climates of Southwest Germany, Meteorol. Z., 24, 39–59,, 2015. 

Wizemann, H.-D., Trebs, I., and Wulfmeyer, V.: Surface energy balance observations at a grassland site in Luxembourg, GFZ Data Services,, 2018. 

Wolf, A., Anderegg, W. R. L., and Pacala, S. W.: Optimal stomatal behavior with competition for water and risk of hydraulic impairment, P. Natl. Acad. Sci. USA, 113, E7222–E7230,, 2016.  

Wulfmeyer, V., Turner, D. D., Baker, B., Banta, R., Behrendt, A., Bonin, T., Brewer, W. A., Buban, M., Choukulkar, A., Dumas, E., Hardesty, R. M., Heus, T., Ingwersen, J., Lange, D., Lee, T. R., Metzendorf, S., Muppa, S. K., Meyers, T., Newsom, R., Osman, M., Raasch, S., Santanello, J., Senff, C., Späth, F., Wagner, T., and Weckwerth, T.: A New Research Approach for Observing and Characterizing Land–Atmosphere Feedback, B. Am. Meteorol. Soc., 99, 1639–1667,, 2018. 

Xiang, T., Vivoni, E. R., Gochis, D. J., and Mascaro, G.: On the diurnal cycle of surface energy fluxes in the North American monsoon region using the WRF-Hydro modeling system, J. Geophys. Res.-Atmos., 122, 9024–9049, 2017. 

Yang, Y., Long, D., Guan, H., Liang, W., Simmons, C., and Batelaan, O.: Comparison of three dual-source remote sensing evapotranspiration models during the MUSOEXE-12 campaign: Revisit of model physics, Water Resour. Res., 51, 3145–3165,, 2015. 

Zhang, Q., Manzoni, S., Katul, G., Porporato, A., and Yang, D.: The hysteretic evapotranspiration – Vapor pressure deficit relation: ET-VPD hysteresis, J. Geophys. Res.-Biogeo., 119, 125–140,, 2014. 

Zheng, H., Wang, Q., Zhu, X., Li, Y., and Yu, G.: Hysteresis Responses of Evapotranspiration to Meteorological Factors at a Diel Timescale: Patterns and Causes, PLOS ONE, 9, e98857,, 2014. 

Zhou, C. and Wang, K.: Evaluation of Surface Fluxes in ERA-Interim Using Flux Tower Data, J. Climate, 29, 1573–1582,, 2016. 

Zhou, C., Wang, K., and Ma, Q.: Evaluation of Eight Current Reanalyses in Simulating Land Surface Temperature from 1979 to 2003 in China, J. Climate, 30, 7379–7398,, 2017. 

Zhou, C., He, Y., and Wang, K.: On the suitability of current atmospheric reanalyses for regional warming studies over China, Atmos. Chem. Phys., 18, 8113–8136,, 2018. 

Short summary
We estimate the phase lag of surface states and heat fluxes to incoming solar radiation at the sub-daily timescale. While evapotranspiration reveals a minor phase lag, the vapor pressure deficit used as input by Penman–Monteith approaches shows a large phase lag. The surface-to-air temperature gradient used by energy balance residual approaches shows a small phase shift in agreement with the sensible heat flux and thus explains the better correlation of these models at the sub-daily timescale.