Using the maximum entropy production approach to integrate energy budget modelling in a hydrological model

Total terrestrial evaporation, also referred to as evapotranspiration, is a key process for understanding the hydrological impacts of climate change given that warmer surface temperatures translate into an increase in the atmospheric evaporative demand. To simulate this flux, many hydrological models rely on the concept of potential evaporation (PET), although large differences have been observed in the response of PET models to climate change. The maximum entropy production (MEP) model of land surface fluxes offers an alternative approach for simulating terrestrial evaporation in a simple way while fulfilling the physical constraint of energy budget closure and providing a distinct estimation of evaporation and transpiration. The objective of this work is to use the MEP model to integrate energy budget modelling within a hydrological model. We coupled the MEP model with HydroGeoSphere (HGS), an integrated surface and subsurface hydrologic model. As a proof of concept, we performed one-dimensional soil column simulations at three sites of the AmeriFlux network. The coupled model (HGS-MEP) produced realistic simulations of soil water content (root-mean-square error – RMSE – between 0.03 and 0.05 m³ m⁻³; NSE – Nash–Sutcliffe efficiency – between 0.30 and 0.92) and terrestrial evaporation (RMSE between 0.31 and 0.71 mm d⁻¹; NSE between 0.65 and 0.88) under semi-arid, Mediterranean and temperate climates. At the daily timescale, HGS-MEP outperformed the stand-alone HGS model where total terrestrial evaporation is derived from potential evaporation, which we computed using the Penman–Monteith equation, although both models had comparable performance at the half-hourly timescale. This research demonstrated the potential of the MEP model to improve the simulation of total terrestrial evaporation in hydrological models, including for hydrological projections under climate change.

1 Introduction

Driven by climate change, warmer surface temperatures are expected to increase the atmospheric evaporative demand (Breshears et al., 2013; Ficklin and Novick, 2017). As such, total terrestrial evaporation (hereinafter terrestrial evaporation), or “evapotranspiration”, is a key process for assessing the impacts of climate change on stream discharge. To simulate this flux, many hydrological models rely on the concept of potential evaporation (PET), which is defined as the maximum evaporation that can occur under ambient meteorological conditions with an unlimited water supply. Two types of PET models are generally used in hydrological models: temperature-based models where PET is estimated via air temperature (Thornthwaite, 1948; Hamon, 1963) and physically based models where PET estimation is based on components of the surface energy budget (Monteith, 1965; Priestley and Taylor, 1972). Various studies have compared PET models to assess their suitable range of application (Xu and Singh, 2001). While the performance of PET models varies across regions, the choice of a PET model appears to have minimal influence on stream discharge simulation under contemporary climate conditions (Andréassian et al., 2004; Oudin et al., 2005; Isabelle et al., 2018). However, when future climate projections are incorporated in PET models, the predicted PET varies significantly from model to model (McKenney and Rosenberg, 1993; Kingston et al., 2009; Donohue et al., 2010; Lofgren et al., 2011; McAfee, 2013; Hosseinzadehtalaei et al., 2016). These differences often translate into large uncertainty in stream discharge projections (Kay and Davies, 2008; Bae et al., 2011; Milly and Dunne, 2011; Prudhomme and Williamson, 2013; Seiller and Anctil, 2016), although some studies have demonstrated opposite results with a relatively low sensitivity of discharge projection to PET model selection (Thompson et al., 2014; Koedyk and Kingston, 2016).

Among PET models, those that are based on temperature tend to mainly reflect changes in mean air temperature, even though temperature is not necessarily the strongest control on PET (Donohue et al., 2010; Shaw and Riha, 2011). For example, the oversensitivity of temperature-based models to surface warming has been linked to an exaggerated assessment of drought severity under climate change (Hoerling et al., 2012). As a result, recent studies have recommended the use of physically based PET models, such as the Penman model, for simulating terrestrial evaporation in climate change assessments (Hobbins et al., 2008; Sheffield et al., 2012; McAfee, 2013). However, additional data inputs (e.g. radiation, humidity and wind speed) are required to implement such models, and, as opposed to the high confidence placed on temperature estimates from global climate models (McMahon et al., 2015), these variables usually suffer from a lower reliability. For example, large differences persist between observations, reanalysis products and global climate model projections of wind speed (McVicar et al., 2008; Pryor et al., 2009). Overall, a dilemma emerges for the simulation of PET under climate change, which involves choosing between a reliable physically based model with greater input data uncertainty versus an overly simplified temperature-based model with reliable input data (Ekström et al., 2007; Kay and Davies, 2008; Kingston et al., 2009).

Given these current limitations with PET models, alternative approaches are needed for the simulation of terrestrial evaporation in climate change assessments. Energy budget modelling offers a physically based approach for simulating terrestrial evaporation (along with other components of the surface energy budget), with the added benefit of conservation of energy that ensures a balance between incoming and outgoing energy. Land surface models provide a means to simulate both water and energy budgets, and they have been used for hydrological modelling, either by (i) coupling hydrological and land surface models (Pietroniro et al., 2001; Maxwell and Miller, 2005; Kunstmann et al., 2008; Zabel and Mauser, 2013) or by (ii) coupling a land surface model to a routing scheme (Gaborit et al., 2017). In land surface models, terrestrial evaporation is generally estimated with a bulk aerodynamic approach (Noilhan and Planton, 1989) that relies on semi-empirical equations and requires information on vertical gradients of air temperature and humidity that can introduce substantial errors. The use of land surface models can yield small gains in hydrological modelling performance (Livneh et al., 2011; Shi et al., 2014), but they remain computationally heavy.

Optimality principles, which are derived from the idea that nature organizes itself to ensure optimal functioning, provide an alternative avenue for modelling terrestrial evaporation. Maximum entropy production (MEP) is one of the optimality principles put forward, and although it has yet to become an established principle, it offers a promising method to improve hydrological modelling (Ehret et al., 2014; Westhoff and Zehe, 2013). The MEP principle has been applied in two ways in hydrology. In the first approach, MEP is defined as a physical principle, which hypothesizes that an open thermodynamic system far from equilibrium can achieve a steady state at which, through self-organization, entropy is produced at the maximum possible rate. Such a dynamic equilibrium can be achieved because, as demonstrated by Paltridge (1975), a gradient drives a flux that, in turn, depletes the initial gradient and further reinforces it. Using this approach, MEP has been used to constrain parameters of hydrological models (Kleidon and Schymanski, 2008; Porada et al., 2011; Westhoff and Zehe, 2013). In the second approach, MEP is defined as a statistical principle and constitutes the most probable state of an open system (Dewar, 2005, 2009). Using the MaxEnt statistical inference algorithm (Jaynes, 1957), the state of MEP can be predicted by maximizing Shannon information entropy while considering constraints imposed by the available information. Both concepts of entropy are linked: the maximization of Shannon information entropy can be used to assess the probability of a given state for any kind of system, and as such, thermodynamic entropy can be considered a special case of Shannon information theory. In hydrology, Wang and Bras (2009, 2011) have used the Shannon information entropy approach to develop the MEP model of land surface fluxes, which is the topic of the present article.

Table 1 Description of study sites.

The MEP model of land surface fluxes (Wang and Bras, 2009, 2011) offers a simple approach for energy budget modelling that ensures the closure of the surface energy balance. In the MEP model developed by Wang and Bras, only three input variables, net radiation, surface temperature and surface specific humidity, are required to model surface heat fluxes, with a different definition of the surface when assessing evaporation (surface is soil) and transpiration (surface is leaf). As it will be presented below, in practice, the following six input variables are needed to operate MEP continuously under all climate conditions: net radiation, soil surface temperature and specific humidity, leaf surface temperature and specific humidity, and a vegetation index. Overall, the MEP model eliminates the need for wind speed and surface roughness data (i.e. input data to the Penman model) as well as the need for data on vertical gradients (i.e. temperature and humidity, input data to the bulk aerodynamic model often implemented in land surface models). Moreover, the predicted surface heat fluxes are constrained by the available energy (i.e. net radiation), which avoids the issue of oversensitivity to temperature and suggests that this approach should be robust in climate change assessments. Recent studies also reported improved performance of the MEP approach for the simulation of terrestrial evaporation when compared against the Penman–Monteith model (Hajji et al., 2018) or a modified Penman–Monteith approach driven by remote sensing and reanalysis data (Xu et al., 2019). The MEP model also had a performance comparable to a complex land surface model (Canadian Land Surface Scheme – CLASS) at snow-free sites at low latitudes (Alves et al., 2019). Overall, the MEP model offers an attractive approach to improve the simulation of terrestrial evaporation in hydrological modelling and to increase the robustness of streamflow projections under climate change. The objective of this study is thus to couple the MEP model of land surface fluxes with a hydrological model and perform proof-of-concept simulations at AmeriFlux sites spanning a range of climatic and vegetation conditions.

2 Study area

We selected three sites of the AmeriFlux network (Baldocchi et al., 2001) where the eddy covariance method is used to measure vertical water and energy fluxes (Table 1; Fig. 1). We selected sites using the following criteria: (i) measurements of volumetric soil water content available at different depths and for an extended period, (ii) absence of snow cover, and (iii) diversity in climates (semi-arid, Mediterranean and temperate) and types of vegetation (grassland, woody savanna and deciduous broadleaf forest).

Figure 1 Location of study sites across the conterminous United States.

The first site, US-Wkg, is characterized by a semi-arid climate where 60 % of precipitation is concentrated in July, August and September. In 2006, a sharp transition occurred in the vegetation and the native grassland species (assemblage of grama grasses, genus Bouteloua) have been supplanted by Lehman lovegrass (Eragrostis lehmanniana), an exotic species spreading throughout the southwestern United States (Moran et al., 2009; Scott et al., 2010). The second site, US-Ton, is characterized by a Mediterranean climate, with dry summers and most precipitation occurring from October to May. The woody savanna is made up of two layers of vegetation, each reaching peak activity at different times of the year (Baldocchi et al., 2004). Trees, mainly blue oaks (Quercus douglasii), are dormant during winter months, reach peak activity in the spring and then carefully regulate their water use during the dry summer period. In contrast, the understory vegetation composed of grasses and forbs is mainly active during the rainy winter period, when water is plentiful. Finally, the US-WBW site is characterized by a temperate climate, with precipitation relatively evenly distributed throughout the year. The mixed deciduous forest is dominated by oak and hickory, and vegetation is active from the spring to early autumn (April to October), with peak activity during the summer (Wilson and Baldocchi, 2000).

3 Methodology

3.1 Maximum entropy production (MEP) model

The MEP model of land surface fluxes uses the optimality principle of maximum entropy production (Dewar, 2005) as an inference tool for nonequilibrium thermodynamic systems. In the MEP model, entropy refers to Shannon entropy, the expected value of information (Shannon, 1948), which is not related to thermodynamic entropy expressed as the ratio of flux to temperature. Indeed, the MEP model is derived from the principle of maximum entropy (MaxEnt) developed by Jaynes (1957) as a general inference tool for assigning probability distributions in statistical mechanics. Wang and Bras (2009) used this statistical approach to develop the MEP model of land surface fluxes. Within the MEP model, land surface fluxes are each expressed in terms of a dissipation or entropy function, and under the constraint of conservation of energy, a unique extremum solution is found.

Dissipation functions have been postulated for land surface fluxes over a dry soil (Wang and Bras, 2009) and over wet soil and vegetation surfaces (Wang and Bras, 2011). For non-vegetated surfaces, bare soil evaporation (E_s) is estimated by solving Eqs. (1) and (2) with the constraint of energy conservation at the surface ($R_{\text{n}}=G+E_{\text{s}}+H$; Wang and Bras, 2011):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{E}_{\text{s}}=B\left(\mathit{\sigma }\right)H,\text{(2)}& {\displaystyle }& {\displaystyle }G={\displaystyle \frac{B\left(\mathit{\sigma }\right)}{\mathit{\sigma }}}{\displaystyle \frac{I_{\text{s}}}{I_{\mathrm{0}}}}H{\left|H\right|}^{-\frac{\mathrm{1}}{\mathrm{6}}},\end{array}$$

where R_n, E_s, H and G are the net radiation and latent, sensible and ground heat fluxes at the soil surface (W m⁻²; positive values indicate a heat flux away from surface); B(σ) is the inverse of the Bowen ratio (see Eq. 7); and σ is a dimensionless parameter (see Eq. 8). I_s is the soil thermal inertia (J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$), calculated with the empirically derived equation from Huang and Wang (2016):

$$\begin{array}{}\text{(3)}& I_{\text{s}}=\sqrt{I_{\text{ds}}^{\mathrm{2}}+\mathit{\theta }{I}_{\text{w}}^{\mathrm{2}}},\end{array}$$

where I_ds is the dry soil thermal inertia (J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$), θ is the volumetric soil moisture (m³ m⁻³) and I_w is the thermal inertia of liquid water (1557 J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$). I₀ is the “apparent thermal inertia of air” (J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$) and was calculated as

$$\begin{array}{}\text{(4)}& I_{\mathrm{0}}={\mathit{\rho }}_{\text{a}}{c}_{\text{p}}\sqrt{C_{\mathrm{1}}kz}{\left(C_{\mathrm{2}}{\displaystyle \frac{kzg}{{\mathit{\rho }}_{\text{a}}{c}_{\text{p}}{T}_{\mathrm{0}}}}\right)}^{\mathrm{1}/\mathrm{6}},\end{array}$$

where ρ_a is the air density (1.22 kg m⁻³), c_p is specific heat of air under constant pressure (1004 J kg⁻¹ K⁻¹), C₁ and C₂ are two constants in the empirical functions representing the effects of stability on the mean profiles of wind speed and temperature within the surface layer (Wang and Bras, 2009), k is the von Kármán constant (0.4), z is the height above ground based on vegetation conditions (m), g is the gravitational acceleration (9.81 m s⁻²), and T₀ is a reference temperature (300 K).

For vegetated surfaces, the MEP model of transpiration considers the energy balance at the vegetation surface. As such, the term G corresponds to the canopy heat flux (not to be confused with the ground heat flux G in the MEP model for non-vegetated surfaces; Eqs. 1–2) and is considered negligible compared with the sensible and latent heat fluxes (Wang and Bras, 2011). Accordingly, the energy balance is simplified, and transpiration (E_t) is estimated by solving Eqs. (5) and (6) under the constraint of energy conservation at the surface ($R_{\text{n}}=E_{\text{t}}+H$):

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{E}_{\text{t}}={\displaystyle \frac{R_{\text{n}}}{\mathrm{1}+B^{-\mathrm{1}}\left(\mathit{\sigma }\right)}},\text{(6)}& {\displaystyle }& {\displaystyle }H={\displaystyle \frac{R_{\text{n}}}{\mathrm{1}+B\left(\mathit{\sigma }\right)}}.\end{array}$$

The same definition of B(σ) and σ applies for non-vegetated (Eqs. 1–2) and vegetated (Eqs. 5–6) surfaces, and the Bowen reciprocal ratio is calculated as

$$\begin{array}{}\text{(7)}& B\left(\mathit{\sigma }\right)=\mathrm{6}\left(\sqrt{\mathrm{1}+{\displaystyle \frac{\mathrm{11}}{\mathrm{36}}}\mathit{\sigma }}-\mathrm{1}\right),\end{array}$$

where σ corresponds to a dimensionless parameter that characterizes the phase change at the evaporating or transpiring surface:

$$\begin{array}{}\text{(8)}& \mathit{\sigma }={\displaystyle \frac{\sqrt{\mathit{\alpha }}{\mathit{\lambda }}^{\mathrm{2}}}{c_{\text{p}}{R}_{\text{v}}}}{\displaystyle \frac{q_{\text{s}}}{T_{\text{s}}^{\mathrm{2}}}},\end{array}$$

where α is the ratio of the eddy diffusivities for water vapour and heat (assumed to be unity; Wang et al., 2014), λ is the latent heat of vaporization of liquid water (J kg⁻¹), C_p is the specific heat of air under constant pressure (J kg⁻¹ K⁻¹), R_v is the gas constant of water vapour (461 J kg⁻¹ K⁻¹), q_s is the surface specific humidity (kg kg⁻¹), and T_s is the surface temperature (K).

In the MEP model of evaporation (MEP-E_s; Eqs. 1–2), the surface specific humidity (q_s) corresponds to the specific humidity at the soil surface (q_ss), and the surface temperature (T_s) corresponds to skin temperature at the soil surface (T_ss). Equation (8) can be modified accordingly:

$$\begin{array}{}\text{(9)}& \mathit{\sigma }={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{c_{\text{p}}{R}_{\text{v}}}}{\displaystyle \frac{q_{\text{ss}}}{T_{\text{ss}}^{\mathrm{2}}}}.\end{array}$$

The specific humidity at the soil surface can be computed as (Huang and Wang, 2016)

$$\begin{array}{}\text{(10)}& q_{\text{ss}}={\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\text{s}}}}\right)}^{\mathit{\beta }}{q}_{\text{sat}},\end{array}$$

where q_sat is the specific humidity at saturation (kg kg⁻¹) at T_ss, which is calculated using the Clausius–Clapeyron equation, θ is the soil water content (m³ m⁻³), θ_s is the soil porosity (m³ m⁻³), and β is an empirical parameter that was set to β=2 based on Huang and Wang (2016).

In the MEP model of transpiration (MEP-E_t; Eqs. 5–6), the surface specific humidity (q_s) corresponds to the specific humidity at the leaf surface (q_ls), and the surface temperature (T_s) corresponds to the leaf temperature (T_ls). A water stress factor (η_s) that varies between 0 and 1 is added to the original formulation of Eq. (8) presented by Wang and Bras (2011) to describe the reduction in plant transpiration under soil water stress (Hajji et al., 2018). For MEP-E_t, Eq. (8) thus translates into

$$\begin{array}{}\text{(11)}& \mathit{\sigma }={\mathit{\eta }}_{\text{s}}{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{c_{\text{p}}{R}_{\text{v}}}}{\displaystyle \frac{q_{\text{ls}}}{T_{\text{ls}}^{\mathrm{2}}}}.\end{array}$$

Various parameterizations of η_s exist in the literature based on soil water potential (Verhoef and Egea, 2014; Ferguson et al., 2016), volumetric soil water content (Feddes et al., 1978; Porporato et al., 2001) or leaf water potential (Tuzet et al., 2003). We used the Wang and Leuning (1998) parameterization of η_s based on soil water content that Hajji et al. (2018) successfully used to apply the MEP model of transpiration under water-limiting conditions:

$$\begin{array}{}\text{(12)}& {\mathit{\eta }}_{\text{s}}=min\left[\mathrm{1},{\displaystyle \frac{\mathrm{10}\left({\mathit{\theta }}_{\text{root}}-{\mathit{\theta }}_{\text{wp}}\right)}{\mathrm{3}\left({\mathit{\theta }}_{\text{fc}}-{\mathit{\theta }}_{\text{wp}}\right)}}\right],\end{array}$$

where θ_root is the weighted average soil water content over the root zone, with weights set according to the vertical root distribution (m³ m⁻³), θ_wp is the soil water content at the wilting point (m³ m⁻³) and θ_fc is the soil water content at field capacity (m³ m⁻³).

The MEP models have been formulated for the two limiting cases of bare soil evaporation and transpiration from a fully vegetated surface. In order to continuously simulate terrestrial evaporation, Hajji et al. (2018) proposed a method to combine the MEP models of evaporation and transpiration using a vegetation index (f_veg) that corresponds to the fraction of the soil covered by vegetation. Assuming that evaporation from intercepted rainfall is negligible, as is the case in this study, the MEP model of terrestrial evaporation (MEP-E) is then defined as

$$\begin{array}{}\text{(13)}& E=\left(\mathrm{1}-f_{\text{veg}}\right)E_{\text{s}}+f_{\text{veg}}{E}_{\text{t}}.\end{array}$$

The vegetation index (f_veg) can be derived from the normalized difference vegetation index (NDVI; Gutman and Ignatov, 1998):

$$\begin{array}{}\text{(14)}& f_{\text{veg}}={\displaystyle \frac{{\text{NDVI}}_t-{\text{NDVI}}_{\min }}{{\text{NDVI}}_{\max }-{\text{NDVI}}_{\min }}},\end{array}$$

where NDVI_t is the NDVI on day t, NDVI_min is the NDVI signal from bare soil and NDVI_max is the NDVI signal from a full vegetation cover.

3.2 Hydrological model

HydroGeoSphere (HGS) is an integrated surface and subsurface hydrologic model that has been used to simulate soil water content at various AmeriFlux sites (Maheu et al., 2018) as well as in a forested headwater catchment (Koch et al., 2016). Maxwell et al. (2014) performed a formal verification of seven integrated surface and subsurface models, and HGS showed good agreement with other models when simulating soil moisture in an idealized test case. HGS is a control-volume finite-element model that simultaneously solves the 2-D diffusion-wave approximation of the Saint-Venant equations and the 3-D form of the Richards equation (Aquanty, 2013). The surface and subsurface are coupled via a first-order exchange coefficient. The commonly used van Genuchten (1980) model is incorporated in HGS to describe the water retention curve. HGS has an adaptive time-stepping procedure where the time step length is controlled by the maximum change allowed in a state variable during any time step. In the present study, we allowed for a maximum change of 0.05 in soil saturation. In the current version of HGS, transpiration and evaporation are simulated as a function of potential evaporation (Kristensen and Jensen, 1975):

$$\begin{array}{}\text{(15)}& E_{\text{t}}\left(z\right)=f_{\mathrm{1}}\left(\text{LAI}\right)f_{\mathrm{2}}\left(\mathit{\theta }\right)\left(E_{\text{p}}-E_{\text{c}}\right)r\left(z\right),\end{array}$$

where E_t(z) is the transpiration rate (m s⁻¹) at depth z, E_p is the potential evaporation rate (m s⁻¹) and E_c is the wet canopy evaporation rate (m s⁻¹). f₁(LAI) is a function of the leaf area index (LAI), representing changes in vegetation over time (dimensionless):

$$\begin{array}{}\text{(16)}& f_{\mathrm{1}}\left(\text{LAI}\right)=max\left\{\mathrm{0},min\left[\mathrm{1},C_{\mathrm{2}}+C_{\mathrm{1}}\text{LAI}\right]\right\},\end{array}$$

where C₁ and C₂ are dimensionless fitting parameters. f₂(θ) is a function describing the water stress on vegetation (dimensionless), which is equivalent to η_s (Eq. 12). r(z) is the root distribution function, which follows a cubic decay distribution between the surface and the maximum root depth. Soil evaporation is described as

$$\begin{array}{}\text{(17)}& E_{\text{s}}\left(z\right)={\mathit{\alpha }}^{*}\left[\mathrm{1}-f_{\mathrm{1}}\left(\text{LAI}\right)\right]\left(E_{\text{p}}-E_{\text{c}}\right)e\left(z\right),\end{array}$$

where E_s(z) is the soil evaporation rate (m s⁻¹) at depth z, and e(z) is the evaporation depth function which follows a cubic decay distribution between the surface and the maximum evaporation depth. α^* is a soil wetness factor (dimensionless):

$$\begin{array}{}\text{(18)}& {\mathit{\alpha }}^{*}=\left\{\begin{array}{ll}\frac{\mathit{\theta }-{\mathit{\theta }}_{e\mathrm{2}}}{{\mathit{\theta }}_{e\mathrm{1}}-{\mathit{\theta }}_{e\mathrm{2}}}& \mathrm{for}{\mathit{\theta }}_{e\mathrm{2}}\le \mathit{\theta }\le {\mathit{\theta }}_{e\mathrm{1}},\\ \mathrm{1}& \mathrm{for}\mathit{\theta }>{\mathit{\theta }}_{e\mathrm{1}},\\ \mathrm{0}& \mathrm{for}\mathit{\theta }<{\mathit{\theta }}_{e\mathrm{2}},\end{array}\right.\end{array}$$

where θ_e1 is the soil water content above which full evaporation occurs and θ_e2 is the soil water content below which evaporation is zero.

Figure 2 Coupling of the HydroGeoSphere (HGS) hydrological model and the maximum entropy production (MEP) model of land surface fluxes. HGS supplies the MEP model with soil water saturations (S_w; in blue), which are used by the MEP model to compute the transpiration (E_t) and evaporation (E_s) rates. These rates are then passed on to HGS to compute the sink terms with the surface (Γ_ex; in red), and HGS removes water from the soil reservoir based on the root and evaporation depth profiles.

Table 2 Soil and vegetation properties at each study site.

^* US-Wkg: Nearing et al. (2005). US-Ton: BADM from the AmeriFlux database. US-WBW: Miller et al. (2007).

3.3 Coupling the HGS and MEP-E models

In the HGS model, the following modified form of Richards equation is used to describe the temporal evolution of soil moisture:

$$\begin{array}{}\text{(19)}& -\mathrm{\nabla }\cdot \mathit{q}+\sum {\mathrm{\Gamma }}_{\mathrm{ex}}\pm Q={\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\theta }}_{\text{s}}{S}_{\text{w}}\right),\end{array}$$

where q is the specific volumetric (Darcy) flux (m³ s⁻¹), Γ_ex is the source or sink term for exchange fluxes with other domains (i.e. surface; m³ s⁻¹), Q is the source or sink term for the porous medium such as pumping or injection (m³ s⁻¹), and S_w is the subsurface water saturation (m³ m⁻³). The MEP-E_s model requires soil moisture information to compute three parameters: the soil thermal inertia (I_s; Eq. 3) and the surface specific humidity of the soil (q_ss; Eq. 9) and leaf (q_ls; Eq. 11) surfaces. At the beginning of a given time step, HGS supplies the soil moisture information from the previous time step to the MEP-E model, which then computes the evaporation and transpiration rates (Fig. 2). The transpiration and evaporation rates computed with the MEP-E model are then transferred to HGS and taken as sinks for the porous medium (Γ_ex in Eq. 19). HGS removes water from the soil reservoir based on the root and evaporation depth profiles, which both follow a cubic decay distribution between the surface and the maximum root or evaporation depth. At the end of the given time step, HGS calculates the saturation throughout the soil column. Given that the MEP-E_s model was not very sensitive to changes in soil thermal inertia, this parameter was set as a constant equal to the dry soil thermal inertia (as detailed in the next section) and was not involved in the coupling procedure.

Table 3 Implementation of the MEP model of total terrestrial evaporation in this study.

4 Model implementation

4.1 HGS-MEP

The coupled HGS-MEP model was not calibrated in the present study. The soil and vegetation properties were therefore not obtained from calibration but were rather defined from pedotransfer functions, observations at AmeriFlux sites or values taken from the literature. We derived soil hydraulic characteristics (Table 2) using soil texture information as an input to the Rosetta model implemented in Hydrus (Šimůnek et al., 2013) to derive the soil hydraulic characteristics. Soil texture information was either taken from the AmeriFlux database (US-Ton) or from the literature (US-Wkg in Nearing et al., 2005, and US-WBW in Miller et al., 2007). At US-Wkg and US-Ton, the observed minimum soil water content was smaller than residual saturation computed by the Rosetta model, and as such, we set the residual saturation in HGS to the observed minimum soil water content. At each site, we defined soil porosity as the observed maximum soil water content closest to the surface, and we used the water retention curve to compute the volumetric water content corresponding to the wilting point (−1.5 MPa) and field capacity (−0.033 MPa). We defined the maximum root depth using site-specific information at US-Ton (Ichii et al., 2009) and US-WBW (Wilson et al., 2001). Given that similar information was unavailable at US-Wkg, we set the maximum root depth based on the root depth reported for Lehman lovegrass in the literature (Gibbens and Lenz, 2001), the main grass species at the site. Given that only little information is available on soil and root properties, parameter allocation for these properties is often an important source of uncertainty in hydrological models. However, previous modelling work at the study sites has shown soil moisture modelling to be relatively robust to variations in saturated hydraulic conductivity and vertical root distribution (Maheu et al., 2018).

Table 3 provides a summary of the implementation of the MEP model of terrestrial evaporation in HGS. In the MEP model of evaporation, skin temperature measurements would ideally be used to set the soil surface temperature (T_ss), but since those measurements were not available, we instead used soil temperature measurements nearest to the surface, which is at a depth of 5 cm at US-Wkg, 2 cm at US-Ton and 2 cm at US-WBW. The height above ground (z in Eq. 4) was set to the flux tower height, which was 6.4 m at US-Wkg, 23 m at US-Ton and 36.9 m at US-WBW. At US-Wkg, the dry soil thermal inertia (I_ds in Eq. 3) was computed according to Wang et al. (2010) as the regression coefficient between diurnal variations in the ground heat flux and surface temperature (830 J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$). Because these data were unavailable at US-Ton and US-WBW and given that the model showed little sensitivity to the soil thermal inertia, the dry soil thermal inertia was set equal to 800 J m⁻² K⁻¹ s${}^{-\mathrm{1}/\mathrm{2}}$ for these two sites. In the MEP model of transpiration, the leaf surface temperature was assumed to be equal to the air temperature measured above the canopy. Since no measurements of the leaf surface specific humidity were available, we used the specific humidity of the air as a proxy and calculated it from air temperature and relative humidity measurements using the Clausius–Clapeyron equation. To combine the MEP models of evaporation and transpiration (Eq. 13), we computed the vegetation index (Eq. 14) using the AVHRR 7 d composite NDVI (https://lta.cr.usgs.gov/NDVI, last access: 4 November 2016). NDVI time series are notoriously noisy because of varying atmospheric conditions and sensor viewing angles (Hird and McDermid, 2009), and we therefore smoothed the time series by applying a 60 d moving average. Based on Montandon and Small (2008), the NDVI signal for bare soil (NDVI₀) was set to 0.2, and the NDVI signal for full vegetation cover varied according to the land cover, with NDVI_∞=0.61 for grassland, NDVI_∞=0.69 for woody savanna and NDVI_∞=0.85 for deciduous broadleaf forest.

4.1.1 HGS with Penman–Monteith (HGS-PM)

We compared HGS-MEP to simulations of terrestrial evaporation and soil water content performed with the stand-alone HGS model, in which the Kristensen and Jensen (1975) model is implemented, as described in Sect. 3.2. We used the same soil and root properties defined for the HGS-MEP simulations (Table 2), and the only difference was the terrestrial evaporation model. In HGS, transpiration and evaporation are set to occur simultaneously. However, this approach led to a large overestimation of terrestrial evaporation. In the present study, we modified the implementation of the model, and evaporation took place only when f₁(LAI)=0. As for the HGS-MEP model implementation, interception was not considered and E_c was set to zero. The maximum and minimum evaporation limiting saturation (Eq. 18) were set to 0.1 (θ_e1=0.1θ_s) and 0.5 (θ_e2=0.5θ_s; Verbist et al., 2012). Potential evaporation (E_p), which is terrestrial evaporation from a saturated land surface, was computed with the Penman–Monteith equation, a physically based model where, similar to the MEP model, the predicted terrestrial evaporation is constrained by available energy at the surface:

$$\begin{array}{}\text{(20)}& E_{\text{p}}={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}{\displaystyle \frac{\mathrm{\Delta }\left(R_{\text{n}}-G\right)+\frac{{\mathit{\rho }}_{\text{a}}{c}_{\text{p}}\left(e_{\text{s}}-e_{\text{a}}\right)}{r_{\text{a}}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\text{s}}}{r_{\text{a}}}\right)}},\end{array}$$

where Δ is the slope of the saturation vapour pressure curve (kPa K⁻¹), e_s is the saturation vapour pressure (kPa), e_a is the actual vapour pressure (kPa), γ is the psychrometric constant (kPa K⁻¹), r_s is the surface resistance (m s⁻¹) and r_a is the aerodynamic resistance (m s⁻¹). When f₁(LAI)>0, the surface resistance (r_s) was set according to the vegetation lookup table in the Noah land surface model (US-Wkg = 40 s m⁻¹, US-Ton = 70 s m⁻¹ and US-WBW = 100 s m⁻¹ in Kumar et al., 2011), and the canopy aerodynamic resistance (r_ac) was computed as (Thom, 1975)

$$\begin{array}{}\text{(21)}& r_{\text{ac}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}u}}\left[\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}m}}}\right)\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}v}}}\right)\right],\end{array}$$

where u is the wind speed (m s⁻¹), z is the wind speed measurement height (m), d₀ is the zero-plane displacement height (m), z_0m is the roughness height for momentum transfer (m) and z_0v is the roughness height for water vapour transfer (m). Equation (21) was derived for neutral atmospheric conditions but has also been successfully used to model terrestrial evaporation over a wide range of conditions (Ershadi et al., 2014). Roughness heights were estimated as a fraction of the vegetation height, h (m; Brutsaert, 1982):

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{d}_{\mathrm{0}}=\mathrm{0.66}h,\text{(23)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{0}m}=\mathrm{0.1}h,\text{(24)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{0}v}=\mathrm{0.1}{z}_{\mathrm{0}m}.\end{array}$$

When f₁(LAI)=0, the surface resistance (r_s) was set to 999 s m⁻¹, the value associated with a barren or sparsely vegetated land cover in the Noah lookup table (Kumar et al., 2011), and the substrate aerodynamic resistance (r_as) was computed as (Shuttleworth and Wallace, 1985)

$$\begin{array}{}\text{(25)}& r_{\text{as}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}u}}\left[\ln \left({\displaystyle \frac{z}{{z^{\prime }}_{\mathrm{0}}}}\right)\ln \left({\displaystyle \frac{d_{\mathrm{0}}+z_{\mathrm{0}m}}{{z^{\prime }}_{\mathrm{0}}}}\right)\right],\end{array}$$

where ${z^{\prime }}_{\mathrm{0}}$ is the roughness length of the soil (m), which was set to 0.01 m (Shuttleworth and Wallace, 1985).

The f₁(LAI) function used to describe temporal changes in vegetation (Eq. 16) was computed by rescaling, between 0 and 1, the vegetation index f_veg, an input to the HGS-MEP model.

4.2 Model setup

As a proof of concept, we set up the coupled HGS-MEP and HGS-PM models to perform one-dimensional soil column simulations to evaluate the capability of both models to simulate water fluxes (terrestrial evaporation) and storage (soil moisture). We represented the soil column with a fine (1 cm) vertical resolution and set the soil column depth to either 1 or 1.5 m in order to capture the entire root zone at each site. We assigned uniform soil properties throughout the soil column, since there were no available data to describe the vertical distribution of soil material with depth. Simulations spanned 5 years at US-Wkg and US-Ton and 2.5 years at US-WBW given that data were not available for a longer period. We used soil water content measurement available at different depths to set initial subsurface conditions, and we used linear interpolation to assign initial conditions at depths without measurements. As for initial surface water conditions, we assumed an initial surface water depth of zero given that the soil was not fully saturated. For boundary conditions, we supplied the model with gap-filled (REddyProc; Reichstein et al., 2005) measurements of precipitation, net radiation, air temperature, relative humidity and soil temperature at a 30 min time step. At the surface, we applied a critical depth boundary condition that allows water to leave the model domain via overland flow. At the bottom of the soil column, we applied a free drainage boundary condition.

4.3 Model performance

We evaluated the performance of models using a series of metrics comparing observed and simulated values of terrestrial evaporation and soil water content at the three sites. First, we computed the root-mean-square error (RMSE) to assess the mean difference between observed and simulated values. Second, we computed the Nash–Sutcliffe efficiency (NSE), where a value of 1 indicates a perfect agreement between the model and observations and a negative value indicates that the average value of observations offers a better predictor than the model. The RMSE and NSE are not independent metrics of performance given that the NSE is a standardized measure of the mean square error. Still, we chose to report the two commonly used metrics, as they provide an assessment of performance in absolute (RMSE) and relative (NSE) terms. Third, we computed the normalized benchmark efficiency (BE), which is analogue to the NSE but compares the model output to a simple benchmark model, in this case being the interannual mean value for every calendar day (Schaefli and Gupta, 2007). Fourth, we computed the coefficient of determination (R²) that describes the proportion of the total variance in observations explained by the model. Finally, we computed the percent bias (PBIAS) to assess the average tendency of simulated values to be larger (positive bias) or smaller (negative bias) than observations. Equations of performance metrics are listed in Table S1 in the Supplement.

Figure 3 (a) Precipitation as well as observed and modelled (b) daily mean soil water content at a depth of 15 cm, and (c) 10 d moving average of terrestrial evaporation at US-Wkg (climate: semi-arid; vegetation: grassland).

We calculated performance metrics at half-hourly and daily timescales for terrestrial evaporation and at a daily timescale for soil water content. When assessing the ability of the models to simulate terrestrial evaporation, we first calculated performance metrics on the entire time series and, second, assessed how well the models performed under water-limited conditions. To do so, we computed the monthly aridity index, which is the ratio between precipitation and PET, where PET was calculated with the Penman–Monteith equation (Eq. 20). Using monthly values of the aridity index, we then calculated performance metrics for dry periods (P ∕ PET < 0.4 at US-Wkg and P ∕ PET < 1 at US-Ton and US-WBW) and wet periods (P ∕ PET ≥ 0.4 at US-Wkg and P ∕ PET ≥ 1 at US-Ton and US-WBW). As US-Wkg is located in a semi-arid climate, we used a different threshold (0.4) to define water-limited conditions, as the monthly aridity index remained below 1 for the entire study period.

5 Results

5.1 Model performance under a semi-arid climate (US-Wkg)

During the study period (2010–2014), the mean annual precipitation at US-Wkg varied between 264 and 415 mm, with 2013 and 2014 being the driest and wettest years, respectively (Fig. 3a). Between January and June, pre-monsoon precipitation was relatively low (≤ 35 mm), with the exception of 2010, which saw 110 mm of rain in the first 6 months of the year. Precipitation was concentrated during the monsoon periods, and, with the exception of 2010, nearly 80 % occurred between July and September. Figure 3b shows that the HGS-MEP model provided a realistic simulation of soil moisture at 15 cm depth (RMSE = 0.04 m³ m⁻³; NSE = 0.30; Table 4) and was able to capture the sharp rise in soil moisture at the start of the monsoon period in July – note that soil moisture observations were missing between 5 May 2010 and 31 December 2011. Overall, soil moisture at 15 cm depth was generally overestimated (PBIAS = 19 %; Table 4) with HGS-MEP, particularly for the lower values outside the monsoon period. For example, the observed annual minimum in soil water content was 0.06 m³ m⁻³ between 2012 and 2014, while modelled soil moisture did not fall below 0.12 m³ m⁻³. A similar overestimation was also observed for the upper (z=5 cm) and lower (z=30 cm) soil layers (Fig. S1 in the Supplement).

Table 4 Performance of the HGS-MEP and HGS with Penman–Monteith (HGS-PM) models when simulating (a) daily mean soil water content (SWC) at a 15 or 20 cm depth and (b) daily mean total terrestrial evaporation (E) as represented by the root-mean-square error (RMSE), Nash–Sutcliffe efficiency (NSE), benchmark efficiency (BE), coefficient of determination (R²) and percentage bias (PBIAS).

Table 5 Performance of the HGS-MEP and HGS with Penman–Monteith (HGS-PM) models when simulating daily mean total terrestrial evaporation (E) during dry and wet periods, as represented by the root-mean-square error (RMSE), Nash–Sutcliffe efficiency (NSE), coefficient of determination (R²) and percentage bias (PBIAS). n represents the number of days.

Figure 3c shows that the HGS-MEP model also performed well in simulating the daily mean terrestrial evaporation (RMSE = 0.31 mm d⁻¹; NSE = 0.88; Table 4). The HGS-MEP model reproduced the seasonal pattern in terrestrial evaporation and captured the sharp increase associated with the increased soil water availability from July to September. Peak terrestrial evaporation was slightly underestimated by HGS-MEP (PBIAS = −10 %; Table 4). For example, the observed annual maximum ranged between 4.8 and 5.8 mm d⁻¹, while the modelled annual maximum did not exceed 4.0 mm d⁻¹. The HGS-MEP model was also able to simulate the increase in terrestrial evaporation following precipitation events outside the monsoon months. However, the pre-monsoon months were particularly wet in 2010, and terrestrial evaporation simulated by HGS-MEP was underestimated, with a modelled average of 0.5 mm d⁻¹ between January and May compared with an observed average of 0.9 mm d⁻¹. At US-Wkg, terrestrial evaporation was dominated by evaporation, particularly during the dry season, during which it represented 69 % of the terrestrial evaporation simulated by HGS-MEP (Fig. 4a). During the monsoon period when vegetation activity is concentrated (July–October), evaporation decreased in importance and, on average, accounted for 60 % of the total terrestrial evaporation simulated by HGS-MEP, while transpiration represented 40 % of total terrestrial evaporation. These proportions varied from year to year. For example, 2014 was the wettest year during the study period, and modelled transpiration represented 46 % of total terrestrial evaporation during the monsoon period.

Figure 4 10 d moving average of observed terrestrial evaporation and modelled terrestrial evaporation partitioned as transpiration and evaporation at (a) US-Wkg (climate: semi-arid; vegetation: grassland), (b) US-Ton (climate: Mediterranean; vegetation: woody savanna) and (c) US-WBW (climate: temperate; vegetation: deciduous broadleaf forest).

Figure 5 Observed and modelled hourly average terrestrial evaporation at (a) US-Wkg (climate: semi-arid; vegetation: grassland), (b) US-Ton (climate: Mediterranean; vegetation: woody savanna) and (c) US-WBW (climate: temperate; vegetation: deciduous broadleaf forest).

Table 6 Performance of the HGS-MEP and HGS with Penman–Monteith (HGS-PM) models when simulating half-hourly mean total terrestrial evaporation, as represented by the root-mean-square error (RMSE), Nash–Sutcliffe efficiency (NSE), coefficient of determination (R²) and percentage bias (PBIAS).

In contrast, modelled transpiration represented 28 % of total terrestrial evaporation during the monsoon period in 2012, although terrestrial evaporation was overall underestimated by the HGS-MEP model that year. Overall, the HGS-MEP model outperformed the HGS-PM model for the simulation of daily terrestrial evaporation and soil moisture (Fig. 3b). Indeed, performance metrics at the daily timescale show a large decline in the performance of HGS-PM during wet periods (NSE = 0.04) compared with dry periods (NSE = 0.64; Table 5). In contrast, the performance of HGS-MEP was comparable for dry (NSE = 0.86) and wet periods (NSE = 0.85). Moreover, in between monsoon periods, terrestrial evaporation was underestimated by HGS-PM (PBIAS = −14 %; Table 4), and the model was unable to catch terrestrial evaporation pulses following rain events (e.g. 2013; Fig. 3c). As a result, soil moisture was generally overestimated (PBIAS = 27 %) throughout the simulation period. At the diurnal scale, the HGS-MEP model reproduced subdaily variations in terrestrial evaporation at US-Wkg well, with an increase in the morning, peak value around 12:00 LT, a decrease in the afternoon, and values close to zero during the night (Fig. 5a). In the HGS-MEP simulation, the daily maximum was generally overestimated, with an average of 0.12 mm h⁻¹ at 12:00 LT compared with an observed average of 0.10 mm h⁻¹. On the other hand, morning and afternoon values were slightly underestimated by about 0.01 mm h⁻¹ (Fig. 5a). While the HGS-PM did not overestimate the daily maximum (Fig. 5a), the HGS-MEP model overall outperformed HGS-PM for the simulation of half-hourly terrestrial evaporation, as shown by the four performance metrics (RMSE, NSE, R² and PBIAS; Table 6).

5.2 Model performance under a Mediterranean climate (US-Ton)

During the study period (2004–2008), mean annual precipitation at US-Ton ranged between 371 and 782 mm, with most precipitation concentrated during winter months, from October to May (Fig. 6a). The first two winters (2004–2005 and 2005–2006) of the study period were particularly wet (precipitation is 717 and 882 mm), while the following years experienced near-normal precipitation (385 and 392 mm). As shown in Fig. 6b, the HGS-MEP model simulated soil moisture exceptionally well at a 20 cm depth (RMSE = 0.03 m³ m⁻³; NSE = 0.92; Table 4) and was able to reproduce the decrease in soil moisture as precipitation stops in summer (Fig. 6b). During the wet winter months, the HGS-MEP model slightly underestimated soil moisture (PBIAS = −5 %) but still captured the increase in soil moisture during this period. Near the surface (z= 0 cm; Fig. S2a), soil moisture was generally underestimated by HGS-MEP during the dry summer months: observed soil water content slowly decreased until it generally reached a minimum of 0.04–0.05 m³ m⁻³ at the end of the season, while modelled soil moisture dropped rapidly to reach a minimum value (0.03 m³ m⁻³) close to the residual water content (0.02 m³ m⁻³; Table 2). For the deeper soil layers (z= 50 cm), the HGS-MEP model generally overestimated soil moisture during dry summer months, with a modelled minimum of 0.17 m³ m⁻³ compared to an observed minimum of 0.14 m³ m⁻³ (Fig. S2c).

When simulating terrestrial evaporation, the HGS-MEP model performed well (RMSE = 0.43 mm d⁻¹; NSE = 0.73), although terrestrial evaporation was underestimated (PBIAS = −14 %; Table 4). The HGS-MEP model performed particularly well for winter months and was able to capture the increase in terrestrial evaporation from about 0.4 mm d⁻¹ (modelled average terrestrial evaporation in September) to 3.3 mm d⁻¹ (modelled average annual maximum terrestrial evaporation) as water became more available. However, following those winter months, as precipitation stopped, terrestrial evaporation simulated by HGS-MEP was underestimated during the first half of the dry period (June–July). For example, the underestimation in terrestrial evaporation during the summer of 2005 and 2006 (i.e. from the last day of precipitation in May or June to the end of September) amounted respectively to a cumulative difference of 8.1 and 6.8 mm between observed and modelled terrestrial evaporation. During the dry summer months, low soil moisture availability near the surface (Fig. S2a) limited evaporation simulated by HGS-MEP, and modelled transpiration accounted on average for 94 % of total terrestrial evaporation (Fig. 4b). Evaporation increased during the monsoon period, and on average, terrestrial evaporation simulated by HGS-MEP was made up of 17 % of evaporation and 83 % of transpiration. Overall, the HGS-PM model performed well at US-Ton at a daily timescale, although the HGS-MEP performed slightly better. As opposed to HGS-MEP, the HGS-PM had difficulties capturing the onset of the wet winter period, and the increase in soil moisture occurred ahead of time (Fig. 6b). Indeed, we observed a decline in the performance of HGS-PM during wet periods (NSE = 0.20) compared with dry periods (NSE = 0.59; Table 5). We observed a similar pattern with HGS-MEP, but it was not as marked, with a NSE of 0.37 during wet periods compared with 0.76 during dry periods (Table 5). With HGS-MEP and HGS-PM, terrestrial evaporation was underestimated, both during wet and dry periods (−27 % ≤ PBIAS ≤ −10 %; Table 5). This underestimation was particularly important in June and July, when precipitation stopped (Fig. 6c). At the diurnal scale, both models performed similarly well (NSE = 0.53 for HGS-MEP and NSE = 0.55 for HGS-PM), although HGS-PM led to a more important bias (PBIAS = −18 %) than HGS-MEP (PBIAS = −3 %; Table 6). Throughout the day, both models underestimated terrestrial evaporation in the morning and afternoon, while the daily maximum was generally overestimated by HGS-MEP and slightly underestimated by HGS-PM (Fig. 5b).

Figure 6 (a) Precipitation as well as observed and modelled (b) daily mean soil water content at a depth of 20 cm, and (c) 10 d moving average of terrestrial evaporation at US-Ton (climate: Mediterranean; vegetation: woody savanna).

5.3 Model performance under a temperate climate (US-WBW)

In the simulation period, 2004 was a relatively wet year (annual precipitation = 1600 mm), while 2005 was relatively dry (annual precipitation = 995 mm; Fig. 7a). In Fig. 7b, the HGS-MEP model provided a realistic simulation of soil moisture (RMSE = 0.05 m³ m⁻³; NSE = 0.61; Table 4) and captured the decrease in soil moisture associated with increased vegetation activity during the summer. However, soil moisture at a 20 cm depth was generally overestimated by HGS-MEP during the summer, and modelled soil water content did not fall below 0.1 m³ m⁻³, while the observed soil water content reached an annual minimum of 0.08 and 0.05 m³ m⁻³ in 2004 and 2005. In contrast, soil moisture during the winter was generally underestimated by HGS-MEP, with a modelled average soil water content of 0.21 m³ m⁻³ between January and April in comparison with an observed average of 0.26 m³ m⁻³. Overall, a similar pattern (overestimation of soil moisture in the summer and underestimation in the winter) was observed in soil moisture simulations near the surface (z= 5 cm; Fig. S3a). As for deeper soil layers (z= 60 cm), soil moisture was systematically underestimated, although this could be due to changes in soil properties given that there is an upper shift in observed soil moisture values compared with upper soil layers (Fig. S3c).

In Fig. 7c, the HGS-MEP model performed well and reproduced the seasonal pattern in terrestrial evaporation (RMSE = 0.71 mm d⁻¹; NSE = 0.65; Table 4). Capturing the onset of vegetation activity was, however, challenging for the model, and early in the summer, terrestrial evaporation was generally overestimated (PBIAS = 11 %). On two occasions, terrestrial evaporation was considerably underestimated by the HGS-MEP model: in June 2005, observed terrestrial evaporation reached a peak of 3.9 mm d⁻¹, while modelled terrestrial evaporation was about 2.5 mm d⁻¹, and in September 2005, observed terrestrial evaporation reached 2.8 mm d⁻¹, while modelled terrestrial evaporation went below 1 mm d⁻¹ (values refer to the 10 d moving average presented in Fig. 4c). This underestimation in terrestrial evaporation is linked to the water stress factor. In both cases, the water stress factor fell below 0.5 for a few days, which largely reduced modelled transpiration. Evaporation simulated by HGS-MEP was relatively constant throughout the year, with an average of 0.4 mm d⁻¹, although in the fall 2005, following a particularly dry summer, modelled evaporation dropped to less than 0.1 mm d⁻¹. Transpiration simulated by HGS-MEP accounted on average for 57 % of total terrestrial evaporation between October and April, but this proportion increased considerably during the summer. Indeed, between May and September, modelled transpiration represented on average 87 % of total terrestrial evaporation. Overall, the HGS-MEP model outperformed HGS-PM, although, as shown by negative BE values, both models had less explanatory power than a simple benchmark model that captures seasonality. The two models performed similarly during dry and wet periods, with only a slight decline in performance at HGS-PM during dry periods in the summer (NSE = 0.42) compared with wet periods in the winter (NSE = 0.66; Table 5). During summer months, the HGS-PM largely underestimated terrestrial evaporation (PBIAS = −23 %), particularly during the dry 2005 year (Fig. 7c). As a result, soil moisture was overestimated by HGS-PM much more so than with HGS-MEP (Fig. 7b). During winter months, soil moisture was generally underestimated, although a similar pattern was observed for HGS-MEP. At the diurnal scale, the HGS-PM (NSE = 0.62) performed better than HGS-MEP (NSE = 0.47) when simulating half-hourly terrestrial evaporation, although both models had an important bias, which was positive for HGS-MEP (PBIAS = 23 %) and negative for HGS-PM (PBIAS = −18 %; Table 6). This bias is particularly reflected in the simulation of the daily maximum, which was overestimated by HGS-MEP and largely underestimated by HGS-PM (Fig. 5c). Indeed, the daily maximum observed mid-day at US-WBW is 0.18 mm h⁻¹, while it reached 0.22 mm h⁻¹ with HGS-MEP and 0.13 mm h⁻¹ with HGS-PM.

Figure 7 (a) Precipitation as well as observed and modelled (b) daily mean soil water content at a depth of 20 cm, and (c) 10 d moving average of terrestrial evaporation at US-WBW (climate: temperate; vegetation: deciduous broadleaf forest).

6 Discussion

6.1 Performance of the HGS-MEP model

6.1.1 HGS-MEP vs. HGS-PM

At the daily timescale, the HGS-MEP model outperformed HGS-PM when simulating terrestrial evaporation at the three study sites, which translated into improved performance for soil moisture modelling as well (Table 4). Notably, we observed a weak performance of HGS-PM during wet periods at US-Wkg (P ∕ PET ≥ 0.4) and US-Ton (P ∕ PET ≥ 1; Table 5). At US-Wkg, this meant that HGS-PM failed to capture the annual maximum terrestrial evaporation, which has important implications for the water budget. In contrast, the weaker performance at US-Ton by HGS-PM, and to a lesser degree by HGS-MEP, meant that the two models struggle to describe the minimum rates of terrestrial evaporation. At US-WBW, both models had a comparable performance, with consistent results for HGS-MEP for wet (P ∕ PET ≥ 1) and dry periods, while HGS-PM performed better during wet periods than dry periods (Table 5).

At the semi-arid site US-Wkg (Fig. 3), a better handling of terrestrial evaporation partitioning can explain the superior performance of HGS-MEP compared with HGS-PM at both daily and half-hourly timescales. Indeed, the HGS-MEP model assesses transpiration and evaporation independently of one another. In the Kristensen and Jensen (1975) model in HGS, transpiration and evaporation are instead derived from a single PET value, although we could realistically expect potential transpiration and evaporation to differ from one another. In fact, Shuttleworth and Wallace (1985) have proposed a two-layer configuration of the Penman–Monteith model which allows the definition of different resistance values for the soil and canopy layers. However, HGS remains limited to the definition of a single PET value at the moment. We encountered a large overestimation of terrestrial evaporation when deriving transpiration and evaporation from a single PET value, and we set up the HGS model to only allow bare soil evaporation in the absence of active vegetation (f₁(LAI)=0). As a result, the HGS-PM model performed weakly at US-Wkg, where soil evaporation makes up more than 50 % of terrestrial evaporation according to observations (Moran et al., 2009). Soil moisture modelling by HGS-MEP also proved to be somewhat challenging at US-Wkg (NSE = 0.30, BE = −0.35; Table 4). However, in the absence of parameter calibration, we consider this performance to be more than acceptable given the challenges of modelling water fluxes under arid climates. Soil moisture was generally overestimated by HGS-MEP during dry periods (Fig. 3), and further work on the definition of the wilting point in arid climates could help improve the performance of the model. Indeed, wilting can occur at a lower threshold than −1.5 MPa (−4 to −2 MPa; Baldocchi et al., 2004) for vegetation having evolved under an arid climate.

At the Mediterranean site US-Ton, daily terrestrial evaporation was underestimated by HGS-MEP (PBIAS = −14 %; Table 4), particularly during the second half of the year, as reduced water supply led to a decline in terrestrial evaporation (Fig. 5). While the HGS-MEP model simulates soil moisture very well at a depth of 20 cm (Fig. 6), it tended to underestimate soil moisture close to the surface (Fig. S2a), where the largest proportion of roots is found according to the vertical root distribution defined by HGS (cubic decay distribution between the surface and the maximum root depth). Given that the water stress factor (η_s) is computed from the weighted average soil water content over the root zone (Eq. 12), this underestimation of soil moisture near the surface translated into an overestimation of the reduction in transpiration resulting from water stress. We also investigated if the issue of water stress overestimation could be due to a mis-definition of the maximum rooting depth, as trees under a Mediterranean climate have been found to access water from deep soil layers or groundwater (Miller et al., 2010). However, we also simulated terrestrial evaporation with the stand-alone MEP model using soil moisture observations, thus avoiding the overestimation of water stress near the surface, and instead found a large overestimation of terrestrial evaporation (Fig. S4). These results suggest that increasing the rooting depth to increase access to water resources would likely not improve the simulation of terrestrial evaporation. Instead, uncertainty relative to the definition of the vertical root distribution (as opposed to the maximum rooting depth) or, as previously discussed, with the definition of water stress points (wilting point and field capacity) may explain the challenge of simulating terrestrial evaporation under water-limited conditions at US-Ton.

At the half-hourly timescale, both HGS-MEP and HGS-PM showed lower performance than at a daily timescale when simulating terrestrial evaporation. For example, the NSE varied between 0.56 and 0.88 at a daily timescale (Table 4), while it varied between 0.47 and 0.65 at a half-hourly timescale (Table 6). At a half-hourly timescale, no model is distinctly superior to another: HGS-MEP performed better at US-Wkg, HGS-PM performed better at US-WBW and both models performed similarly at US-Ton (Table 6). Still, we observed a distinct pattern where peak terrestrial evaporation during the day was generally overestimated by HGS-MEP and underestimated by HGS-PM (Fig. 5). A known issue of energy imbalance, particularly at subdaily timescales, is associated with eddy covariance measurements (Leuning et al., 2012). This issue typically leads to an underestimation in observations of terrestrial evaporation and could in part explain the apparent issue of overestimation of peak values by HGS-MEP. On the other hand, the negative bias of HGS-PM at the half-hourly timescale (−18 % ≤ PBIAS ≤ −9 %; Table 6) may actually be even more important when taking the energy imbalance issue into account.

Overall, the predictive ability of the MEP model is particularly noteworthy given that we did not rely on calibration and instead used a priori estimation of parameters describing the soil and vegetation. The MEP model thus offers a promising alternative to model hydrologic fluxes without relying on calibration (Wagener, 2007). We chose the Penman–Monteith model as a benchmark against which to compare the MEP model, as it is a physically based model that allows for a detailed parameterization of vegetation. However, studies have shown that the Penman–Monteith model leads to an underestimation of terrestrial evaporation under a contemporary climate (Ershadi et al., 2014) and an overestimation under climate change (Milly and Dunne, 2016). Other models could have been considered, although Hajji et al. (2018) demonstrated the superior performance of the MEP model compared with other models such as the modified Priestley–Taylor Jet Propulsion Laboratory (PT-JPL) and the air-relative-humidity-based two-source model (ARTS).

6.1.2 HGS-MEP vs. HGS using observed terrestrial evaporation as a forcing

Maheu et al. (2018) assessed HGS' skills to model soil moisture at the same AmeriFlux sites used in the present study. In Maheu et al. (2018) and the present study, the same input values were used to define soil and vegetation properties. Both studies also considered the same periods, with the exception of the Mediterranean site (US-Ton), where simulations were performed for different periods (2004–2008 in the present study vs. 2008–2012 in Maheu et al., 2018) due to data availability. When modelling soil moisture with HGS, Maheu et al. (2018) used observed terrestrial evaporation as a forcing, thus reducing the uncertainty of simulating this flux and focusing the assessment of the model on subsurface processes and properties that control soil water content. These simulations thus offer a benchmark against which to compare the results of the present study and assess how much uncertainty is introduced by the simulation of terrestrial evaporation by the MEP model. At the three sites, the HGS simulations of soil moisture that used observed terrestrial evaporation as a forcing performed slightly better than simulations with HGS-MEP. Still, the performance was overall comparable between HGS with observed terrestrial evaporation (RMSE = 0.04 m³ m⁻³ and NSE between 0.39 and 0.86) and HGS-MEP (RMSE between 0.03 and 0.05 m³ m⁻³ and NSE between 0.30 and 0.92). At the semi-arid site (US-Wkg), the HGS-MEP simulation of soil moisture had a greater overestimation bias (19 %; Table 4) than the HGS simulation (10 %; Table 4 in Maheu et al., 2018), which could be due to the underestimation of peak (2011 and 2012) or pre-monsoon (2010, 2012 and 2013) terrestrial evaporation (Fig. 3). At the temperate site (US-WBW), simulations with HGS-MEP and HGS in Maheu et al. (2018) both showed the same pattern of underestimation of soil moisture in the winter and overestimation in the summer. This suggests that these biases are in good part associated with the definition of soil (hydraulic parameters) or vegetation (root distribution) properties rather than with the MEP model itself.

6.1.3 Partitioning of total terrestrial evaporation by HGS-MEP

In the present study, evaporation simulated by HGS-MEP represented on average 60 % of total terrestrial evaporation at the semi-arid site (US-Wkg) during the monsoon period (Fig. 4a). These results are in line with experimental results from Moran et al. (2009), who estimated that, following the Lehmann lovegrass invasion, evaporation accounted for 55 % of the total terrestrial evaporation during the growing season of a year with average precipitation. Results from the present study are also concordant with those of Scott and Biederman (2017), who found that, on average, evaporation represented 54 % of growing-season total terrestrial evaporation at US-Wkg between 2004 and 2015. Using overstory and understory flux tower measurements at the Mediterranean site US-Ton, Miller et al. (2010) found that transpiration from trees dominate during the dry summer months and that understory evaporation (i.e. bare soil evaporation and transpiration from grasses and forbs) is close to zero. Our results are consistent with these observations, and according to the HGS-MEP simulation, transpiration accounted on average for 94 % of total terrestrial evaporation (Fig. 4b). At the temperate site (US-WBW), soil evaporation, as measured by an understory flux tower, accounted for 16 % of total terrestrial evaporation on an annual basis and for generally less than 8 % of total terrestrial evaporation during the growing season (Wilson et al., 2001). In the present study, soil evaporation simulated by HGS-MEP amounted to 13 % of total terrestrial evaporation during the growing season between May and September (Fig. 4c), which agrees with these experimental estimates. Overall, the HGS-MEP model showed a good capability for the partitioning of total terrestrial evaporation under various climates (semi-arid, Mediterranean and temperate); this is a conclusion that was also reached for a humid, energy-limited environment (Wang et al., 2017).

6.2 Using the MEP model to integrate the energy budget in hydrological modelling: strengths and limitations

The MEP model of land surface fluxes offers an effective means to implement coupled water and energy budget modelling for hydrological applications. First, the MEP model of terrestrial evaporation requires six input variables: net radiation, soil surface temperature, leaf surface temperature and specific humidity, vegetation index and soil water content, the last of which can be supplied by a hydrological model. Thus, the MEP model eliminates the need for wind speed and surface roughness (input to the Penman model) or for vertical gradients of temperature and humidity (inputs to the aerodynamic method often implemented in land surface models). Second, the MEP model ensures, by design, the closure of the energy balance. As such, terrestrial evaporation simulated by the MEP model is constrained by available energy, which avoids issues of overestimation associated with the use of temperature-based PET models for hydrological projection. Third, the MEP model relies on a small set of equations, making it straightforward to implement, with minimal computational needs. Land surface models also offer a means of implementing coupled water and energy budget modelling, but contrary to the MEP model, these models are generally complex and computationally heavy. Last, the explicit partitioning of total terrestrial evaporation into evaporation and transpiration is also a strength of the MEP model given that the water feeding these two fluxes is drawn from different pools. Partitioning has particularly important implications for terrestrial evaporation modelling under water-limiting conditions (e.g. arid environments; Kurc and Small, 2004) or under changing land cover conditions (Huxman et al., 2005), and the MEP model thus offers a tool for better representing these conditions in hydrological modelling.

While these four highlighted features make the MEP model a promising approach for coupling water and energy budget modelling, certain limitations also need to be considered. First, the MEP model has mainly been tested with input data at a half-hourly time step. Although meteorological and climate data are increasingly available at a subdaily time step, with the continuing increase in the temporal resolution of reanalysis and climate projection datasets, additional tests would be needed to assess the applicability of the MEP model at a daily timescale. Second, soil water content is the key coupling variable between the MEP and hydrological models, which limits the choice of hydrological models which the MEP model can be coupled to. For the moment, the choice of a hydrological model appears to be limited to physical models, with the drawback being that these models are often computationally intensive as well as challenging in terms of parameterization. Indeed, the simulation time for the one-dimensional soil columns in this study was greater than an hour with the HGS-MEP model. HGS is a relatively complex model, and the MEP model has also been coupled to a soil moisture force–restore model with satisfying results (Huang and Wang, 2016; bias < 0.01 m³ m⁻³ and R²=0.81). As for conceptual hydrological models, further investigation is needed to assess the possibility of a coupling with the MEP model. In the absence of information on soil water content, a method would be needed to derive the specific humidity at the soil surface (q_ss) as well as the water stress factor (η_s) from the subsurface storage component of the conceptual model. Third, in its current form, the HGS-MEP model is still driven by dependent variables. For example, both net radiation and surface temperature are inputs to the model, although incoming long-wave radiation, a component of net radiation, is largely dependent on air temperature (used as a proxy for surface temperature T_ls). Moreover, soil surface temperature (T_ss) is an input to the model even though it is a function of the heat fluxes predicted by the model. In the present project, we have focused on the coupling of water fluxes between HGS and MEP, although in the future, the two models could also be more closely coupled, since thermal transport modelling is implemented within HGS (Brookfield et al., 2009). In addition to soil moisture information, the HGS model could supply information on soil surface temperature to the MEP model and thus eliminate its need as an input variable. For example, Huang and Wang (2016) simulated surface soil temperature and moisture using a force–restore model that relies on the MEP model to simulate the heat budget.

7 Conclusion

Using the MEP model of land surface fluxes, we proposed a simple approach to integrate energy budget modelling in hydrological models in order to improve the simulation of terrestrial evaporation. The MEP model requires six input variables (net radiation, soil surface temperature and specific humidity, leaf surface temperature and specific humidity and vegetation index) and ensures energy budget closure, which imparts a strong physical basis and avoids issues of oversensitivity to air temperature associated with certain PET models. We coupled the MEP model to HGS, an integrated surface and subsurface hydrologic model. Without calibration, the coupled HGS-MEP model performed well in simulating soil water content and terrestrial evaporation at three AmeriFlux sites with varying climates (semi-arid, Mediterranean and temperate). For both the simulation of daily soil moisture and terrestrial evaporation, HGS-MEP outperformed the stand-alone HGS model where, as defined by the Kristensen and Jensen (1975) model, terrestrial evaporation is derived from potential evaporation, which we computed using the Penman–Monteith equation. Overall, results indicate that, through a simple coupling procedure, the MEP model offers a physically constrained approach to simulate terrestrial evaporation in hydrological models. This approach may offer a tool for better assessing climate change impacts on water resources, although the predictive ability of the MEP model under environmental change, may it be natural (e.g. wildfires) or anthropogenic (e.g. land cover change and climate change), still needs to be assessed. This study focused on the simulation of vertical water fluxes, but to use HGS-MEP for flow simulation and projection, lateral fluxes will need to be considered in further work. Various routing models are available and could be used in conjunction with HGS-MEP to simulate lateral fluxes in a computationally efficient way. Finally, the present study focused on the application of the HGS-MEP model at snow-free sites, and the MEP model has undergone little testing in cold regions, with tests limited to the snow-free period (Wang et al., 2017). An MEP model for snow surfaces (Wang et al., 2014) is available and could also be integrated to hydrological models to allow energy budget modelling throughout the year in northern environments.