The Modiﬁed Soil Moisture Constraint Scheme Signiﬁcantly Enhances the Evapotranspiration Simulation Accuracy of the MOD16 Model

: Remotely sensed (RS) evapotranspiration (ET) models can make full use of the land surface information retrieved using remote sensing and are therefore widely used in large-scale ET estimates. The MODIS Global Evapotranspiration model (MOD16) is one of the most commonly used remote sensing ET models. MOD16 parameterizes the moisture constraints on soil evaporation (Es) using atmospheric vapor pressure deﬁcit (VPD) and relative humidity (RH). This moisture constraint algorithm has been criticized by many studies due to the weak correlation between soil moisture and VPD or RH over short timescales (e.g., hourly and daily). In this study, we introduce a modiﬁed moisture constraint algorithm of ET, based on the ratio of antecedent accumulated precipitation to soil equilibrium evaporation, in order to improve the ET simulation capabilities of the MOD16 model. The original and modiﬁed MOD16 models are evaluated at 14 ChinaFlux sites and 286 basins in China, using eddy covariance measurements and water-balance-based ET estimates. The results show that the modiﬁed MOD16 model outperforms the original MOD16 model at both the site and basin scales. Compared with the original model, the modiﬁed model increases the value of KGE by an average of 0.17 at the ﬂux site scale and by 0.01 at the basin scale. Using soil moisture measurements from ﬂux sites as a reference, we further found that the modiﬁed MOD16 model also has a better soil moisture simulation capacity than the original model. This study highlights the importance of reliable soil moisture constraints in remotely sensed ET models.


Introduction
Evapotranspiration (ET) is an important component of the global water cycle [1,2], and the largest water flux beyond precipitation in terrestrial water budgets [3]. About 60% of annual land precipitation is returned to the atmosphere via ET, with the percentage increasing to as much as 90% in arid and semiarid regions [4]. Moreover, ET is a vital energy flux determined by the partition of available energy into sensible and latent heat fluxes [5,6]. Therefore, it is essential to precisely quantify ET estimation in order to comprehensively understand regional water cycles and to better manage water resources for agriculture, industry, and other hydrological practices.
Traditional ET estimation methods mainly rely on in situ measurements, such as eddy covariance (EC) and the lysimeter and Bowen ratios [7]. However, these techniques may not adequately meet the demand for ET estimates at regional and large scales due to the variability of ET across different types of land cover [8]. In addition, the instruments required for ET measurements are often located in only several sparsely distributed locations, limiting their spatial representativeness and causing high maintenance costs [9]. In comparison to conventional methods, remote sensing (RS) techniques are generally regarded as effective and reliable tools for monitoring regional ET [10]. The RS methods can be used to where s is the slope of the curve relating saturated water vapor pressure to temperature (kPa • C −1 ); A c and A s (MJ m −2 day −1 ) are the available energy allocated to the canopy and soil surfaces, respectively; ρ is the air density (kg m −3 ); C p is the specific heat capacity of air (MJ kg −1 • C −1 ); F c is the fractional vegetation cover (dimensionless); VPD is the atmospheric vapor-pressure deficit (kPa); F wet is the relative surface wetness (dimensionless); λ is the latent heat of evaporation (MJ kg −1 ); P a is the atmospheric pressure (kPa); ε is the ratio of the molecular weight of water to dry air (i.e., 0.622); γ is the psychrometric constant (kPa • C −1 ); F sm is soil water constraint (dimensionless).
In the MOD16 model, F sm follows the soil moisture constraint equation detailed by Fisher et al. [29] and is defined as the function of VPD and RH. It is based on the complementary relationship hypothesis [30] assuming that soil moisture is reflected in the adjacent atmospheric moisture [29]. For ease of presentation, hereafter, F sm is referred as to f VPD : f VPD = (RH/100) (VPD/β) (5) where β is the relative sensitivity to VPD. In the old algorithm by Mu et al. [20], the value of β is 0.1 kPa and has been revised as 0.2 kPa in the modified algorithm by Mu et al. [6]. We set the parameter β at 0.25 kPa in alignment with the latest version by Running et al. [31].
The changes of soil moisture have been proved not only to be impacted by the atmospheric factors (i.e., temperature and humidity), but also the variability of precipitation. According to the method proposed by Morillas et al. [32], we introduce a soil moisture constraint algorithm based on the ratio of previous precipitation and soil equilibrium evaporation to improve the ET simulation of the MOD16 model. Equation (6) is used to calculate f drying during the effective precipitation days (P i > P min = 0.5 mm day −1 ).
where f drying is the modified soil water constraint (dimensionless); P i is the accumulated daily precipitation (mm); f LP is the value for the last effective precipitation day (dimensionless); α (day −1 ) is a parameter controlling the rate of soil drying, higher α values reflecting higher soil drying speed; ∆t is the number of days between this and the current day i. Apart from the variation in the F sm calculation method, other parameters and algorithms of the MOD16 model remain unchanged and are consistent with the settings of Running et al. [31].

