Different Vegetation Information Inputs Signiﬁcantly Affect the Evapotranspiration Simulations of the PT-JPL Model

: Evapotranspiration (ET) is an essential part of the global water cycle, and accurate quan-tiﬁcation of ET is of great signiﬁcance for hydrological research and practice. The Priestley-Taylor Jet Propulsion Laboratory (PT-JPL) model is a commonly used remotely sensed (RS) ET model. The original PT-JPL model includes multiple vegetation variables but only requires the Normalized Difference Vegetation Index (NDVI) as the vegetation input. Other vegetation inputs (e.g., Leaf Area Index (LAI) and Fraction of Absorbed Photosynthetically Active Radiation (FAPAR)) are estimated by the NDVI-based empirical methods. Here we investigate whether introducing more RS vegetation variables beyond NDVI can improve the PT-JPL model’s performance. We combine the vegetation variables derived from RS and empirical methods into four vegetation input schemes for the PT-JPL model. The model performance under four schemes is evaluated at the site scale with the eddy covariance (EC)-based ET measurements and at the basin scale with the water balance-based ET estimates. The results show that the vegetation variables derived by RS and empirical methods are quite different. The ecophysiological constraints of the PT-JPL model constructed by the former are more reasonable in spatial distribution than those constructed by the latter. However, as vegetation input of the PT-JPL model, the scheme derived from empirical methods performs best among the four schemes. In other words, introducing more remotely sensed vegetation variables beyond NDVI into the PT-JPL model degrades the model performance to varying degrees. One possible reason for this is the unrealistic ET partitioning. It is necessary to re-parameterize the biophysical constraints of the PT-JPL model to ensure that the model obtains reasonable internal process simulations, that is, “getting the right results for right reasons.”


Introduction
Evapotranspiration (ET) is an essential part of the global water cycle, which is the second-largest water flux behind precipitation. Approximately 60% of global land precipitation returns to the atmosphere through ET [1]. ET even accounts for more than 90% of annual precipitation in some arid regions [2]. In this sense, the ratio of ET to precipitation directly determines the regional water availability on the mean annual scale. Therefore, accurately quantifying ET is of great importance for hydrological practices such as water resources management, crop yield prediction, and drought forecasting [3][4][5][6]. However, ET is a complex natural process, regulated by multiple interacting factors such as climate, topography, soil, and vegetation, and often has great variability in space and time [7,8]. Accurately quantifying the spatial-temporal change in ET at large spatial scales is challenging. Unlike other water balance components (e.g., precipitation and runoff), ET cannot be measured directly due to its invisibility. Existing indirect measurement methods (e.g., the the model performance under the four schemes is systematically evaluated at the site and basin scales. Note that empirical methods are not entirely free of RS vegetation information, and they require RS NDVI to estimate other vegetation variables (LAI, FAPAR, and FIPAR).

Description of the PT-JPL Model
The Priestley-Taylor equation has been widely used for calculating potential ET under unstressed conditions [25]. The general form of the Priestley-Taylor equation is expressed as: where α is the Priestley-Taylor coefficient over wet surfaces, ∆ is the slope of saturated vapor pressure curve (kPa • C −1 ), γ is the psychrometric constant (~0.066 kPa • C −1 ), R n and G are net radiation and ground heat flux (W m −2 ), respectively. G is often ignored in some ET algorithms since it is small compared to R n , especially when the land surface is fully covered by vegetation or the calculation time steps are daily or longer [19]. For the underlying surface of bare soil or sparse vegetation, however, G may play an important role in ET estimation. Here, an empirical formula based on vegetation coverage is used for G estimation [32]. Based on the Priestley-Taylor equation, Fisher et al. [26] developed a RS-based ET algorithm that uses multiple biophysical constraints to reduce the potential ET to actual ET. The model assumes that NDVI can well capture the chlorophyll changes of a canopy because these vegetation indices based on NIR differences (e.g., NDVI) are sensitive to chlorophyll changes. Based on this assumption, other vegetation variables (e.g., LAI, FAPAR, and FIPAR) in the model can be calculated by NDVI. It balances the complexity and accuracy of the model well and has been widely used for regional or global-scale ET simulations [30,33]. The PT-JPL model runs at the daily timescale, and it partitions ET into canopy transpiration (ET c ), soil evaporation (ET s ), and interception evaporation (ET i ). Detailed formulas of the PT-JPL model are shown in Table 1.
Unlike the original PT-JPL model, here we use the daily average (rather than midday) relative humidity (RH) and vapor pressure deficit (VPD) to calculate the soil moisture constraint (f SM ) because of the better model performance obtained from the daily average values (see Table S1). The optimum plant growth temperature (T opt ) used in the function plant temperature constraint (see Table 1) is set to 25 • C following [2], which is different from the original setting of the PT-JPL model. T opt = 25 • C can provide better ET simulations than the original parameterization of T opt [26] in the PT-JPL model at both site and basin scales (see Table S2). The LAI used in the original PT-JPL model is the total LAI that includes the green and non-green leaf area of a canopy per unit ground area [26], while RS LAI is essentially the green LAI. Strictly speaking, RS LAI does not meet the input requirements of the original PT-JPL model. However, it is difficult to obtain reliable total LAI globally, either by RS or empirical methods [34]. In practice, RS LAI data are widely used for allocating energy availability to the canopy and soil surface [35].
Priestley and Taylor [25] suggested an average of 1.26 for α. In practical applications, the coefficient α was found to vary from 0.72 to 1.57, depending on surface vegetation and microclimatic conditions [36]. Here, we insist on using the original coefficient of the Priestley-Taylor model since it has commanded substantial experimental support, especially in humid regions [37,38]. Table 1. Summary of parameters and equations for the PT-JPL model. R n is the net radiation, R nc is net radiation to the canopy (R n − R ns ), R ns is net radiation to the soil (R n exp(−k Rn LAI)), RH is relative humidity (%), T max is maximum air temperature, T opt is optimum plant growth temperature, and the default value of T opt is 25 • C [2], VPD is saturation vapor pressure deficit, FAPAR max is the maximum value of FAPAR in a year, k Rn is extinction coefficient (k Rn = 0.60) [39], k PAR = 0.5 [40], [41], m 2 = 1.0, b 2 = −0.05 [26], SAVI is soil adjusted vegetation index (NDVI × 0.45 + 0.132) [28].

Variables Description Equation References
ET

The Setting of Vegetation Input Schemes for the PT-JPL Model
The PT-JPL model contains three vegetation variables: LAI, FAPAR, and FIPAR. LAI and FAPAR can be obtained from multiple RS data sources [43][44][45][46][47], while global RS FIPAR product is not available. Here we combine LAI and FAPAR estimates from empirical methods and RS products into four vegetation input schemes for the PT-JPL model (  The specific empirical equations used for LAI and FAPAR estimation can be found in Table 1. RS-based LAI and FAPAR refer to the GLASS LAI and FAPAR products [44]. The GLASS LAI product is generated from general regression neural networks, and it showed high accuracy during comparison with similar products and validation with in situ observations [48]. The GLASS FAPAR product is derived from the GLASS LAI product through a series of algorithms [45]. The two products have a spatial resolution of 0.05 • and a time span from 1981 to 2018. The NDVI dataset is obtained from the Global Inventory Monitoring and Modeling System (GIMMS) project [49]. This dataset is generated from multiple satellite sensors and considers multiple factors affecting data quality, such as calibration loss, orbital drift, and volcanic eruptions [50].

The Forcing Data and Evaluation Data for the PT-JPL Model
The forcing data of the PT-JPL model includes meteorological data and albedo data in addition to vegetation input data. Meteorological data includes daily temperature, relative humidity, wind speed, and sunshine duration. We use a professional meteorological interpolation software, namely the Anusplin [51], to interpolate the gauge-based data to gridded data. The Anusplin's algorithm can consider the effect of terrain on interpolated variables. The albedo and sunshine duration data are used to calculate the R n , following the method of Allen et al. [19]. We validate R n estimates with flux site observations. The results indicate that there is a good agreement between our estimated R n and observed R n (see Figure S1). The resolution of original NDVI data is 1/12 degree, and the data is resampled to 0.05 • × 0.05 • with the nearest neighbor resampling approach [52]. The 8-day or half-month RS vegetation data are linearly interpolated to the daily data to match the time step at which the model runs. Natural runoff data from 1981 to 2000 are provided by the Hydrological Bureau of the Ministry of Water Resources in China. More information on the forcing and validation data of the PT-JPL model can be found in Table 2. The performance of the PT-JPL model under four vegetation input schemes is evaluated at the site and basin scales. At the site scale, the EC-based ET measurements from 10 ChinaFlux sites are used as the benchmark data to evaluate the model performance (see Table 3). These flux sites cover different types of ecosystems in China, including three forest sites, four grassland sites, one shrubland site, and two cropland sites ( Figure 2 and Table 3). The main crop types at the two cropland sites are winter wheat and summer maize. The observation data from rainy days (daily precipitation greater than 0.5 mm) are excluded from model evaluation given the potentially large observation errors on rainy days. At the basin scale, the water balance-based ET estimates are used to validate the PT-JPL model's performance. The water balance method assumes that the total storage change (∆S) can be ignored on the mean annual scale. Then, mean annual ET can be estimated as the difference between precipitation (P) and runoff (R): ET = P-R. This method has been widely used for the basin-scale ET evaluation [53,54]. A total of 286 basins are used here (Figure 2), and these basins cover diverse climatic and landscape conditions, with the aridity index ranging from 0.47 to 2.73 [55]. The basin-scale model evaluation is only performed for the period 1982-2000 because runoff data outside this period are not available. In addition, the Global Evaporation Amsterdam Model (GLEAM) ET product (version 3.5) [56,57] is used to evaluate the reliability of ET component partitioning for the PT-JPL model.

Model Performance Assessment
Three model performance indicators are used: coefficient of determination (R 2 ), percent relative bias (Bias, %), and Kling-Gupta efficiency (KGE) [58]. The value of R 2 ranges from 0 to 1, and R 2 = 1 is the optimal value. Bias measures the degree to which simulations are overestimated (positive value) or underestimated (negative value) relative to observations, and it ranges from −∞ to +∞, with an optimal value of zero. KGE is a comprehensive indicator to measure the consistency between simulations and observations. It integrates three independent indicators (bias, variance, and correlation) into one function [59]. The value of KGE is in the range between −∞ and 1.0, with an optimal value of 1.0.

Comparison of Vegetation Variables Estimated by Empirical and Remote Sensing Methods
We compare LAI and FAPAR estimated by empirical methods and RS (GLASS) products. Figure 3 shows the spatial patterns of mean annual LAI (FAPAR) estimates from the two methods and their difference across China. The empirically based LAI shows a similar spatial pattern to RS LAI, with a decreasing trend from southeast to northwest China (Figure 3a,b). However, the two methods have significant differences in LAI estimates, particularly in well-vegetated areas (Figure 3c). On a national scale, mean annual LAI estimates from the empirical method and RS product are 0.81 and 0.92 m 2 /m 2 , respectively. The empirically based LAI has lower values than RS LAI in most areas. The FAPAR estimated by the two methods shows a similar spatial pattern to the LAI (Figure 3d,e), and empirically based FAPAR in most areas is higher than RS FAPAR (Figure 3f). We also compare the interannual changes of national-scale LAI (Figure 4a) and FA-PAR (Figure 4b) estimated by the two methods from 1982 to 2015. The empirically based LAI is apparently less than RS LAI, but the empirically based FAPAR is always larger than RS FAPAR. The annual trends of the two vegetation variables from RS products ((LAI: 2.9 × 10 −3 , FAPAR: 7 × 10 −4 ) are larger than those from empirical methods (LAI: 2.7 × 10 −3 , FAPAR: 4 × 10 −4 ).

The Difference in ET Estimates under Four Vegetation Input Schemes
We compare mean annual ET estimates of the PT-JPL model under four vegetation input schemes ( Figure 5). Mean annual ET estimates from the four schemes show a consistent spatial pattern: a decreasing trend from southeast to northwest China. However, there are pronounced differences in national average ET estimates under the four schemes. Mean annual ET estimates from schemes I and III are close (440.6 mm/year versus 445.7 mm/year), and their values are significantly larger than those obtained from schemes II (394.4 mm/year) and IV (416.6 mm/year). Spatial differences in mean annual ET estimates under the four schemes are mainly reflected in areas with higher ET values, such as southeast China.
We also compare the interannual variability and seasonal cycle of national average ET estimates under the four schemes from 1982 to 2015 ( Figure 6). National-scale annual ET estimates show similar interannual variability but different trends under the four schemes. The trends in annual ET under the four schemes (from I to IV) are 0.05, 0.02, 0.08, and −0.14 mm/year, respectively. The four schemes produce a reverse V-shape intra-annual variation of ET, and the highest values occur in July (Figure 6b). Mean monthly ET estimates from schemes I and III are close to each other, and their values are apparently higher than those from schemes II and IV.

ET Assessment at the Site and Basin Scales
We use the EC-based ET observations from ten flux sites to evaluate the performance of the PT-JPL model under four vegetation input schemes (Figure 7). In terms of the KGE score, scheme I performs best and slightly outperforms scheme III, followed by schemes   We also evaluate mean annual ET estimates under the four schemes using the water balance-based ET estimates ( Figure 8). As shown in Figure 8, scheme I is generally superior to other schemes, which performs best in terms of R 2 and KGE. The R 2 values obtained by the four schemes (from I to IV) are 0.72, 0.65, 0.71, and 0.65, respectively. The KGE values under the four schemes are 0.81, 0.80, 0.80, and 0.80, respectively. In addition, scheme III is slightly inferior to scheme I, but it generally outperforms the other two schemes.

Why Does Introducing More RS Vegetation Variables beyond NDVI Degrade the Performance of the PT-JPL Model?
Some RS ET models require multiple vegetation variables as inputs to allocate the energy and construct environmental constraints for ET [2,26,30,60,61]. There are usually two strategies to obtain these vegetation variables. The first strategy is to use the existing RS data as much as possible [2,62,63]. The second strategy is to use a single RS vegetation variable (e.g., NDVI) in combination with empirical equation(s) to estimate other vegetation variables (e.g., LAI and FAPAR). Both strategies have their advantages and shortcomings. The first strategy can make full use of multi-source RS vegetation information, but the uncertainties of RS vegetation variables from different sources may also be propagated into ET simulations [35]. The second strategy simplifies the model's vegetation inputs but may introduce the uncertainty of empirical equation(s) into ET simulations [64]. This study evaluated the performance of the PT-JPL model under the above two strategies and their combinations. We first compared the differences in two vegetation inputs (LAI and FAPAR) estimated by empirical methods and RS products. The results indicated that there are considerable differences in the LAI and FAPAR estimated by the two methods (Figures 3 and 4). Considering that the RS LAI and FAPAR products are generated based on complex algorithms and have been extensively verified globally [48,65], we have reason to believe that the two vegetation variables estimated by empirical methods have larger uncertainties than RS products. Figure 9 confirms our hypotheses that the RS-derived green canopy fraction (f g ) and plant moisture constraint (f M ) (see Table 1) are more reasonable on spatial patterns than those estimated by empirical equations. The RS-derived f g and f M show similar spatial patterns to LAI and aridity index (Figure 9a,c), respectively, which match the physical meaning represented by the two state variables. In contrast, the empirically based f g and f M are much higher than those derived by RS products, and their spatial patterns are unreasonable in terms of the physical meaning they represent. However, why do better f g and f M estimates not yield better ET simulations? To answer this question, we further examined the proportion of ET components (ET i , ET s , and ET c ) across China under the four vegetation input schemes. An extensively validated global ET product, namely the Global Land Evaporation Amsterdam Model (GLEAM v3.5) [56,57], was used as benchmark data to evaluate the reliability of ET component partitioning for the PT-JPL model. Under the four vegetation input schemes, the PT-JPL model yields a much larger contribution from ET i (13-15%) than the GLEAM (6%) product in China (Figure 10), and the contribution of ET i to ET even exceeds 40% in parts of southern China (see Figure S2). Some studies have also reported the overestimation of ET i by the PT-JPL model [66,67]. The reason for the overestimation of ET i in the PT-JPL model is most likely the unreasonable parameterization of relative surface wetness (f wet , it is also called the faction of wet canopy). The PT-JPL model uses relative humidity (RH) to parameterize f wet (f wet = RH 4 ) if RH > 70%. This implies that ET i occurs as long as RH > 70%, even in the absence of precipitation. This assumption is problematic since RH and precipitation are often decoupled, especially in humid regions of China, where there are many days during a year when RH > 70% but no precipitation occurs.  In addition to the overestimation of ET i , the proportion of transpiration (ET c ) to ET from the PT-JPL model is likely to be underestimated (see Figure 10). Three potential reasons contribute to the low proportion of ET c to ET across China. First, the energy allocated to the canopy layer was underestimated due to the underestimation of LAI (see Figure 4a). The LAI defined in the original PT-JPL model is the total LAI, while GLASS LAI is essentially the green LAI. In theory, the GLASS LAI should be less than or equal to total LAI if a reasonable LAI parameterization is used. However, the NDVI-derived LAI is significantly smaller than the RS LAI (see Figure 4a). Second, ET c may be over-constrained in the PT-JPL model: three biophysical constraints (f g , f T , and f M ) are used in the model. A low value for any of these constraints would make ET c much smaller than its potential value (the value under unstressed conditions). By contrast, the GLEAM uses the same method of potential evaporation as the PT-JPL but only uses a soil water constraint to reduce potential evaporation to actual evaporation [56,57]. Third, the proportion of ET s to ET may be overestimated in the PT-JPL model to offset the underestimation of ET c , thereby ensuring the rationality of ET simulations. The overestimation of ET s is achieved by the underestimation of LAI and the overestimation of soil moisture constraint (f SM ). The model calculates the f SM with daily relative humidity and VPD. This constraint may well reflect the soil moisture status in humid areas but not in arid areas [68]. For example, the desert and Gobi areas of northwest China have extremely dry climates, where the mean annual aridity index exceeds 5.0 (see Figure S3). The real f SM in these areas should be close to zero, but f SM estimates in these areas range from 0.3 to 0.6. Overall, the PT-JPL model provides an efficient tool to simulate ET, but the rationality of ET partitioning needs to be revisited. The unrealistic ET portioning is largely attributable to inadequate model constraint schemes, that is, "getting the right results for the wrong reasons." This explains why introducing more remote sensing variables other than NDVI degrades the model's performance to varying degrees.
From the perspective of the model application, it is redundant to introduce more RS vegetation variables beyond NDVI, if we only pursue the rationality of ET simulations. However, model users sometimes need to divide ET components reasonably. In this case, it is necessary to re-parameterize the biophysical constraints of the PT-JPL model to achieve a reasonable ET partitioning. Specific model improvement suggestions are summarized as follows: (i) increasing the exponent of the f wet function (f wet = RH 4 ) to reduce the fraction of ET i to ET, (ii) introducing RS LAI and FAPAR data, and (iii) re-parameterizing the biophysical constraints based on multi-source ET observations.

The Model Sensitivity to Vegetation Variables
The differences in ET simulations between different vegetation input schemes essentially reflect the model sensitivity to vegetation input variables. For example, the LAI estimates from schemes I and III differ greatly (Figure 4a), but little difference in ET estimates was found between the two schemes. This means that the model is not sensitive to the input of LAI. To fully demonstrate the model sensitivity to vegetation input variables (LAI, FAPAR, and FIPAR), we exerted four perturbations (−10%, −5%, 5%, 10%) to each vegetation variable one by one and then calculated the relative difference in ET simulations (∆ET, %) from baseline simulations (the results of the scheme I). The larger the ∆ET, the more sensitive the model is to the change in vegetation variable. The results show that the PT-JPL model is most sensitive to the change in FAPAR, followed by FIPAR and LAI ( Figure 11). This result explains why schemes I and III behave similarly, while considerable differences exist between the scheme I and schemes II and IV.

Potential Uncertainty Sources
Many factors can cause uncertainty in the evaluation results. First, the EC-based ET measurements have two inherent deficiencies: the non-closure of energy balance components and the scale mismatch during the "site-to-pixel" model evaluation [4]. The issue of energy balance non-closure has been observed on almost all flux stations globally [69][70][71]. The non-closure degree of surface energy components in some European sites is generally between 10-30%, as reported by Mauder et al. [72]. Numerous factors affect the closure degree of the energy balance components, including the sampling errors, instrument bias, neglected energy sinks, and horizontal and/or vertical advection of heat and water vapor [69]. In addition, the spatial representation of an EC flux tower is about several hundred square meters, depending on the measured height above the canopy layer and wind speeds [9]. However, the grid area of the model calculation is about 18-29 km 2 (depending on the altitude). The scale mismatch also has the capacity to skew evaluation results, given the nature of considerable spatial heterogeneity in ET [4]. The uncertainty in the water balance-based ET estimates may also affect the evaluation results. The water balance method assumes that the changes in water storage (S) can be ignored on the mean annual timescale, and then the basin-average ET is estimated as the difference between precipitation and runoff. However, the assumption of S = 0 may not hold, especially in small basins [53,73]. Small basins often have frequent water exchange with neighboring basins, and ignoring this part of the water may result in a large bias in ET estimation from the water balance method [74].
The uncertainty of the forcing data is an important error source for ET simulations [75]. In this study, the meteorological variables from 824 meteorological stations were interpolated into gridded data to drive the PT-JPL model. However, the spatial distribution of these stations is sparse and uneven, especially in the western regions (see Figure 2). This may affect the accuracy of the interpolated meteorological data. In addition, R n is the key input for the PT-JPL model, and its error may greatly affect the output of the PT-JPL model [26,75]. However, R n is not a routine observation item of meteorological stations in China. Here we use the sunshine duration, RS albedo, and FAO-56 PM method to calculate the daily R n [19]. R n is calculated as the difference between incoming and outgoing radiation of both short and long wavelengths. This method itself has some flaws. It uses the air temperature instead of the land surface temperature to calculate the outgoing longwave radiation since the land surface temperature data before 2000 are not available. This may result in bias in R n estimates [31]. Nevertheless, the validation results at flux sites indicated that the R n estimates from the FAO-56 PM method are generally reliable.
The PT-JPL model under the four vegetation input schemes systematically underestimates the ET at the two cropland sites (see Figure 7). This may be due to the effect of agricultural irrigation, which can significantly increase the regional ET rate. However, the PT-JPL model does not take into account the effect of irrigation on ET, thus underestimating ET estimates at these sites. Previous studies have also reported the underestimation of ET at these cropland sites [76,77]. In addition, a single RS vegetation data source (i.e., GLASS) was used here, and these data inevitably have varying degrees of uncertainty. Replacing GLASS with other RS data sources may have the effect of skewing the evaluation results.

Conclusions
The major purpose of this study is to investigate whether introducing more vegetation variables beyond NDVI into the PT-JPL model can improve the model's performance. The vegetation variables estimated by empirical methods and RS (GLASS) products were combined into four different vegetation input schemes for the PT-JPL model. We then evaluated the performance of the PT-JPL model under four schemes at the site and basin scales. The main conclusions are as follows.
(1) Introducing more RS vegetation information beyond NDVI degrades the accuracy of ET simulations for the PT-JPL model to varying degrees. (2) A possible reason for this is the misinterpretation of ET components caused by unreasonable parameterization schemes of biophysical constraints. This study highlights the importance of vegetation variable inputs from different sources to ET simulations. 21 May 2022). The GLEAM v3.5 datasets are available at https://www.gleam.eu/ (accessed on 16 February 2022).

Conflicts of Interest:
The authors declare no conflict of interest.