Bayesian performance evaluation of evapotranspiration models based on eddy covariance systems in an arid region

Evapotranspiration (ET) is a major component of the land surface process involved in energy fluxes and energy balance, especially in the hydrological cycle of agricultural ecosystems. While many models have been developed as powerful tools to simulate ET, there is no agreement on which model best describes the loss of water to the atmosphere. This study focuses on two aspects, evaluating the performance of four widely used ET models and identifying parameters, and the physical mechanisms that have significant impacts on the model performance. The four tested models are the Shuttleworth–Wallace (SW) model, Penman–Monteith (PM) model, Priestley–Taylor and Flint–Childs (PT–FC) model, and advection–aridity (AA) model. By incorporating the mathematically rigorous thermodynamic integration algorithm, the Bayesian model evidence (BME) approach is adopted to select the optimal model with half-hourly ET observations obtained at a spring maize field in an arid region. Our results reveal that the SW model has the best performance, and the extinction coefficient is not merely partitioning the total available energy into the canopy and surface but also including the energy imbalance correction. The extinction coefficient is well constrained in the SW model and poorly constrained in the PM model but not considered in PT–FC and AA models. This is one of the main reasons that the SW model outperforms the other models. Meanwhile, the good fitting of SW model to observations can counterbalance its higher complexity. In addition, the detailed analysis of the discrepancies between observations and model simulations during the crop growth season indicate that explicit treatment of energy imbalance and energy interaction will be the primary way of further improving ET model performance.

1 Introduction

Surface energy fluxes are an important component of Earth's global energy budget and a primary determinant of surface climate. Evapotranspiration (ET), as a major energy flux process for energy balance, accounts for about 60 %–65 % of the average precipitation over the surface of the Earth (Brutsaert, 2005). In agricultural ecosystems, more than 90 % of the total water losses are due to ET (Morison et al., 2008). Therefore, robust ET estimation is crucial to a wide range of problems in hydrology, ecology, and global climate change (Xu and Singh, 1998). In practice, mtuch of our understanding of how land surface processes and vegetation affect weather and climate is based on numerical modeling of surface energy fluxes and the atmospherically coupled hydrological cycle (Bonan, 2008). Several models are commonly used in agricultural systems to evaluate ET. The Penman–Monteith (PM) and Shuttleworth–Wallace (SW) models are physically sound and rigorous (Zhu et al., 2013) and thus widely used to simulate ET for seasonally varied vegetation. The models consider the relationships between net radiation, all kinds of heat flux (such as latent heat, sensible heat, and heat from soil and canopy), and surface temperature. The Priestley–Taylor and Flint–Childs (PT–FC) model (based on radiation) and the advection–aridity (AA) model (based on meteorological variables) have also been widely used because they only require a small number of ground-based measurements to set up the models (Ershadi et al., 2014).

Comparing the performance of the competing ET models and evaluating and understanding the discrepancies between simulations of the models and corresponding observed surface–atmosphere water flux remain challenging problems (Legates and McCabe, 1999). Both non-Bayesian analysis (Szilagyi and Jozsa, 2008; Vinukollu et al., 2011; Li et al., 2013; Ershadi et al., 2015) and Bayesian analysis have been used to evaluate the performance of ET models (Zhu et al., 2014; Chen et al., 2015; Liu et al., 2016; Zhang et al., 2017; Elshall et al., 2018; Samani et al., 2018;). Li et al. (2013) compared the ET simulations of the PM, SW and adjusted SW models under film-mulching conditions of maize growth in an arid region of China. They found that the half-hourly ET was overestimated by 17 % in the SW model. In contrast, the PM and adjusted SW models underestimated the daily ET by 6 % and 2 %, respectively. Therefore, the performances of PM and adjusted SW models are better than those of the SW model in their case study. Ershadi et al. (2014) evaluated the surface energy balance system (SEBS), PM, PT–JPL (JPL is the Priestley–Taylor Propulsion Laboratory Model; a modified Priestley–Taylor model), and AA models. Based on the average value of model efficiency (EF) and RMSE, the model ranking from worst to best was AA, PM, SEBS, and PT–JPL. Ershadi et al. (2015) also compared the response of the models to different formulations of aerodynamic and surface resistances with global FLUXNET data. Their results showed considerable variability in model performance among and within biome types. Currently, ET model selection and comparison have been still conducted using traditional error metrics. It is known that error metrics are not adequate for providing a reasonable result of model ranking for disregarding model complexity (Marshall et al., 2005; Samani et al., 2018). The focus of this study is to use a Bayesian approach to evaluate the performance of the PM, SW, PT–FC, and AA models, which is a novelty contribution of this study. In ET models, the land surface energy system is governed by presumably infinitely dimensional physics. However, considering the ET models to be finitely dimensional can be more precise by covering all relevant relations. Therefore, employing consistent criteria for model selection might be justified when the aim is to better understand the processes involved (Höge et al., 2018). When using consistent model selection, Bayesian model evidence (BME), also known as marginal likelihood, measures the average fit of model simulations to their corresponding observations over a model's prior parameter space. This feature enables BME to consider model complexity (in terms of the number of model parameters) for model performance evaluation. When comparing several alternative conceptual models, the model with the largest marginal likelihood is selected as the best model (Lartillot and Philippe, 2006). BME can thus be used for evaluating the model fit (over the parameter space) and for comparing alternative models. In previous studies, the Bayesian information criterion (BIC; Schwarz, 1978) and the Kashyap information criterion (KIC; Kashyap, 1982) have been used to approximate BME by using maximum likelihood theories to reduce the computational cost of evaluating BME (Ye et al., 2004). However, these approximations have theoretical and computational limitations (Ye et al., 2008; Xie et al., 2011; Schöniger et al., 2014), and a numerical evaluation (not a likelihood approximation) of BME is necessary, especially for complex models (Lartillot and Philippe, 2006). Lartillot and Philippe (2006) advocated the use of thermodynamic integration (TI) for estimating BME, also known as path sampling (Gelman and Meng, 1998; Neal, 2000), in order to avoid sampling solely in the prior or posterior parameter space. TI uses samples that are systematically generated from the prior to the posterior parameter space by conducting path sampling with several discrete power coefficient values (Liu et al., 2016). It is more numerically accurate than the generally used harmonic mean method (Xie et al., 2011).

Most applications of Bayesian methods have focused on the calibration of individual models, while the comparison of alternative models continues to be performed using traditional error metrics. More generally, Bayesian approaches to model calibration, comparison, and analysis have been used far less in the evaluation of ET models than in other areas of environmental science. In this study, the Bayesian approach is used to calibrate and evaluate the four ET models (PM, SW, PT–FC, and AA) based on an experiment over a spring maize field in an arid area of northwestern China from 3 June to 27 September 2014. The objectives of the study are as follows: (1) to calibrate ET model parameters using the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm (Vrugt et al., 2008, 2009), (2) to identify which parameters had a greater impact on the model performance and to explain why the selected optimal model performed best, (3) to evaluate the performance of the models using traditional error metrics and BME, and (4) to analyze discrepancies between model simulations and observation data in order to better understand model performance and identify ways to improve these models. We expect that the study will not only boost the development of model parameterization and model selection but also contribute to the improvement of the ET models.

2 Data and methodology

2.1 Description of the study area

The experiment of maize growth was conducted at the Daman Superstation, located in Zhangye, Gansu province, northwestern China. Daman oasis is located in the middle Heihe River basin, which is the second largest inland river basin in the arid region of northwestern China. The midstream area of the Heihe River basin is characterized by oases with irrigated agriculture and is a region that consumes a large amount of water for both domestic and agricultural uses. The annual average precipitation and temperature are 125 mm and 7.2^∘ (1960–2000), respectively. The annual accumulated temperature (>10^∘) is 3234^∘, and the annual average potential evaporation is about 2290 mm. The average annual duration of sunshine is 3106 h, with 148 frost-free days. The predominant soil type is silty-clay loam, and the depth of the frozen layer is about 143 mm. The study area is a typical irrigated agricultural region, and the major source of water is snowmelt from the Qilian Mountains. Maize and spring wheat are the principal crops grown in the region. Maize is generally sown in late April and harvested in mid-September and is planted with a row spacing of 40 cm and a plant spacing of 30 cm. The plant density is about 66 000 plants ha⁻¹ in the study area.

2.2 Measurements and data processing

Our data were collected from the field observation systems of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) project as described in Li et al. (2013). The observation period was from DOY (day of the year) 154 to DOY 270 in 2014. An open-path eddy covariance (EC) system was installed in a maize field, with the sensors at a height of 4.5 m. Maize is the main crop in the study region and thus covers sufficient planting area to set the EC measurements. The EC data were logged at a frequency of 10 Hz and then processed with an average time interval of 30 min. Sensible and latent heat fluxes were computed by the EC approach of Baldocchi (2003). Flux data measured by EC were controlled by traditional methods, including three-dimensional rotation (Aubinet et al., 2000); Webb–Penman–Leuning (WPL) density fluctuation correction (Webb et al., 1980); frequency response correction (Xu et al., 2014); and spurious data removal caused by rainfall, water condensation, and system failure. About 85 % of the energy balance closure was observed in the EC data (Liu et al., 2011).

Standard hydro-meteorological variables, including rainfall, air temperature, wind speed, and wind direction, were continuously measured at the heights of 3, 5, 10, 15, 20, 30, and 40 m above the ground. Soil temperature and moisture were measured at heights of 2, 4, 10, 20, 40, 80, 120, and 160 cm. Photosynthetically active radiation was measured at a height of 12 m. Net radiation, including downward, upward, and longwave radiation, was measured by a four-component net radiometer. An infrared thermometer was installed at a height of 12 m. The leaf area index (LAI) was measured approximately every 10 d during the growing season.

2.3 Model description

In this section, we summarize the mathematical definitions forming the basis of each of the four models. Appendix A contains a summary of the names and physical meanings of the model parameters.

2.3.1 Penman–Monteith (PM) model

The PM model can be formulated in the following way (Monteith, 1965):

$$\begin{array}{}\text{(1)}& \mathit{\lambda }E={\displaystyle \frac{\mathit{\epsilon }A+\left(\mathit{\rho }{C}_{p}t/\mathit{\gamma }\right)D_{\mathrm{a}}{g}_{\mathrm{a}}}{\mathit{\epsilon }+\mathrm{1}+g_{\mathrm{a}}/g_{\mathrm{s}}}},\end{array}$$

where $\mathit{\epsilon }=\mathrm{\Delta }/\mathit{\gamma }$, and A is defined as $A=R_{\mathrm{n}}-G$.

In the present study, g_a is parameterized in the way suggested by Leuning et al. (2008), and g_s is defined as

$$\begin{array}{}\text{(2)}& g_{\mathrm{s}}=g_{\mathrm{s}}^{\mathrm{c}}\left[{\displaystyle \frac{\mathrm{1}+\frac{\mathit{\tau }{g}_{\mathrm{a}}}{(\mathit{\epsilon }+\mathrm{1})g_{\mathrm{s}}^{\mathrm{c}}}\left[f-\frac{(\mathit{\epsilon }+\mathrm{1})(\mathrm{1}-f)g_{\mathrm{s}}^{\mathrm{c}}}{g_{\mathrm{a}}}\right]+\frac{g_{\mathrm{a}}}{\mathit{\epsilon }{g}_i}}{\mathrm{1}-\mathit{\tau }\left[f-\frac{(\mathit{\epsilon }+\mathrm{1})(\mathrm{1}-f)g_{\mathrm{s}}^{\mathrm{c}}}{g_{\mathrm{a}}}\right]+\frac{g_{\mathrm{a}}}{\mathit{\epsilon }{g}_i}}}\right],\end{array}$$

where 1−τ and τ are the fractions of the total available energy absorbed by the canopy and by the soil, τ=exp (K_aLAI), and g_i and $g_{\mathrm{s}}^{\mathrm{c}}$ are defined in Eqs. (3) and (4), respectively (Monteith, 1965),

$$\begin{array}{}\text{(3)}& g_i={\displaystyle \frac{A}{\left(\mathit{\rho }{C}_p/\mathit{\gamma }\right)D_{\mathrm{a}}}},\end{array}$$

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle }{g}_{\mathrm{s}}^{\mathrm{c}}=& {\displaystyle \frac{g_{\max }}{K_q}}In\left[{\displaystyle \frac{Q_{\mathrm{h}}+Q_{\mathrm{50}}}{Q_{\mathrm{h}}\exp \left(-K_q\mathrm{LAI}\right)+Q_{\mathrm{50}}}}\right]\\ {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+D_{\mathrm{a}}}}/D_{\mathrm{50}}\right]f\left(\mathit{\theta }\right),\end{array}\end{array}$$

where f(θ) represents water stress and is expressed as

$$\begin{array}{}\text{(5)}& f\left(\mathit{\theta }\right)=\left\{\begin{array}{ll}\mathrm{1}& \mathit{\theta }>{\mathit{\theta }}_a,\\ {\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_b}{{\mathit{\theta }}_a-{\mathit{\theta }}_b}}& {\mathit{\theta }}_b<\mathit{\theta }<{\mathit{\theta }}_{a}l,\\ \mathrm{0}& \mathit{\theta }<{\mathit{\theta }}_{b}l,\end{array}\right.\end{array}$$

and θ_a is set as θ_a=0.75 θ_b. Aerodynamic conductance g_a is calculated as

$$\begin{array}{}\text{(6)}& g_{\mathrm{a}}={\displaystyle \frac{k^{\mathrm{2}}{u}_{\mathrm{m}}}{\ln \left[\left(z_{\mathrm{m}}-d\right)/z_{\mathrm{0}\mathrm{m}}\right]\ln \left[\left(z_{\mathrm{m}}-d\right)/z_{\mathrm{0}\mathrm{v}}\right]}},\end{array}$$

where the quantities d, z_0 m, and z_0 v are calculated using $d=\mathrm{2}h/\mathrm{3}$, z_0 m=0.123 h, and $z_{{\mathrm{0}}_{\mathrm{v}}}=\mathrm{0.1}{z}_{\mathrm{0}\mathrm{m}}$ (Allen et al., 1998).

2.3.2 Shuttleworth–Wallace (SW) model

The SW model comprises a one-dimensional model of plant transpiration and a one-dimensional model of soil evaporation. The two terms are calculated by the following equations:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }\mathrm{ET}=\mathit{\lambda }E+\mathit{\lambda }T=C_{\mathrm{s}}{\mathrm{ET}}_{\mathrm{s}}+C_{\mathrm{c}}{\mathrm{ET}}_{\mathrm{c}},\text{(8)}& {\displaystyle }& {\displaystyle }{\mathrm{ET}}_{\mathrm{s}}={\displaystyle \frac{\mathrm{\Delta }A+\left\{\mathit{\rho }{C}_p\left(e_{\mathrm{s}}-e_{\mathrm{a}}\right)-\mathrm{\Delta }{r}_{\mathrm{a}}^{\mathrm{s}}\left(A-A_{\mathrm{s}}\right)\right\}/\left(r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{s}}\right)}{\mathrm{\Delta }+\mathit{\gamma }\left\{\mathrm{1}+r_{\mathrm{s}}^{\mathrm{s}}/\left(r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{s}}\right)\right\}}},\text{(9)}& {\displaystyle }& {\displaystyle }{\mathrm{ET}}_{\mathrm{c}}={\displaystyle \frac{\mathrm{\Delta }A+\left\{\mathit{\rho }{C}_p\left(e_{\mathrm{s}}-e_{\mathrm{a}}\right)-\mathrm{\Delta }{r}_{\mathrm{a}}^{\mathrm{c}}{A}_{\mathrm{s}}\right\}/\left(r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{c}}\right)}{\mathrm{\Delta }+\mathit{\gamma }\left\{\mathrm{1}+r_{\mathrm{s}}^{\mathrm{c}}/\left(r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{c}}\right)\right\}}},\end{array}$$

where the available energy input above the soil surface is defined as $A_{\mathrm{s}}=R_{\mathrm{ns}}-G$.

R_ns can be calculated using the Beer's law relationship:

$$\begin{array}{}\text{(10)}& R_{\mathrm{ns}}=R_{\mathrm{n}}\exp \left(-K_{\mathrm{a}}\mathrm{LAI}\right).\end{array}$$

The coefficients C_s and C_c are obtained as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{s}}={\left\{\mathrm{1}+R_{\mathrm{s}}{R}_{\mathrm{a}}/R_{\mathrm{c}}\left(R_{\mathrm{s}}+R_{\mathrm{a}}\right)\right\}}^{-\mathrm{1}},\text{(12)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{c}}={\left\{\mathrm{1}+R_{\mathrm{c}}{R}_{\mathrm{a}}/R_{\mathrm{s}}\left(R_{\mathrm{c}}+R_{\mathrm{a}}\right)\right\}}^{-\mathrm{1}},\end{array}$$

where

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{a}}=(\mathrm{\Delta }+\mathit{\gamma })r_{\mathrm{a}}^{\mathrm{a}},\text{(14)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{s}}=(\mathrm{\Delta }+\mathit{\gamma })r_{\mathrm{a}}^{\mathrm{s}}+\mathit{\gamma }{r}_{\mathrm{s}}^{\mathrm{s}},\text{(15)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{c}}=(\mathrm{\Delta }+\mathit{\gamma })r_{\mathrm{a}}^{\mathrm{c}}+\mathit{\gamma }{r}_{\mathrm{s}}^{\mathrm{c}}.\end{array}$$

Soil surface resistance is expressed as

$$\begin{array}{}\text{(16)}& r_{\mathrm{s}}^{\mathrm{s}}=\exp \left(b_{\mathrm{1}}-b_{\mathrm{2}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{s}}}}\right).\end{array}$$

In this study, we consider the reciprocal of bulk stomatal resistance, known as canopy conductance. The calculation of $g_{\mathrm{s}}^{\mathrm{c}}$ is the same as in the PM model. The two aerodynamic resistances ($r_{\mathrm{a}}^{\mathrm{a}}$ and $r_{\mathrm{a}}^{\mathrm{s}}$) and the boundary layer resistance ($r_{\mathrm{a}}^{\mathrm{c}}$) are modeled following the approach proposed by Shuttleworth and Gurney (1990).

2.3.3 Priestley–Taylor and Flint–Childs (PT–FC) model

The Priestley–Taylor model (Priestley and Taylor, 1972) was introduced to estimate evaporation from an extensive wet surface under conditions of minimum advection (Stannard, 1993; Sumner and Jacobs, 2005). The ET is expressed as

$$\begin{array}{}\text{(17)}& \mathit{\lambda }\mathrm{ET}={\mathit{\alpha }}_{\mathrm{PT}}{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}\left(R_{\mathrm{n}}-G\right),\end{array}$$

where α_PT is a unitless coefficient. The Priestley–Taylor model was modified by Flint and Childs (1991) in order to scale the Priestley–Taylor potential ET to actual ET for nonpotential conditions (hereafter the PT–FC model):

$$\begin{array}{}\text{(18)}& \mathit{\lambda }\mathrm{ET}=\mathit{\alpha }{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}\left(R_{\mathrm{n}}-G\right),\end{array}$$

where α is as a function of the environmental variables, which could be related to any process that limits ET (e.g., soil hydraulic resistance, aerodynamic resistance, and stomatal resistance); however, only soil moisture status was considered for simplifying ET estimation in the PT–FC model (Flint and Childs, 1991). In this model, α is defined as

$$\begin{array}{}\text{(19)}& \mathit{\alpha }={\mathit{\beta }}_{\mathrm{1}}\left[\mathrm{1}-\exp \left(-{\mathit{\beta }}_{\mathrm{2}}\mathrm{\Theta }\right)\right],\end{array}$$

where $\mathrm{\Theta }=\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$.

2.3.4 Advection–aridity (AA) model

The AA model was first proposed by Brutsaert and Stricker (1979) and further improved by Parlange and Katul (1992). The model relies on the feedback between actual ET (λET) and potential ET, which assumes that actual potential ET should converge to wet surface ET at wet surface conditions. Its general form is

$$\begin{array}{}\text{(20)}& \mathit{\lambda }\mathrm{ET}=\left(\mathrm{2}{\mathit{\alpha }}_{\mathrm{PT}}-\mathrm{1}\right){\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}\left(R_{\mathrm{n}}-G\right)-{\displaystyle \frac{\mathit{\gamma }}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{\mathit{\rho }\left(q^{*}-q\right)}{r_{\mathrm{a}}}},\end{array}$$

where α_PT is the Priestley–Taylor coefficient, usually taken as 1.26 (Priestley and Taylor, 1972), and r_a is similar to that used for the Penman–Monteith model (Brutsaert and Stricker, 1979; Brutsaert, 2005; Ershadi et al., 2014). This model is based mainly on meteorological variables and does not require any information related to soil moisture, canopy resistance, or other measures of aridity (Ershadi et al., 2014). In this study, we changed α_PT to α, which is calculated using the same equation as in the PT–FC model.

2.4 BME Estimation

The Bayesian model evidence (BME) of a model, M, is defined as (Schöniger et al., 2014)

$$\begin{array}{}\text{(21)}& \mathrm{BME}=p\left(D\right|M)=\int p(D|\mathit{\theta },M)p\left(\mathit{\theta }\right|M)d\mathit{\theta },\end{array}$$

where D is observed or estimated data, θ is the vector of parameters associated with model M, and p(θ|M) is the prior density of θ under model M; $p\left(D\right|\mathit{\theta },M)$ is the joint likelihood of model M and its parameters θ. Estimating BME using power posterior estimators such as thermodynamic integration (Lartillot and Philippe, 2006) depends mainly on the calculation of the marginal likelihood p(D|M). The main idea of power posterior sampling is to define a path that links the prior to the unnormalized posterior. Thus, using an unnormalized power posterior density,

$$\begin{array}{}\text{(22)}& q_{\mathit{\beta }}\left(\mathit{\theta }\right)=p\left(D\right|\mathit{\theta },M{)}^{\mathit{\beta }}p\left(\mathit{\theta }\right|M),\end{array}$$

the power coefficient β∈[0, 1] is a scalar parameter for discretizing a continuous and differentiable path linking two unnormalized power posterior densities. The unnormalized power posterior density q_β(θ) in Eq. (22) uses the normalizing constant Z_β to yield the normalized power posterior density,

$$\begin{array}{}\text{(23)}& p_{\mathit{\beta }}\left(\mathit{\theta }\right)={\displaystyle \frac{q_{\mathit{\beta }}\left(\mathit{\theta }\right)}{Z_{\mathit{\beta }}}},\end{array}$$

such that

$$\begin{array}{}\text{(24)}& Z_{\mathit{\beta }}=\int q_{\mathit{\beta }}\left(\mathit{\theta }\right)d\mathit{\theta }.\end{array}$$

The above integral takes a simplified form by the potential

$$\begin{array}{}\text{(25)}& U\left(\mathit{\theta }\right)={\displaystyle \frac{\partial \ln q_{\mathit{\beta }}\left(\mathit{\theta }\right)}{\partial \mathit{\beta }}}.\end{array}$$

Thus, the integral can be directly estimated in the following way:

$$\begin{array}{}\text{(26)}& p\left(D\right|M)={\displaystyle \frac{Z_{\mathrm{1}}}{Z_{\mathrm{0}}}}=\exp \left\{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{E}_{\mathit{\theta }}\left[\ln p\right(D|\mathit{\theta },M)]d\mathit{\beta }\right\}.\end{array}$$

The one-dimensional integral with respect to β is evaluated by using numerical methods by discretizing β into a set of β_k values. Since there is no theoretical method for selecting β_k values (Liu et al., 2016), we determined these values using an empirical but straightforward method. Following Xie et al. (2011), a schedule of the power posterior coefficients β_k is generated by

$$\begin{array}{}\text{(27)}& {\mathit{\beta }}_k=(k/K{)}^{\mathrm{1}/\mathit{\epsilon }},\end{array}$$

for k=0, 1, 2 …, K. Using ε=0.3 and K=20 is a reasonable initial choice. By using the trapezoidal rule of numerical integration, Eq. (26) is evaluated via

$$\begin{array}{}\text{(28)}& p\left(D\right|M)=\exp \left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{y}_{\mathit{\beta }}d\mathit{\beta }\right)=\exp \left(\sum _{k=\mathrm{0}}^K{r}_{\mathrm{TI},k}\right),\end{array}$$

such that

$$\begin{array}{}\text{(29)}& r_{\mathrm{TI},k}=\left({\mathit{\beta }}_k-{\mathit{\beta }}_{k-\mathrm{1}}\right)\left[{\displaystyle \frac{y_k+y_{k-\mathrm{1}}}{\mathrm{2}}}\right],\end{array}$$

and

$$\begin{array}{}\text{(30)}& y_k=E_{{\mathit{\beta }}_k}\left[\ln p\left(D|{\mathit{\theta }}_k,M\right)\right]={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\ln p\left(D|{\mathit{\theta }}_{k,i},M\right),\end{array}$$

where n is the number of random samples of θ_k corresponding to β_k, and θ_k,i is the ith sample.

The random samples, θ_k,i, are drawn by using the Markov chain Monte Carlo (MCMC) method implemented in the DREAM code. See Appendix B for further details on Bayesian inference and the DREAM algorithm. In the DREAM-based calculation, the Metropolis acceptance ratio is ${\mathit{\alpha }}_k=\min (\mathrm{1},[{\mathit{\alpha }}_{k,\mathrm{power}\text{-}\mathrm{posterior}}{\mathit{\alpha }}_{k,\mathrm{prior}}\left]\right)$ with the power posterior ratio given by ${\mathit{\alpha }}_{k,\mathrm{power}\text{-}\mathrm{posterior}}=({\mathit{\alpha }}_{k,\mathrm{posterior}}{)}^{{\mathit{\beta }}_k}$. The prior probability ratio ${\mathit{\alpha }}_{k,\mathrm{prior}}=\mathrm{Pr}\left({\mathit{\theta }}_{k,\mathrm{new}}\right|M)/\mathrm{Pr}({\mathit{\theta }}_{k,\mathrm{old}}\left|M\right)$ is the ratio of the probability of the newly proposed sample θ_k,new and the probability of the previously accepted sample θ_k,old. The posterior probability ratio ${\mathit{\alpha }}_{k,\mathrm{posterior}}=L\left(D\right|{\mathit{\theta }}_{k,\mathrm{new}},M)/L(D|{\mathit{\theta }}_{k,\mathrm{old}}$, M) is the likelihood ratio of samples θ_k,new and θ_k,old, and β_k  is the power posterior coefficient. Thus, to use the DREAM algorithm to sample any power posterior distribution, the regular Metropolis acceptance ratio $\mathit{\alpha }=\min (\mathrm{1},[{\mathit{\alpha }}_{\mathrm{posterior}}{\mathit{\alpha }}_{\mathrm{prior}}\left]\right)$ is changed to ${\mathit{\alpha }}_k=\min (\mathrm{1},[{\mathit{\alpha }}_{k,\mathrm{power}\text{-}\mathrm{postrior}}{\mathit{\alpha }}_{k,\mathrm{prior}}\left]\right)$ in DREAM.

2.5 Traditional statistical metrics of evaluating model performance

The traditional error metrics for evaluating model performance include R² and slope (correlation-based measures), the index of agreement (IA) and EF (relative error measures), and the root-mean-square error (RMSE) and mean bias error (MBE; Poblete-Echeverria and Ortega-Farias, 2009). The definitions of the listed metrics are

$$\begin{array}{}\text{(31)}& {\displaystyle }& {\displaystyle }\mathrm{IA}=\mathrm{1}-{\displaystyle \frac{\sum _{t=\mathrm{1}}^n\left[O\right(t)-M(t){]}^{\mathrm{2}}}{\sum _{t=\mathrm{1}}^n\left[\right|O\left(t\right)-\stackrel{\mathrm{‾}}{O\left(t\right)}|+|O\left(t\right)-\stackrel{\mathrm{‾}}{M\left(t\right)}|{]}^{\mathrm{2}}}},\text{(32)}& {\displaystyle }& {\displaystyle }\mathrm{EF}=\mathrm{1}-{\displaystyle \frac{\sum _{t=\mathrm{1}}^n\left[O\right(t)-M(t){]}^{\mathrm{2}}}{\sum _{t=\mathrm{1}}^n\left[O\right(t)-\stackrel{\mathrm{‾}}{O\left(t\right)}{]}^{\mathrm{2}}}},\text{(33)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{t=\mathrm{1}}^n\left[O\right(t)-M(t){]}^{\mathrm{2}}},\text{(34)}& {\displaystyle }& {\displaystyle }\mathrm{MBE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{t=\mathrm{1}}^n\left[O\right(t)-M(t\left)\right],\end{array}$$

where O(t) is the observation and $\stackrel{\mathrm{‾}}{O\left(t\right)}$ is the mean observation at time t, M(t) is the modeled value and $\stackrel{\mathrm{‾}}{M\left(t\right)}$ is the mean modeled value estimated by the posterior median parameter values, and n is the total number of the observed values.

Figure 1 Trace plots of the G–R statistic using DREAM for the PM model (a) and (b) the SW model. Different parameters are coded with different colors. The dashed line denotes the default threshold used to diagnose convergence to a limiting distribution.

Table 1 Prior distributions and parameter limits for the PM, SW, PT–FC and AA models. The values are derived from the literature.

3 Results

3.1 Parameter estimation

The PM model has five parameters, g_max, D₅₀, Q₅₀, K_q, and K_a; the SW model has seven parameters – the five used in the PM model and parameters b₁ and b₂. The PT–FC and AA models each include two parameters, denoted by β₁ and β₂ (Table 1). The prior probability density of each parameter is specified as an uniform distribution with the ranges listed in Table 1. A total of 50 000 realizations were generated with the DREAM algorithm, which was used to estimate the posterior probability density function of each parameter with the calibration period data from DOY 154 to DOY 202. In the calculations, the chain number, N, was equal to the number of parameters in the associated model. Therefore, N is equal to 5, 7, 2, and 2 for the PM, SW, PT–FC, and AA models, respectively. For each model, the first 10 000 samples were discarded as burn-in data, and the remaining 40 000 samples were used for calibration. In total, 40 000×N realizations were used to set up posterior density functions for each model. To illustrate the efficiency and convergence of DREAM for the ET models, Fig. 1 shows the trace plots of the G–R (Gelman and Rubin, 1992) statistic for each of the different parameters in the PM and SW models using a different color. The algorithm required about 8000 generations to make the G–R statistic close to 1.0 for the two models. The acceptance rates for the PM and SW models were about 15.3 % and 18.9 %, respectively.

Figure 2 (a–e), (f–l), (m–n), and (o–p) show histograms for the PM (black), SW (cyan), PT–FC (magenta), and AA (orange) models, respectively. These histograms are constructed from all chains for each model, and a total of 40 000×N realizations are simulated using DREAM. The x axes represent the prespecified limits of the parameters.

Table 2 Maximum likelihood estimates (MLEs), mean estimates, 95 % high-probability intervals (lower limit, upper limit).

^a PT–FC model; ^b AA model.

Histograms of the DREAM-derived marginal distributions of the parameters are presented in Fig. 2 and summarized in Table 2 by maximum likelihood estimates (MLEs), posterior medians, and 95 % probability intervals. Figure 2a–e, f–l, m–n, and o–p show histograms of the PM, SW, PT–FC, and AA models, respectively. Parameter g_max (Fig. 2a) in the PM model; parameters g_max, K_a, and b₂ (Fig. 2f, j, l) in the SW model; and parameter β₁ (Fig. 2m) in the PT–FC model and AA model (Fig. 2o) were well constrained and occupied a relatively small range. These parameters displayed a unimodal distribution and appeared approximately Gaussian. In contrast, the distributions of the other parameters differed significantly from a Gaussian distribution, as shown by the corresponding histograms. The distributions of all but one of these parameters concentrated most of the probability mass at their upper limits. The exception was parameter b₁ for the SW model (Fig. 2k), which clearly does not follow a normal distribution, with most of the mass concentrated in the lower bounds. In contrast, Q₅₀ was not only poorly constrained (Fig. 2g) but was also the upper-edge-hitting parameter in the SW model. Moreover, the corresponding distributions of the same parameter in different models were slightly different. For example, the mean of g_max in the PM model (0.04 mm s⁻¹) was less than that in the SW model (0.01 mm s⁻¹; Fig. 2a and f; Table 2) except that D₅₀ in the PM and SW models and β₂ in the PT–FC and AA models exhibited similar regions. It is interesting to observe that the distribution of K_a in the PM model (Fig. 2e) has a truncated distribution with highest probability mass at the upper bound, whereas the distribution of K_a in the SW model (Fig. 2j) tends to become approximately normal. Overall, the marginal posterior probability density function of most of the individual parameters occupied only a relatively small region compared with the uniform prior distributions and exhibited relatively large uncertainty reduction.

3.2 Performance of the models

The performance of each of the four ET models was evaluated over the course of the whole season in 2014. The calibrated parameters of the four models were used, and individual ET models were run to estimate the half-hourly λET values. Table 3 summarizes the statistical results for the performance of the models using the regression line slope, R², RMSE, MBE, IA, and EF. The regressions between measured and modeled λET values and the MBE are shown in Figs. 3 and 4, respectively.

Figure 3 Regressions between measured and modeled half-hourly ET values produced by different models from DOY 154 to DOY 270: (a) PM, (b) SW, (c) PT–FC, and (d) AA. The regressions are Y=0.99X (R²=0.76), Y=1.05X (R²=0.82), Y=0.91X (R²=0.75), and Y=0.92X (R²=0.75) for the PM, SW, PT–FC, and AA models, respectively.

Figure 4 Mean bias error (MBE) of predicted and observed ET values for (a) PM, (b) SW, (c) PT–FC, and (d) AA models from DOY 154 to DOY 270. Parameters used for prediction are estimated by DREAM with the dataset for the calibration period from DOY 154 to DOY 202.

Table 3 Slope and coefficient of determination (R²) of regression between measured and modeled half-hourly evapotranspiration values, and statistics of root mean square error (RMSE), mean bias error (MBE), index of agreement (IA), model efficiency (EF) and Logarithm of BME for the four ET models.

In general, the four models produced slightly better fits to the measured λET for all the seasons with R² larger than 0.75 (Fig. 3). However, obvious discrepancies in the predictions made by the models were detected by comparing measured and modeled λET. According to the regression line slope and MBE, the PM model overestimated ET by 1 %, with an MBE of −9.52 W m⁻², and the SW model overestimated ET by 5 %, with a relatively higher MBE of −19.07 W m⁻² compared to the PM model. The PT–FC and AA models tended to underestimate λET by 9 % and 8 %, with an MBE of 25.42 and 23.29 W m⁻², respectively. From a comparison between the slope and MBE, the PM model performance was higher than that of the other three models, with a slope almost equal to 1 and a relatively lower MBE. The SW model was ranked second, while the performance of the AA model was slightly higher than that of the PT–FC model. However, if R², RMSE, IA, and EF were used to evaluate performance, the SW model had the best overall performance, with R²=0.83, RMSE = 76.34 W m⁻², IA = 0.95, and EF = 0.79. The second-best model was the PM model, the PT–FC model was ranked third, and the AA model ranked fourth. Based on the analysis of these traditional error metrics, the PT–FC and AA models yielded similar results. The observed and modeled λET for the four ET models were tightly grouped along the regression lines (Fig. 3), and the PT–FC and AA models had similar modeled ET values with a similar degree of point scattering along the regression lines (Fig. 3 and d).

Figure 4 shows that large seasonal variations arise in the MBE for the four ET models. From the variations in the MBE, the estimated λET values for all models were generally lower than the measured values before the early jointing stage of maize growth (DOY 154–177; left dashed line) and after the late maturity stage (DOY 256–265; right dashed line) with the corresponding LAI < 2.5 m² m⁻². More positive MBE values for the PT–FC and AA models after the late maturity stage indicate their underestimated performances; however, these estimations appeared to be even more consistent, with a symmetrical scattering of points along the 0–0 line (Fig. 4c and d) during DOY 177–256 and with LAI > 2.5 m² m⁻².

3.3 Comparison of the models using BME

Since there is currently no theoretical method for selecting power posterior β values, we determined these values using empirical but straightforward methods. For any power coefficient of β∈[0, 1], a sample was drawn from the distribution p_β (Eq. 23) through running DREAM. Although adding more β_k values might improve the BME estimation, this was not done because of the computational cost. For each β_k value, at least 150 000 DREAM simulations were large enough to ensure convergence. Figure 5 shows the evolution of ln p(D|θ, M) for the four models as a function of β for a dataset covering the entire period. The BME for the SW model was substantially larger than that for the other three models, and the BME for the AA model was the smallest. The BME-based model ranking (from the best to the worst) is SW, PM, PT–FC, and AA. The PT–FC and AA models, which consist of the same number of parameters, had similar potential patterns of evolution with respect to the coefficient β_k. The results illustrate that with the addition of parameters, the model complexity and the model performance are both increased.

Figure 5 Variation of the mean posterior expectation of the potential y_k with β_k for the PM, SW, PT–FC, and AA models.

4 Discussion

4.1 Parameter uncertainty analysis

With regard to the efficiency of the DREAM algorithm, the acceptance rates of the PM (15.3 %) and SW (18.9 %) models were much higher than those obtained by some MCMC algorithms that have been used in previous studies (Sadegh and Vrugt, 2014). The posterior parameter bounds exhibit a larger reduction using the DREAM algorithm compared with other studies using the Metropolis–Hastings algorithm. This demonstrates that DREAM could efficiently handle problems involving high dimensionality, multimodality, and nonlinearity.

The results showed that the assumed prior uncertainty ranges from most parameters in the four models were significantly reduced. This indicates that the observed ET data contained sufficient information for estimating these parameters. Surface conductance g_s and modeled ET in the PM model are relatively insensitive to Q₅₀, D₅₀, and K_q. Hence, these parameters could not be well constrained, and further relaxing the ranges for these parameters could not result in physically realistic behavior of the model. The calculation of $g_{\mathrm{s}}^{\mathrm{c}}$ in the SW model is the same as in the PM model, and thus $g_{\mathrm{s}}^{\mathrm{c}}$ and modeled ET in the SW model are also insensitive to the parameters of Q₅₀, D₅₀, and K_q. Therefore, these three parameters were also not well constrained in the SW model. In addition, the uncertainties present in the edge-hitting parameters may be the outcome of model biases or EC-measured ET data errors or the characteristic timescale of parameters governing the processes affecting ET being not exactly on the order of half-hours (Braswell et al., 2005). For example, Q₅₀ and D₅₀ govern changes in visible radiation flux and the humidity deficit at which stomatal conductance is half its maximum value, respectively, and these parameters may change over a shorter or longer timescale than half-hours.

The ecophysiological parameter g_max is a variable in the $g_{\mathrm{s}}^{\mathrm{c}}$ equation in both the PM and SW models, but this parameter is sensitive to $g_{\mathrm{s}}^{\mathrm{c}}$ and has a significant impact on the evaluated ET. Its effect is relatively independent compared to the other meteorological parameters in the models, and therefore this parameter was well specified in the PM and SW models. The posterior mean value of g_max (0.04 m s⁻¹) in the PM model from our study was close to that (0.05 m s⁻¹) reported in northwestern China (Li et al., 2013; Zhu et al., 2014), but g_max (0.01 m s⁻¹) in the SW model was less than the reported value. Parameter β₁ was well constrained in the PT–FC and AA models because it was relatively independent and did not directly relate to other observed variables.

Parameter K_a implicitly appears in the surface conductance equation (Eq. 2) in the PM model, and K_a is insensitive to g_s and modeled ET (Leuning et al., 2008). In contrast, K_a is contained in the equation of net radiation flux into the substrate (Eq. 10) in the SW model. This parameter can explicitly partition the total available energy into that absorbed by the canopy and by the soil in the SW model. An analysis of Eq. (10) found that the variation of K_a could not only account for the extinction effect but also correct the energy forcing data errors. This also meant that the estimated value of K_a using calibration data was actually not just the true extinction coefficient but also included the energy imbalance correction in the SW model. From this analysis, we could see that K_a not only involved the distribution of energy between the canopy and the soil surface but also the energy imbalance. Therefore, parameter K_a has a great influence on the performance of the SW model. This is why K_a is poorly constrained in the PM model but well constrained in the SW model. To further illustrate the insights regarding the influence of parameter K_a on the performance of the SW model, we calibrated the SW model again and reran the model with a constant value of K_a. The results showed a significant reduction in model performance when K_a was held constant. This implied that the main reason for the SW model outperforming the PM model in our study was not only the more physically rigorous structure of the SW model but also the key parameter K_a being well constrained in the SW model.

In general, parameters related to soil surface resistance in the SW model were well evaluated, while parameters related to canopy surface resistance in PM and SW models were poorly estimated. Therefore, using a reliable canopy surface resistance equation in the ET model was crucial for improving its performance. In addition, in our study, the traditional approach was used to quantify the uncertainty, which assumed that the uncertainty mainly arose because of the parameter uncertainty. However, this method cannot explicitly consider errors in the input data and model structural inadequacies. This is unrealistic for real applications, and it is desirable to develop a more reliable inference method to treat all sources of uncertainty separately and appropriately (Vrugt et al., 2008). Moreover, simultaneous direct measurement by the micro-lysimeter of sap flow and daily soil evaporation will further help to constrain the model parameters.

4.2 Evaluation and selection of the models

In this study, the traditional statistical measures and BME were chosen to evaluate and compare the performance of four ET models. From the respective composition of these measures, the statistical measures can be divided into residual-based metrics (such as regression slope and MBE) and squared-residual-based measures (such as R², RMSE, IA, and EF). The rankings of the models obtained using the same type of metric (residual-based or squared-residual-based) are similar. Slope and MBE, for example, which are both residual-based measures, produce identical rankings. However, the rankings produced by metrics of different types are not the same. For example, the PM model outperforms the SW model according to the residual-based metrics, but the performance of the PM model is worse than the SW model based on the squared-residual-based measures. The comparative analysis shows consistency between BME and the squared-residual-based metrics (hence the residual-based metrics disagreed with the BME measures). This reveals that the more complex SW model is the best model based on BME and squared-residual-based statistics. The rank order of overall performance of the models from best to worst is SW, PM, PT–FC, and AA.

Previous studies showed that BME evaluated by TI provided estimates similar to the true values and selected the true model if the true model was included within the candidate models (Marshall et al., 2005; Lartillot and Philippe, 2006). Meanwhile, some have argued that Bayesian analysis would choose the simplest model (Jefferys and Berger, 1992; Xie et al., 2011) because of the best trade-off between good fit with the data and model complexity (Schöniger et al., 2014). In this case, the most complex SW model had the highest BME and was chosen as the model with the best performance. This probably resulted from the fact that the complex SW model is indeed the most reliable model among the alternative ET models and can provide a good fit to justify its higher complexity. The SW model is a two-layer model and simulates soil evaporation and plant transpiration separately, whereas the PM model is a single-layer model in which the plant transpiration and soil evaporation cannot be separated (Monteith, 1965). The PT–FC model is a simplified version of the PM model and only requires meteorological and radiation information (Priestley and Taylor, 1972), whereas the AA model only relies on the feedback between actual ET and potential ET (Brutsaert and Stricker, 1979).

The results indicate that the squared-residual-based measures yielded the same rank order as the BME consistently, which makes the squared-residual-based metrics seem to identify a reasonable rank order. However, this has not been the general case, since the error metrics and BME belong to different types of model selection and because there are differences in the behavior and optimality of the two types of model selection. BME is a consistent model selection that tries to identify which of the models produced the observed data. Conversely, nonconsistent model selection uses the available data to estimate which of the models might be best in predicting future data. In fact, the error metrics are essentially nonparsimonious model selection, which is a special case of nonconsistent model selection. The simple traditional statistical measures were known to usually provide a biased view of the efficacy of a model (Kessler and Neas, 1994; Legates and McCabe, 1999), where only the goodness of fit is used for rating models without penalizing the model complexity, thus lacking consistency for the selected model (Höge et al., 2018). In addition, sensitivity to outliers is associated with these metrics and leads to relatively high values due to the squaring of the residual terms (Willmott, 1981). Furthermore, these traditional statistical metrics ignore the priors, which are in fact used in Bayesian analysis. The PT–FC and AA models provide identical estimates of R² and IA. This is most likely because both models had the same dimension and a similar model structure. Marshall et al. (2005) argued that EF would provide an incorrect conclusion, and Samani et al. (2018) suggested that RMSE would select the complex model as the best performing model. As for the slope and MBE, the rankings produced by these residual-based metrics were in obvious disagreement with the one based on BME. Part of the lower simulation values could be counterbalanced by the higher values of that in the slope and MBE methods; thus these criteria provide an erroneous and unreliable evaluation of the models. Therefore, the squared-residual-based and residual-based measures were not certain in providing reasonable results in terms of model ranking. The consistency between BME and the squared-residual-based metrics only indicates that the optimal model evaluated by BME would also provide the best predictions, and thus consistent model selection should also be asymptotically efficient (Leeb and Pötscher, 2009; Shao, 1997).

4.3 Analysis of model–data mismatch

Conceptual and structural inadequacies of the hydrological model together with measurement errors of the model input (forcing) and output (calibration) data introduce errors in the estimated parameters and model simulations (Laloy et al., 2015). Hydrological systems are indeed heavily input-driven, and errors in forcing data can dramatically impair the quality of calibration results and model output (Bardossy and Das, 2008; Giudice et al., 2016). Measurement errors occur for a variety of reasons, including unreasonable gap filling in rainy days; dew and fog; inadequate areal coverage of point-scale soil water measurement; mechanical limitations of the EC system; and inaccurate measurements of wind speed, soil water, radiation, and vapor pressure deficit. The ET process is described using equations that can only capture parts of the complex natural processes, and any ET model is an inherent simplification of the real system. These inadequacies can thus lead to biased parameters and implausible predictions.

In our study, the results indicated that the PM and SW models overestimated the half-hourly ET compared to the measured ET. Several studies also indicated that ET was overestimated by the PM model (Fisher et al., 2005; Ortega-Farias et al., 2006; Li et al., 2015) and the SW model (Li et al., 2013, 2015; Zhang et al., 2008). Possible reasons for the inaccurate estimates included the following: (1) anisotropic turbulence with weak vertical and strong horizontal fluctuation leads to energy imbalance. The total turbulent heat flux was lower by ∼10 %–30 % compared to the available energy in many land surface experiments (Tsvang et al., 1991; Beyrich et al., 2002; Oncley et al., 2007; Foken et al., 2010) and influx networks (Franssen et al., 2010). Liang et al. (2017) also showed an energy imbalance result in the semiarid area in China and indicated that the energy balance closure ratio ranged from 0.52 to 0.90 during the day, whereas it was about 0.25 at night. However, the measured ET only included vertical flux and not horizontal flux, leading to the measured ET being lower than that of ET predicted by the PM and SW models using the available energy. (2) The absence of a mechanistic representation of the physiological response to plant hydrodynamics makes it difficult for the available ET models to resolve the dynamics of intradaily hysteresis, producing patterns of diurnal error, while the imbalance or lack of between-leaf water demand and soil water supply imposes hydrodynamic limitations on stomatal conductance (Thomsen et al., 2013; Zhang et al., 2014; Matheny et al., 2014). Li et al. (2015) also concluded that neglecting the restrictive effect of the soil on water transport in empirical canopy resistance equations can result in large errors in the partial canopy stage. However, these equations can estimate ET accurately under the full canopy stage (Alves and Pereira, 2000; Katerji and Rana, 2006; Katerji et al., 2011; Rana et al., 2011). Li et al. (2015) showed that the PM model combined with canopy resistance overestimated maize ET during the partial and dense canopy stages by 16 % and 13 %, respectively. Moreover, in a study of ET in vineyards, Leuning et al. (2008) found that the PM model coupled with canopy resistance overestimated ET during the entire growth stage by 29 %.

The estimates for ET produced by the PT–FC and AA models were generally lower than the measured values during the entire season. In addition, the four models also underestimated ET during periods of partial cover (LAI < 2.5 m² m⁻²). The PT–FC and AA models consistently underestimated ET, especially during the late maturity stage. The underestimation probably resulted from the following: (1) nonclassical situations, such as the oasis effect, may occur in the study area. Strong evaporation from the moist ground and plants results in latent heat cooling. However, this upward latent heat flux was opposed by a downward sensible heat flux from the warm air to the cool ground, and thus the latent heat flux was positive while the sensible heat flux is negative. Therefore, the latent heat flux can be greater in magnitude than the solar heating because of the additional energy extracted from the warm air by evaporation (Stull, 1988). (2) The lack of mechanistic representation of rainfall interception in ET models probably led to inaccurate simulation for periods soon after rainy days. Bohn and Vivoni (2016) found that evaporation of canopy interception accounted for 8 % of the annual ET across the North American monsoon region. Comparing the AA and PT–FC models, the former includes forcing data of available radiation, soil water content, and relative humidity, but the PT–FC model only requires available radiation and soil water content and is independent of relative humidity. However, the similar statistical results and similar degrees of MBE scatter indicate that relative humidity has little influence on the AA model simulation. The consistent and consecutive underestimations of ET by the PT–FC and AA models during the late maturity stage show that the model–data disagreement is not caused by regional advection and rainfall interception because atmospheric processes and thermally induced circulation can only occur at certain times and during certain days. Therefore, we think that the consistent underestimation of ET by the PT–FC and AA models results primarily from conceptual and structural inadequacies, energy imbalance, and soil water stress. Although the PM and SW models share a common theoretical basis and the PT–FC model is a simplification of the PM model, these models perform significantly differently. Part of the overestimation of ET by the PM and SW models, caused by coupling with the canopy resistance, may be offset by underestimation caused by energy imbalance and soil water stress. However, underestimation of ET by the PT–FC and AA models cannot be counterbalanced by overestimation during the later maturity stage because the PT–FC and AA models are independent of the canopy resistance. Consequently, the half-hourly patterns of errors in the estimates of ET by the PM and SW models are characterized by symmetry and a low degree of scatter, but the PT–FC and AA models exhibit consistently asymmetrical error patterns. By contrast, other studies showed that the PM model (Kato et al., 2004) and the SW model (Chen et al., 2015) underestimated half-hourly ET. As for the PT–FC and AA models, some studies reported that the PT–JPL (Zhang et al., 2017) and the AA model showed an overall poor performance (Zhang et al., 2017), while other studies have indicated that the AA method performed well for both maize and canola crops (Liu et al., 2012). Therefore, the performance of the four ET models appears to vary not only for different crops and locations but also for different meteorological, physiological, and soil conditions. Moreover, the performance is also related to the stage of crop growth. Note that these conclusions about the ET models evaluation are derived from traditional error metrics rather than those based on BME model selection. It would be desirable to use available data from other study areas or from other crops for BME-based model selection to confirm whether the SW model is the optimal model under other conditions. Overall, combined with the parameter uncertainty analysis described in Sect. 4.1, we conclude that energy imbalance and energy interaction between the canopy and soil surface have a greater impact on the model performance; thus, explicitly treating energy error and incorporating the elements of existing hydrologic theory about energy interaction between the canopy and surface or conceptually correcting the energy interaction are a practicable option for model improvement and application.

5 Conclusions

This study illustrated the application of the Bayesian approach on the statistical analysis and model selection of four widely used ET models. The results showed that the DREAM algorithm successfully reduced the assumed prior uncertainties for most of the parameters in the four models. In the model calibration, the key parameters which had a significant influence on ET simulations were well constrained. The main reasons for the outperforming of SW model were its physically rigorous structure and the extinction coefficient parameter, which is sensitive and has a significant impact on the performance of the model, being well constrained. BME is a consistent model selection for identifying the best fitting to the observed data. Although the squared-residual-based metrics, including R², IA, RMSE, and EF, produced a ranking identical to that of BME, it must be noted that these squared-residual-based metrics do not allow using prior information and do not penalize the model complexity when comparing the models. Therefore, some caution is needed when using these statistical methods to compare different models.

The model–data discrepancies were analyzed to facilitate model improvement after Bayesian model calibration and comparison. The results indicate that the discrepancies arose mainly as a result of energy imbalance caused by anisotropic turbulence, additional energy induced by advection processes, the absence of a mechanistic representation of the physiological response to plant hydrodynamics and the energy interaction between the canopy and surface. Among these causes, energy imbalance and additional energy are related to forcing data errors rather than to an unreasonable model structure. Thus, understanding the process of the physiological response to plant hydrodynamics, and the interaction between the canopy and surface is essential for improving the performance of evapotranspiration models. Overall, the applications of Bayesian calibration, Bayesian model evaluation, and analysis of model–data discrepancies in our study provide a promising framework for reducing uncertainty and improving the performance of ET models. It would be desirable to confirm whether the SW model is the optimal model using data of other crops or other climate regions.

Appendix A: List of symbols and physical characteristics in ET models

A	Available energy for the whole canopy (W m−2)
As	Available energy for the soil surface (W m−2)
Rn	Net radiation fluxes into the canopy (W m−2)
Rns	Net radiation flux into the substrate (W m−2)
G	Soil heat flux (W m−2)
λET	Sum of the latent heat flux from the crop (λT) and soil (λE; W m−2)
ETc	Canopy transpiration (W m−2)
ETs	Soil evaporation (W m−2)
Cc	Canopy resistance coefficient (dimensionless)
Cs	Soil surface resistance coefficient (dimensionless)
LAI	Leaf area index
Q50	Visible radiation flux when stomatal conductance is half its maximum value (W m−2)
D50	Vapor pressure deficit at which stomatal conductance is half its maximum value (kPa)
Da	Vapor pressure deficit at the reference height ($D_{\mathrm{a}}=e_{\mathrm{s}}-e_{\mathrm{a}}$; kPa)
Qh	Flux density of visible radiation at the top of the canopy (W m−2)
Kq	Extinction coefficient
Ka	Extinction coefficient
f	Fraction of evaporation soil and total evaporation
λ	Latent heat of water evaporation (MJ kg−1)
Δ	Slope of the saturated vapor pressure curve (Pa K−1)
γ	Psychrometric constant (kPa K−1)
ρ	Density of air (kg m−3)
k	Karman constant (0.41)
es	Saturated vapor pressure (kPa)
ea	Actual vapor pressure (kPa)
q*	Saturation-specific humidity at air temperature (kg kg−1)
q	Specific humidity of the atmosphere (kg kg−1)
b1	Empirical constant (s m−1)
b2	Empirical constant (s m−1)
β1	Empirical constant
β2	Empirical constant
θ	Soil water content (m3 m−3)
θa	Critical water content at which plant stress starts (m3 m−3)
θb	Water content at the wilting point (m3 m−3)
θr	Residual soil water content (m3 m−3)
θs	Saturated water content (m3 m−3)
Θ	Relative water saturation
d	Zero-plane displacement height (m)
zm	Height of the wind speed and humidity measurements (3 m)
z0 m	Roughness length governing the transfer of momentum (m)
z0 v	Roughness length governing the transfer of water vapor (m)
h	Canopy height (m)
uz	Wind speed at height zm (m s−1)
ga	Aerodynamic conductance (m s−1)
gs	Surface conductance (m s−1)
gmax	Maximum stomatal conductance of leaves at the top of the canopy (m s−1)
$g_{\mathrm{s}}^{\mathrm{c}}$	Canopy conductance (m s−1)
ra	Aerodynamic resistance (s m−1)
$r_{\mathrm{a}}^{\mathrm{a}}$	Aerodynamic resistance between canopy source height and a reference level (s m−1)
$r_{\mathrm{a}}^{\mathrm{s}}$	Aerodynamic resistance between the substrate and the canopy source height (s m−1)
$r_{\mathrm{a}}^{\mathrm{c}}$	Bulk boundary layer resistance of the vegetation element in the canopy (s m−1)
$r_{\mathrm{s}}^{\mathrm{s}}$	Surface resistance of the canopy (s m−1);
$r_{\mathrm{s}}^{\mathrm{c}}$	Bulk stomatal resistance of the canopy (s m−1)

Appendix B: Bayesian inference and the DREAM algorithm

The posterior probability distribution of the parameter is calculated by Bayes' theorem:

$$\begin{array}{}\text{(B1)}& \mathit{\pi }\left(\mathit{\theta }\right|D,M)={\displaystyle \frac{\mathit{\pi }\left(\mathit{\theta }\right|M\left)p\right(D|\mathit{\theta },M)}{p\left(D\right|M)}},\end{array}$$

where π(θ|M) represents the prior density of θ under model M, p(D|θ, M) is the joint likelihood of model M and its parameters θ, and the marginal likelihood, or Bayesian model evidence (BME), is

$$\begin{array}{}\text{(B2)}& p\left(D\mathrm{|}M\right)=\int p(D\mathrm{|}\mathit{\theta },M)p\left(\mathit{\theta }\mathrm{|}Mt\right)d\mathit{\theta }).\end{array}$$

The likelihood function, p(D|θ, M), used for parameter estimation, is specified according to the distributions of observation errors. Error e(t) in each observation D(t) at time t is expressed by

$$\begin{array}{}\text{(B3)}& e\left(t\right)=D\left(tt\right)-f\left(t\right).\end{array}$$

Assuming that e(t) follows a Gaussian distribution with a zero mean, the likelihood function can be expressed as

$$\begin{array}{}\text{(B4)}& p\left(D\right|\mathit{\theta })=\prod _{t=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }}}{e}^{-\frac{\left[e\right(t){]}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

where n is the number of observations and σ represents the error variances.

In this study, we used the DREAM algorithm (Vrugt et al., 2008, 2009) to explore the ET models' parameter space and to estimate BME. The DREAM sampling scheme is an adaptation of the global optimization algorithm of a shuffled complex evolution Metropolis (SCEM-UA). This algorithm was described in more detail in Vrugt et al. (2008, 2009).