Data Sources and Processing
There are two types of data used in this study: meteorological data and surface flux data. Meteorological data were collected from approximately 800 national standard meteorological stations ( Figure 1) and provided by the China Meteorological Administration. Then, these collected data were interpolated into raster data with a resolution of 0.05 • × 0.05 • using the Anusplin method [33]. The Anusplin method is an effective interpolation approach that can accurately model the impact of terrain on interpolated variables [34]. Several studies have demonstrated that this method has a higher interpolation accuracy than other interpolation methods, such as inverse distance weighting (IDW) or ordinary kriging [35,36]. To assess the performance of ET estimates from the modified and original MOD16 methods, evaluations were conducted on both the site and basin scales. The performance of the two methods was evaluated using daily evapotranspiration (ET) observations at 14 ChinaFlux sites and soil moisture observations at 7 out the 14 ChinaFlux sites (see Table  1). It is worth noting that soil moisture observations were only available at 7 of 14 ChinaFlux sites. We use the symbol * to denote the stations at which soil moisture observations were made (see Table 2). These flux sites include multiple land cover types, including forests, grasslands, and croplands (see Figure 1 and Table 2).  Land surface data, including albedo and leaf area index (LAI), were obtained from the GLASS AVHRR product, which provides a resolution of 0.05 degrees and 8-day compositing periods in the simulation projection. To be consistent with the time step of the MOD16 model, these land surface data were then converted into daily data from 1981 to 2018. The albedo data were used to obtain the net radiation based on the methodology proposed by Allen et al. [37]. In this study, the land cover type was classified based on the International Geosphere-Biosphere Programme (IGBP) land cover classification, and land cover data were obtained from the MCD12C1 product. The MOD16 algorithm cannot calculate ET from persistent wetlands and water bodies, as it does not provide relevant parameters for the two land cover types. Here, ET estimates from persistent wetlands and water bodies were obtained from a Chinese open water evaporation dataset (https://osf.io/qd28m/, (accessed on 23 June 2023)). This dataset was developed based on the modified Penman model, considering the impact of heat storage and area on evaporation estimations [38].
To assess the performance of ET estimates from the modified and original MOD16 methods, evaluations were conducted on both the site and basin scales. The performance Sustainability 2023, 15, 12460 5 of 13 of the two methods was evaluated using daily evapotranspiration (ET) observations at 14 ChinaFlux sites and soil moisture observations at 7 out the 14 ChinaFlux sites (see Table 1). It is worth noting that soil moisture observations were only available at 7 of 14 ChinaFlux sites. We use the symbol * to denote the stations at which soil moisture observations were made (see Table 2). These flux sites include multiple land cover types, including forests, grasslands, and croplands (see Figure 1 and Table 2). On the basin scale, the performance of the two soil moisture methods was evaluated using a water balance approach. This approach is based on the assumption that the total storage change (∆S) can be ignored when estimating annual values [39]. Therefore, the mean annual evapotranspiration (ET) can be calculated following the water balance equation: where P is precipitation (mm) and R is runoff. In total, the 286 test basins used in this study ( Figure 1) cover various climate conditions and land cover types. The natural runoff data of the study basins were provided by The soil moisture product generated from the China Land Data Assimilation System (CLDAS) was used to evaluate the performance of the two soil moisture constraint algorithms. The dataset has a spatial resolution of 1/16 degree and a temporal resolution of 1 day [40]. The NCAR/CLM3.5 land surface models generated the CLDAS soil moisture products by integrating the ground observations, satellite observations, and numerical model products as forcing data [41]. The product was evaluated for consistency with in situ observations of soil moisture and is widely used for hydrological applications.

Model Performance Assessment
To assess the performance of the two soil moisture constraint methods, the coefficient of determination (R 2 , Equation (8)), the percent relative bias (bias, Equation (9)), and the Kling-Gupta efficiency (KGE, Equation (10)) were used. R 2 is a metric that can evaluate the goodness of fit between predictions and observations. The range of R 2 lies between 0 and 1, with a value of 1 indicating a perfect match between the predicted and observed values [42]. Bias can measure the deviation between predictions and observations, with positive (negative) values indicating overestimation (underestimation) [43]. The Kling-Gupta efficiency (KGE) [44] is a comprehensive indicator of model performance which considers the correlation, bias, and variability between simulations and observations [45].
Negative KGE values indicate a bad model performance and positive KGE values indicate a good model performance, with an optimal value of 1.0.
where x i and y i are the observation and simulation data of grid i, respectively; r is the liner correlation between the observation and simulation data; σ obs and σ sim are the standard deviation in observations and simulations, respectively; µ obs is the mean observation value; and µ sim is the mean simulation value.

The Spatiotemporal Feature of ET Estimates
We compared the spatial patterns of ET estimations and their difference calculated by the modified and original MOD16 methods in China from 1982 to 2018. As is shown in Figure 2a  (p < 0.05) (Figure 3a). However, the trend determined by the modified method is larger than that of the original method, at 1.56 mm/yr 2 for the modified model and 0.43 mm/yr 2 for the original model. On the seasonal scale, the ET estimates from the two methods both show an inverted V-shape variation, with peaks in July at 89.8 mm/month (modified) and 70.6 mm/month (original) (Figure 3b). From April to September, the modified monthly ET is higher than the original one, and the largest difference of 19.2 mm is observed in July. In other months, the modified monthly ET is lower than the original ET.
value; and is the mean simulation value.

The Spatiotemporal Feature of ET Estimates
We compared the spatial patterns of ET estimations and their difference calculated by the modified and original MOD16 methods in China from 1982 to 2018. As is shown in Figure 2a,b, both methods have a similar spatial pattern of mean annual ET estimates, with higher values in the southern and eastern regions of China and lower values in the northern and western regions of China. On a national scale, the difference in the mean annual ET estimates between the modified and original algorithms indicates that the modified algorithm leads to generally higher mean annual ET estimates, especially in arid and semiarid regions (Figure 2c)    methods both show an inverted V-shape variation, with peaks in July at 89.8 mm/month (modified) and 70.6 mm/month (original) (Figure 3b). From April to September, the modified monthly ET is higher than the original one, and the largest difference of 19.2 mm is observed in July. In other months, the modified monthly ET is lower than the original ET.  Table 3 shows the performance of the two methods at the 14 flux sites. At all flux sites, the modified method produced higher R 2 and KGE values than the original method. Relative to the original method, the average R 2 of the modified method increased from 0.52 to 0.68 and the average KGE increased from 0.54 to 0.71. These results demonstrate that the modified method provides more accurate ET estimates than the original algorithm. Regarding different land cover types at the 14 flux sites, the modified method performs particularly well for ET estimates at grassland and cropland sites. For the five grasslands sites, the improvement in R 2 and KGE values between the two methods is larger than 0.20 and 0.10, respectively. At the DX site, the modified method resulted in the greatest increase in R 2 and KGE, from 0.28 to 0.67 and 0.22 to 0.58, respectively. For the six cropland areas, the improvements in the average R 2 and KGE between the two methods are 0.09 and 0.12, respectively. At the HL and DaXing site, the improvement in R 2 and KGE between the two methods was larger than 0.20 for both, indicating that the modified method had a more significant impact at these sites than at other cropland sites. For forest areas, however, there is no significant difference in R 2 and KGE between the two methods compared to other land cover type sites. For instance, the improvement in R 2 and KGE values between the two methods is 0.02 and 0.01, respectively, at the CBS site. At the QYZ site, the improvements of R 2 and KGE are 0.05 and 0.03, respectively. Table 3. Validation of ET estimates from two methods at 14 flux sites. MF: mixed forest; ENF: evergreen needleleaf forests; GRA: grasslands; CRO: croplands; MAP: mean annual precipitation; MAT: mean annual temperature. fVPD represents the performance of the original method and fdrying represents the performance of the modified method. Bold columns indicate better results regarding the three assessment indictors.  Table 3 shows the performance of the two methods at the 14 flux sites. At all flux sites, the modified method produced higher R 2 and KGE values than the original method. Relative to the original method, the average R 2 of the modified method increased from 0.52 to 0.68 and the average KGE increased from 0.54 to 0.71. These results demonstrate that the modified method provides more accurate ET estimates than the original algorithm. Regarding different land cover types at the 14 flux sites, the modified method performs particularly well for ET estimates at grassland and cropland sites. For the five grasslands sites, the improvement in R 2 and KGE values between the two methods is larger than 0.20 and 0.10, respectively. At the DX site, the modified method resulted in the greatest increase in R 2 and KGE, from 0.28 to 0.67 and 0.22 to 0.58, respectively. For the six cropland areas, the improvements in the average R 2 and KGE between the two methods are 0.09 and 0.12, respectively. At the HL and DaXing site, the improvement in R 2 and KGE between the two methods was larger than 0.20 for both, indicating that the modified method had a more significant impact at these sites than at other cropland sites. For forest areas, however, there is no significant difference in R 2 and KGE between the two methods compared to other land cover type sites. For instance, the improvement in R 2 and KGE values between the two  We also evaluated the mean annual ET estimates of the two methods using the water balance approach (Figure 4). The modified method is generally superior to the original method, although the advantages of the modified model over the original model are not as pronounced as on the site scale. Compared to the original method, the KGE value of the modified method increased from 0.72 to 0.73 and the bias value decreased from 12.36% to 8.51%. The red line is clearly nearer to the 1:1 dashed line than the black line. The results further confirm the advantages of the modified method over the original model in ET estimation. We also evaluated the mean annual ET estimates of the two methods using the water balance approach (Figure 4). The modified method is generally superior to the original method, although the advantages of the modified model over the original model are not as pronounced as on the site scale. Compared to the original method, the KGE value of the modified method increased from 0.72 to 0.73 and the bias value decreased from 12.36% to 8.51%. The red line is clearly nearer to the 1:1 dashed line than the black line. The results further confirm the advantages of the modified method over the original model in ET estimation.

Comparison of Different Soil Moisture Constraint Schemes
We compared the spatial pattern of the soil moisture constraint between the two methods and the CLDAS product, as shown in Figure 5. Both the original and modified methods simulated similar spatial patterns of soil moisture consistent with the CLDAS

Comparison of Different Soil Moisture Constraint Schemes
We compared the spatial pattern of the soil moisture constraint between the two methods and the CLDAS product, as shown in Figure 5. Both the original and modified methods simulated similar spatial patterns of soil moisture consistent with the CLDAS product, with a higher soil moisture in vegetated regions and a lower soil moisture in semiarid and arid regions. However, the two methods exhibited different average soil moisture constraints compared to the CLDAS product. Compared to the CLDAS product, the modified method had a smaller difference in the average soil moisture constraint value at 0.04, while the original method had a larger difference at 0.13. These results suggest that the modified method simulates a more accurate soil moisture constraint compared to the original method.   Figure 6 shows the correlation between the soil moisture constraint derived from in situ measurements and the soil moisture constraint estimated by the two model methods. The observations are better correlated with the modified method than the original method. The correlation coefficient between the observations and the modified method ranges from 0.22 to 0.71, with an average value of 0.52. In contrast, the correlation coefficient between the observations and the original method ranges from −0.47 to 0.39, with a smaller average of 0.02. Figure 6. Correlation of the soil moisture constraint (Fsm) derived from measurements and the orig-    Figure 6 shows the correlation between the soil moisture constraint derived from in situ measurements and the soil moisture constraint estimated by the two model methods. The observations are better correlated with the modified method than the original method. The correlation coefficient between the observations and the modified method ranges from 0.22 to 0.71, with an average value of 0.52. In contrast, the correlation coefficient between the observations and the original method ranges from −0.47 to 0.39, with a smaller average of 0.02.

The Reason for the Improved Performance of the Modified Method over the MOD16 Model
The results above indicate that the modified method provides more accurate ET estimates and is better correlated with the observations of soil moisture than the original

The Reason for the Improved Performance of the Modified Method over the MOD16 Model
The results above indicate that the modified method provides more accurate ET estimates and is better correlated with the observations of soil moisture than the original method. To investigate the reasons for different ET estimates from the two methods, we compared the interannual variability of the soil moisture constraint between the two methods and observations at three different land cover sites (Figure 7). The interannual variability of the soil moisture constraint at three different land cover sites reveals that precipitation consistently varies with the observed soil moisture constraint, indicating a high correlation between precipitation and soil moisture. The interannual variability of the modified soil moisture constraint is consistent with that of observations (Figure 7b). However, there are large deviations in the interannual variability of the soil moisture constraint between the original method and observations (Figure 7a). These findings suggest that the modified method is superior to the original method because it accounts for the effects in the variability on precipitation to the soil moisture. variability of the soil moisture constraint at three different land cover sites reveals that precipitation consistently varies with the observed soil moisture constraint, indicating a high correlation between precipitation and soil moisture. The interannual variability of the modified soil moisture constraint is consistent with that of observations ( Figure 7b). However, there are large deviations in the interannual variability of the soil moisture constraint between the original method and observations (Figure 7a). These findings suggest that the modified method is superior to the original method because it accounts for the effects in the variability on precipitation to the soil moisture.

Uncertainties
This study introduced a modified soil moisture constraint algorithm to the MOD16 model and then evaluated the performance of the modified MOD16 method on ET estimations. The results show that the modified MOD16 method outperforms the original MOD16 method in ET estimations, as evidenced by higher R 2 and KGE values. However, Figure 7. Comparison of the interannual variability of the soil moisture constraint (F sm ) between the two methods and observations at three different land cover sites. Red and blue shadows represent the standard deviation from the original method (a,c,e) and the modified method (b,d,f), respectively. Pr: precipitation (mm/day).

Uncertainties
This study introduced a modified soil moisture constraint algorithm to the MOD16 model and then evaluated the performance of the modified MOD16 method on ET estimations. The results show that the modified MOD16 method outperforms the original MOD16 method in ET estimations, as evidenced by higher R 2 and KGE values. However, many uncertainties still exist in the original and the modified MOD16 methods. The biases in the two soil moisture constraint methods arise from several major causes.
First, the MOD16 model forcing data are a major source of uncertainty. Meteorological forcing data collected from 824 national meteorological stations are used in this study. However, the distribution of meteorological stations across China is uneven, with fewer stations located in the western and border regions compared to other regions (Figure 1). The uneven distribution of stations may potentially affect the accuracy of the interpolation of meteorological variables.
Uncertainty may arise from the inaccuracy of ET observations via eddy covariance (EC) measurements. There are two main sources of deficiencies in the EC-based ET measurements: the energy balance closure problem and the spatial mismatch between the situ EC observations and ET estimates of the model. The non-closure energy balance problem is a long-standing issue in EC-based ET estimates. Many studies have indicated that EC measurements tend to underestimate turbulent heat fluxes, and thus the energy balance is not closed [46,47]. Errors in EC-based ET measurements might affect the performance assessment of the modified MOD16 model. In addition, the spatial mismatch between in situ EC observations and ET estimates of the MOD16 model is another issue. Usually, the pix footprint of EC-based ET sites is constrained to several kilometers around the flux station, depending on the measurement height above the canopy layer and the wind speed [48]. Therefore, EC measurements over heterogeneous land surfaces may not accurately represent the grid scale average ET, causing errors in the validation results [49]. Considering the feature of spatial heterogeneity in ET estimates, the issue of scale mismatch has the potential to skew evaluation results.

Conclusions
To develop the accuracy of ET estimates, we introduced a modified soil moisture constraint algorithm to the MOD16 model. Then, we evaluated the performance of the original and modified method at the site and basin scales. We further compared the accuracy of the two methods in capturing soil moisture variability. The main conclusions are as follows: