Estimation of Evapotranspiration and Its Components across China Based on a Modiﬁed Priestley–Taylor Algorithm Using Monthly Multi-Layer Soil Moisture Data

: Although soil moisture ( SM ) is an important constraint factor of evapotranspiration ( ET ), the majority of the satellite-driven ET models do not include SM observations, especially the SM at different depths, since its spatial and temporal distribution is difﬁcult to obtain. Based on monthly three-layer SM data at a 0.25 ◦ spatial resolution determined from multi-sources, we updated the original Priestley Taylor–Jet Propulsion Laboratory (PT-JPL) algorithm to the Priestley Taylor–Soil Moisture Evapotranspiration (PT-SM ET ) algorithm by incorporating SM control into soil evaporation ( E s ) and canopy transpiration ( T ). Both algorithms were evaluated using 17 eddy covariance towers across different biomes of China. The PT-SM ET model shows increased R 2 , NSE and reduced RMSE , Bias, with more improvements occurring in water-limited regions. SM incorporation into T enhanced ET estimates by increasing R 2 and NSE by 4% and 18%, respectively, and RMSE and Bias were respectively reduced by 34% and 7 mm. Moreover, we applied the two ET algorithms to the whole of China and found larger increases in T and E s in the central, northeastern, and southern regions of China when using the PT-SM algorithm compared with the original algorithm. Additionally, the estimated mean annual ET increased from the northwest to the southeast. The SM constraint resulted in higher transpiration estimate and lower evaporation estimate. E s was greatest in the northwest arid region, interception was a large fraction in some rainforests, and T was dominant in most other regions. Further improvements in the estimation of ET components at high spatial and temporal resolution are likely to lead to a better understanding of the water movement through the soil–plant–atmosphere continuum.


Introduction
Understanding water and heat exchange information of the soil-plant-atmosphere continuum (SPAC) is imperative to better manage more restricted water resources in the future. As one of the most important connecting terms in SPAC [1], evapotranspiration (ET), which consists of evaporation from wet surfaces (E) and transpiration through plants (T), is critical for both the energy budget and water balance in the earth-atmospheric system. Normally, vegetation T is a physiological process that is closely associated with CO 2 assimilation and has a great influence on the gross productivity of terrestrial ecosystems, while evaporation from wet canopy surfaces (E i ) or soil (E s ) is often regarded as a wasteful loss of water and does not directly contribute to the ecosystem productivity [2,3]. Therefore, ET partitioning is of prime importance for the water resources and forests management where α is the PT coefficient of 1.26, ∆ is the slope of the saturated vapor pressure curve (kPa • C −1 ), γ is the psychrometric constant (~0.066 kPa • C −1 ), λ is the latent heat of vaporization (kPa • C −1 ), R n is the net radiation at the surface (MJ m −2 ), and G is the ground heat flux (MJ m −2 ). Subsequently, Fisher et al. (2008) developed a modified PT model (i.e., the PT-JPL model) to scale the PET down to the actual ET using eco-physiological constraints and soil evaporation partitioning, which are driven by atmospheric moisture (vapor pressure deficit and relative humidity, VPD and RH) and vegetation indices (normalized and soil adjusted vegetation indices, NDVIs, and the enhanced vegetation index, EVI) [37]. The PT-JPL algorithm estimate ET at the monthly scale is defined as follows: R nc (5) where f wet is the relative surface wetness (RH 4 ); f sm is a soil moisture constraint (RH min VPD min /β , where VPD min is calculated by T max and RH min ); f C is the green canopy fraction (f APAR /f IPAR ) with f APAR being the fraction of absorbed photosynthetically active radiation (PAR) absorbed by green vegetation cover (m 1 EVI + b 1 ) and f IPAR the fraction of intercepted PAR (m 2 NDVI + b 2 ); f T is a plant temperature constraint (exp(−((T max − T opt )/T opt ) 2 )); f M is a plant moisture constraint (the ratio between f APAR and f APARmax ); T max is the maximum air temperature; and T opt is T max at max (R n T max EVI/VPD). R ns is the net radiation at the soil surface calculated as: R ns = R n exp(−k R n LAI) where k R n = 0.60, and LAI is the leaf area index; R nc is the canopy net radiation (R n − R ns ).

PT-SM Algorithm
In the PT-JPL model, the soil water deficit index, f sm (RH min VPD min /β ) is employed to constrain the rate of evaporation from the soil surface according to the Bouchet's complementary hypothesis that the surface moisture status is linked to the evaporative demand of the atmosphere [37]. In the original PT algorithm, the REW is usually used as an index that normalizes the impact based on the soil moisture and soil properties [41]. Note that E s is driven by moisture within the top few centimeters of the soil profile. Therefore, with the aid of observed surface soil moisture (SM, m 3 m −3 ), the scalar f REW can be calculated from the following equation: where SM 1 is the 0-5 cm soil moisture (m 3 m −3 ); SM 1fc and SM 1wp are the soil field capacity and soil moisture at the wilting point, respectively. Therefore, the new E s equation is conducted by replacing f sm with f REW : In the original PT-JPL formulation, the plant moisture stress is inferred from the deviation from maximum greenness ( f M ), which is calculated by normalizing phenological parameters by the maximum observed value per pixel. Actually, at the monthly scale, the latent responses from vegetation to soil moisture deficits can also affect the vegetation T to some extent, especially for arid and semi-arid regions [23]. To evaluate the influences of plant access on atmospheric demand and of deeper water storage intensification or mitigating vegetation sensitivity on water availability, the transpiration response curve to soil water availability can be drawn (the relation between transpiration stress and soil moisture conforms to the exponential form), and the results indicate that higher canopy heights show less sensitivity to soil water availability [24,42]. Furthermore, the root system developed in different vegetations is different. Generally, the depth of root system is related to plant species, but it is also related to the climate [24]. For example, in the subtropical region, there is plenty of rain, tree species can get enough water from the surface soil, so they should not grow deeper into the soil. Therefore, in this study, vegetation is divided into four categories, forest, shrub, grass, and crop, and a new eco-physiological scalar is formulated as follows: where SRH k is the soil relative humidity at k-layer. Three layers (0-5, 5-20, and 20-100 cm) are used in this study, k = 1, 2, 3 for frost and shrub covers, and k = 1, 2 for grass and crop covers. While the transpiration constraint f TREW is supposed to be related to the type and height of vegetation, and the climate condition, which can be calculated by: where SM kCR is the critical k-layer soil moisture at which soil water availability limits ET; CH scalar = √ CH is a canopy height (CH) scalar set to range from 1 to 5, and the square root is formed from Martens et al. (2017) [40]; VT scalar is a vegetation type (VT) scalar that represents the impact of vegetation sensitivity on soil water availability; and CL scalar is a climate (CL) scalar, e.g., the deep-rooted tree species can gradually change into the shallow-rooted species in the subtropical environment. Both the vegetation type and Remote Sens. 2021, 13, 3118 5 of 24 climate scalars are set to 0-1, which are determined by comparing the modeled ET series to the observations. Thus, the algorithm of vegetation transpiration can be described as follows:

Evaluation
Eddy covariance (EC) observations of latent heat flux (LE) across different land cover types are often used for ET model evaluation. The performances of PT-JPL and PT-SM algorithms for estimating monthly ET were evaluated by comparing the simulated monthly ET series against ET values observed by EC sites. Four statistics, namely the coefficient of determination (R 2 ), the Nash-Sutcliffe Efficiency (NSE), the root mean square error (RMSE), and the amount of system deviation (Bias), were used to evaluate the performances of these two PT models. The four statistics are defined in the forms of with ET obs,i = n ∑ i=1 ET obs,i /n as the average of the observed ET i , where n is the number of time-steps.

Eddy Covariance (EC) Flux Data
Ranging from 3.8 to 53.6 • N and from 73.6 to 135.1 • E (Figure 1), China is the third largest national territorial area in the world with a total area of over 9.73 million km 2 . With such a large area, China is characterized by a complex climate and diverse land cover types and soil types ( Figure 1B,C). Ground-measured LE data of 17 eddy covariance sites from China FLUX, USCCC, WATER, and Asia flux, distributed in different land cover types of China, were collected to validate the PT-JPL and PT-SM algorithms ( Figure 1A and Table 1). The land cover types of the flux towers include croplands (Crop; four sites), deciduous needle leaf forests (DNF; one site), evergreen broadleaf forests (EBF; two sites), open shrubland (OSH; one site), evergreen needle leaf forests (ENF; one site), grasslands (GRA; seven sites), and mixed forests (MF; one site). Our analysis is based on EC site data during the period of 2003-2011 (see Table 1). It has been noted that the sum of sensible heat (H) and LE as measured by the EC method is generally less than the available energy [43]. Therefore, the original ground-measured LE was corrected using the method put forward by Twine et al. (2000) and Jung et al. (2010) [19,44].

Observed Meteorological and Hydrological Data
A high spatial-temporal resolution (with a daily temporal resolution and a spatial resolution of 0.1 • ), gridded, near-surface meteorological dataset from the China Meteorological Forcing Dataset (CMFD) was used in this study. The dataset was made through fusion of remote sensing products, reanalysis of the dataset, and in situ observations at weather stations [45]. Monthly precipitation, relative humidity, radiation, and maximum      The time series of EVI and NDVI were acquired from moderate resolution imaging spectroradiometer (MODIS) products (MOD13A3), which provide 1 km spatial and monthly temporal resolutions. We used an average of four surrounding pixels around the EC flux sites to acquire the EVI and NDVI values. The LAI and FPAR products, potential factors influencing total ET segmentation, were extracted from MCD15A2H at 500 m spatial and 8-day temporal resolution. The 8-day temporal EVI, NDVI, LAI, and FPAR records were integrated into a monthly scale.

Canopy Height and Soil Moisture Data
Observations of canopy height were used to model plant sensitivity to water content at different depths. These canopy height (CH) and soil relative humidity (SRH) datasets of 376 agricultural gas stations spanning from 2003 to 2015 on the China meteorological data network (http://data.cma.cn (accessed on 11 June 2021)) were selected. Since the datasets only include the CH and SRH data of the cropland, the CH from the Geoscience Laser Altimeter System (GLAS) on the Ice, Cloud, and Land Elevation Satellite [46] and soil moisture obtained from the outputs of the GLDAS-2 Noah Land Surface Model (LSM) L4 model [47], which here-after is referred to as LDAS-derived SM, were employed to enrich the basic study data. Considering the different spatial-temporal resolution of each SM data product (see Table 2), in this study, using the multiple linear regression and trapezoidal method [48], monthly SM data from 2003 to 2015 of three layers (0-5, 5-20, and 20-100 cm) were reconstructed through integrating the observed precipitation, farmland SRH observation, SM dataset from Yang et al. (2020) [49], GRACE-derived ∆TWS [50], and GLDAS SM outputs at different depths (0-10, 10-40, and 40-100 cm).
The multiple linear regression used in our study is in terms of where SM is the soil moisture; SRH is the soil relative humidity, SRH = SM/SM field moisture capacity × 100%. During January 2002-December 2011, we define SM reconstructed = SM ITPLDAS , then the coefficients (a 1 , a 2 , a 3 , and a 4 ) are estimated, and SM reconstructed during January 2012-December 2015 are calculated. Note that the relationship between SRH observations and products was transplanted to the neighbored grid where there was a lack of observations. The trapezoidal method [51] used in this study can be described as: where the subscript ofθ refers to the soil layer (from i = 1 to k); k corresponds to the final depth of measurement; and ∆z is the depth interval between two successive measurements.

Spatial Distribution Characteristics of Soil Moisture
The estimation and prediction of SM, especially the root zone SM, have gained the interest of researchers [52,53]. Overall, many kinds of methods, including land surface modeling and the data assimilation technique, have been applied to estimate the multilayer SM. In this study, the monthly SM data from 2003 to 2015 of three layers (0-5, 5-20, and 20-100 cm) were reconstructed from multi-sources. Particularly, SM data collected at three observational networks (i.e., Ngari, Maqu, and Naqu; detailed information can be found in Yang et al. (2020)) across the Tibetan Plateau were employed to validate the reconstructed SM dataset [49]. Due to the lack of observation from deeper soils, Figure 2 only shows the comparison of the LDAS-derived and reconstructed monthly SM with the observed ones in the three layers of Naqu and surface layer of Ngari and Maqu. Both the observed station-averaged SM and its standard deviation are given in this figure. Only the LDAS-derived SM at a depth of 5 cm in the Ngari network is higher than the in situ data, as the LDAS-derived SM in the Maqu and Naqu networks are lower than the in situ data. This phenomenon might be due to the fact that the LDAS-derived SM only denoted liquid water and did not account for soil ice. In general, the reconstructed SM is much closer to the observation in all three networks. However, the Ngari network still exhibits underestimation in summer ( Figure 3A-C, we can also see that the SM values of individual depth are different, for example, the SM values of 0-5 cm in depth in southwest China are higher than those of 5-20 and 20-100 cm in depth. This region belongs to the tropical rainforest climate, the forest coverage is high, and the surface layer of soil has a strong capacity of water fixation [54]. A significant positive correlation between SM and P in eastern China has been identified in many previous studies [55,56]. Despite this, the relation between SM and P is a complex non-local process, and it is still unclear whether it represents the discovery of relevant phenomena or the explanation of the possible mechanism.  Figure 3A-C, we can also see that the SM values of individual depth are different, for example, the SM values of 0-5 cm in depth in southwest China are higher than those of 5-20 and 20-100 cm in depth. This region belongs to the tropical rainforest climate, the forest coverage is high, and the surface layer of soil has a strong capacity of water fixation [54]. A significant positive correlation between SM and P in eastern China has been identified in many previous studies [55,56]. Despite this, the relation between SM and P is a complex non-local process, and it is still unclear whether it represents the discovery of relevant phenomena or the explanation of the possible mechanism.

Comparison of Algorithm Parameters
The differences between the parameters of the PT-JPL algorithm and PT-SM algorithm lie in the core of the performances of these two algorithms. Therefore, before we evaluate the performances of the PT-JPL algorithm and PT-SM algorithm at EC sites and hydrological catchments, the model parameters, including soil moisture constraint on soil evaporation fsm, fREW, soil water availability constraint on transpiration fM, and fTRM are firstly examined and compared. As a typical site with different underlying surfaces, Fig

Comparison of Algorithm Parameters
The differences between the parameters of the PT-JPL algorithm and PT-SM algorithm lie in the core of the performances of these two algorithms. Therefore, before we evaluate the performances of the PT-JPL algorithm and PT-SM algorithm at EC sites and hydrological catchments, the model parameters, including soil moisture constraint on soil evaporation f sm , f REW , soil water availability constraint on transpiration f M , and f TRM are firstly examined and compared. As a typical site with different underlying surfaces, Figure 4 shows the variation plot between f sm and f REW , f M and f TRM calculated from the reconstructed SM data at the CN-Qia (forest) flux tower, CN-Ha2 (shrub) flux tower, CN-Du2 (grass) flux tower, and YC (crop) flux tower. It can be seen that for all biomes, almost all of the f TRM values are larger than the f M values, and the values of f REW are larger than those of f sm , especially in the valley values. Comparatively, the values of f REW are much larger than those of f sm for the forest and the shrub biomes, while those of f TRM are much larger than those of f M for the forest and the crop land cover types. This result is reasonable, as several studies have shown that in comparison with shorter vegetation, taller vegetation is less sensitive to soil water deficits since the deeper root can absorb water from deep soil or groundwater, where there is a delay in response to drought [23,57]. For instance, Figure 4G,H indicates that for the crop biome, the f REW series closely follows f sm , while in the peak values, the values of f REW are larger than those of f sm . For f TRM , it can also be clearly seen that the new soil water availability constraints on transpiration are always higher than the original f M , especially at the f M valley values (mostly occurring in the period from October to next March). Moreover, f TRM being considerably higher could be due to the fact that Yucheng tower is an agricultural experiment station that monitors areas in which wheat and corn crops are mainly planted in winter and summer, respectively.

Comparison of Incorporating Soil Moisture into Soil Evaporation and Transpiration
Since the algorithm models individual ET components, we can separately quantify the added value from incorporating SM into soil E and canopy T. Consequently, we evaluated the updated PT-SM algorithm with SM incorporation into soil E (PT_Es), the updated algorithm with SM incorporation into T (PT_T), and the updated algorithm with

Comparison of Incorporating Soil Moisture into Soil Evaporation and Transpiration
Since the algorithm models individual ET components, we can separately quantify the added value from incorporating SM into soil E and canopy T. Consequently, we evaluated the updated PT-SM algorithm with SM incorporation into soil E (PT_Es), the updated algorithm with SM incorporation into T (PT_T), and the updated algorithm with SM incorporation into both soil E and T (PT-SM) compared with the original PT-JPL model. Comparisons among the four EC sites, i.e., CN-Qia (forest) flux tower, CN-Ha2 (shrub) flux tower, CN-Du2 (grass) flux tower, and YC (crop) flux tower, are shown in Figure 5. It can be clearly seen that incorporating explicit SM can improve estimates of monthly ET series. Especially at the Yucheng EC site, the updated PT-SM algorithm shows greater  Figure 5D and Table 3). At the Yucheng EC site, by only replacing soil evaporation E s , R 2 and NSE increased by 2% and 7%, respectively, and RMSE and Bias decreased by 16% and 2 mm, respectively. By only replacing canopy T, R 2 and NSE increased by 4% and 15%, respectively, and RMSE and Bias decreased by 54% and 7 mm, respectively. According to our results, the performance of the updated PT-SM algorithm with SM incorporation into both E s and canopy T considerably improved (i.e., R 2 and NSE increased by 5% and 15%, respectively, and RMSE and Bias decreased by 56% and 9 mm, respectively; Table 3) when compared with the original PT-JPL-modelled ET. Moreover, observations between 60 and 120 mm are better represented by the updated model with a scatter closer to the 1:1 line from raised underestimation in ET (see Figure 5). Among all of the EC towers, we noted that, by only replacing soil evaporation E s , the R 2 and NSE increased by 2% and 7%, respectively, and the RMSE and Bias decreased by 20% and 3 mm, respectively. By only replacing canopy T, R 2 and NSE increased by 4% and 18%, respectively, and RMSE and Bias decreased by 34% and 7 mm, respectively (see Table 3). The improvement to the updated PT-SM algorithm is more obvious at the EC towers in the arid area (PET/P > 1.7) than that at the EC towers in the humid regions (PET/P < 1.7), with the R 2 and NSE being increased by 8% and 25%, respectively, and RMSE and Bias being decreased by 40% and 9 mm, respectively. In addition, both Figure 5 and Table 3 show that this improvement is strongly dependent on the simulation accuracy of vegetation T; the result also confirms that our approach of incorporating the soil moisture constraints into a simple algorithm can provide good estimates of ET series at different land cover regions.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 26 respectively; Table 3) when compared with the original PT-JPL-modelled ET. Moreover, observations between 60 and 120 mm are better represented by the updated model with a scatter closer to the 1:1 line from raised underestimation in ET (see Figure 5).

Model Evaluation with In Situ Forcing
Both the PT-JPL and PT-SM algorithms were executed for in situ evaluation. The in situ modeled ET from both the PT-JPL and PT-SM show good consistency with the observations (Figure 6, Table 3). We noted a site-wide average improvement in R 2 , NSE, RMSE, and absolute Bias as a result of the model improvements. Sites are ordered from wet to dry (left to right; top to bottom) based on the aridity index (PET/P). The PT-SM algorithm demonstrates greater performance than that of the original PT model at waterlimited sites (Table 3). On average, the PT-SM algorithm shows an apparent decrease in RMSE and absolute Bias (36% and 8 mm, respectively), a weak increase in R 2 (5%), and a large increase in NSE (20%) when compared with the original PT-JPL model. An interesting feature is that the more overall statistical improvements are observed at CN-Cng and YC sites, two sites that both have relatively dry conditions (PET/P = 2.38 and 1.82, respectively; Figure 6M,I) and a large fraction of ET comes from soil evaporation. It can also be seen that at YC, AR, and CN-Cng sites, the PT-SM model demonstrates considerable improvements to the model estimation of ET during the seasonal dry term ( Figure 6).
Compared with the observed 17 EC sites' monthly ET series, the R 2 and NSE of the PT-JPL monthly ET series are 0.94 and 0.90, respectively, and the RMSE and Bias are 10.05 and −6.34, respectively ( Figure 7A). The R 2 and NSE of the PT-SM monthly ET series are 0.97 and 0.96, respectively, and the RMSE and Bias are 5.90 and 2.13, respectively, when compared with the observed ET ( Figure 7B). Among all of the land covers (biomes), both models represent the observed monthly ET series well, especially at the EC sites across the MF and Crop biomes. Only at the EC sites in the GRA land cover type, the two models represent the observed monthly ET series not so well. Comparing Figure 7B to Figure 7A, we can see that the PT-SM model provides much better improvements than those of the original PT-JPL model at the EC sites across the EBF and GRA land cover types. This phenomenon illustrates that at the EBF and GRA land cover regions, the PT ET model, considering the soil water availability constraint on transpiration, is better at reflecting the monthly ET series.

Evapotranspiration Partitioning
As an important component of ET, the interception evaporation (E i ) of vegetation plays a vital role in water resources at both global and regional scales [23,57,58]. The ratios of canopy interception evaporation-to-precipitation (E i /P) values of different land covers are notably different, with the averages for all biomes ranging from 0.11 to 0.27, but all values are below 0.30 ( Figure 8 and Table 4). The average E i /P for the ENF land cover is higher than that of other types, and the OSH biome has the lowest value (0.11). This may be related to different raining conditions and LAIs across land cover types. The E i /P variation in the Crop biome is larger than that of other biomes, with values ranging from 0.11 to 0.27 and an average across the four Crop sites of 0.18. Research on the E i /P ratio has attracted much attention. Under Mediterranean climate conditions, Llorens and Domingo (2007) discovered that for an annual rainfall of 90-800 mm, the mean relative interception was approximately 18% for trees and 31.6% for shrubs [59]. Miralles et al. (2010) reviewed 42 studies from different periods and found average E i /P values for the EBF biome of 0.17 with ranges from 0.08 to 0.29, and the ENF biome average was 0.23 with ranges from 0.16 to 0.42 [57]. Gu et al. (2018) obtained the E i /P from 11 terrestrial biomes based on 75 eddy covariance towers and found average E i /P values for the Crop biome of 0.12 with ranges from 0.06 to 0.16, and the GRA biome average was 0.13 with ranges from 0.02 to 0.17 [58]. Compared with these previous studies, the overall variation is relatively consistent, but the E i /P of the individual biome is somewhat different, which is most likely due to the fact that the EC sites used in this study are relatively few in number (e.g., Crop: 11 sites in Gu et al. (2018) vs. four sites in our study; GRA: 12 sites in Gu et al. (2018) vs. seven sites in our study) [58]. Overall, we found that the E i /P variation in the land covers was large, ranging from 0.11 to 0.27, with an average of 0.16 among 17 sites (see Figure 8A), which is relatively reasonable compared with those of other similar studies (shown in Table 4). The canopy interception evaporation-to-evapotranspiration (E i /ET) ratios across different land cover types are shown in Figure 8B. The average E i /ET values across biomes range from 0.11 (OSH biome) to 0.32 (GRA biome). Large variations are noted within each land cover type, with the GRA biome having the greatest range (0.19-0.32) followed by the EBF biome (ranges of 0.20-0.28, respectively).  models represent the observed monthly ET series well, especially at the EC sites across the MF and Crop biomes. Only at the EC sites in the GRA land cover type, the two models represent the observed monthly ET series not so well. Comparing Figure 7B to Figure 7A, we can see that the PT-SM model provides much better improvements than those of the original PT-JPL model at the EC sites across the EBF and GRA land cover types. This phenomenon illustrates that at the EBF and GRA land cover regions, the PT ET model, considering the soil water availability constraint on transpiration, is better at reflecting the monthly ET series.

Evapotranspiration Partitioning
As an important component of ET, the interception evaporation (Ei) of vegetation plays a vital role in water resources at both global and regional scales [23,57,58]. The ratios of canopy interception evaporation-to-precipitation (Ei/P) values of different land covers are notably different, with the averages for all biomes ranging from 0.11 to 0.27, but all values are below 0.30 ( Figure 8 and Table 4). The average Ei/P for the ENF land cover is higher than that of other types, and the OSH biome has the lowest value (0.11). This may be related to different raining conditions and LAIs across land cover types. The Ei/P variation in the Crop biome is larger than that of other biomes, with values ranging from 0.11 to 0.27 and an average across the four Crop sites of 0.18. Research on the Ei/P ratio has attracted much attention. Under Mediterranean climate conditions, Llorens and Domingo (2007) discovered that for an annual rainfall of 90-800 mm, the mean relative interception was approximately 18% for trees and 31.6% for shrubs [59]. Miralles et al. (2010) reviewed 42 studies from different periods and found average Ei/P values for the EBF biome of 0.17 with ranges from 0.08 to 0.29, and the ENF biome average was 0.23 with ranges from 0.16 to 0.42 [57]. Gu et al. (2018) obtained the Ei/P from 11 terrestrial biomes based on 75 eddy covariance towers and found average Ei/P values for the Crop biome of 0.12 with ranges from 0.06 to 0.16, and the GRA biome average was 0.13 with ranges from 0.02 to 0.17 [58]. Compared with these previous studies, the overall variation is relatively consistent, but the Ei/P of the individual biome is somewhat different, which is most likely due to the fact that the EC sites used in this study are relatively few in number (e.g. our study) [58]. Overall, we found that the Ei/P variation in the land covers was large, ranging from 0.11 to 0.27, with an average of 0.16 among 17 sites (see Figure 8A), which is relatively reasonable compared with those of other similar studies (shown in Table 4). The canopy interception evaporation-to-evapotranspiration (Ei/ET) ratios across different land cover types are shown in Figure 8B. The average Ei/ET values across biomes range from 0.11 (OSH biome) to 0.32 (GRA biome). Large variations are noted within each land cover type, with the GRA biome having the greatest range (0.19-0.32) followed by the EBF biome (ranges of 0.20-0.28, respectively).     [58,60], that used 75 eddy covariance towers across a wide range of biomes. Previous studies [6,13,16] also indicate that a large variability occurs in T/ET, which is partly due to different ET measurements, model structures, and inconsistent spatial and temporal resolutions of forcing data. Table 4. Annual ratios (average) of canopy evaporation losses to the total precipitation based on land cover types.

Landcover
Interception Ratio Current Study Literature

Spatial Distributions of Mean Annual Evapotranspiration and Its Components
We applied the verified PT-SM algorithm to the whole of China for the period of 2003-2015 at a spatial resolution of 0.25° using CMFD meteorological data, MODIS products, and integrated SM datasets as described in Sections 3.2-3.4. Figure 9 shows the maps of mean annual ET, Ei, T, and Es over the period of 2003-2015. The mean annual ET distribution presented an intricate spatial structure with relatively high ET values in the southern China and low values located in the northwest. In most regions of the Xinjiang, Tibet, Qinghai, and Inner Mongolia provinces, the ET values are lower than 200 mm/year. In contrast, in the south of the Yangtze River and the Huai River, the ET values are higher than 800 mm/year. Generally, the mean annual ET increased from the northwest to the southeast ( Figure 9A). Almost all of the Ei values in Xinjiang, Qinghai, and Inner Mongolia provinces are lower than 50 mm/year, especially in the northwestern part of China (barren lands) ( Figure 9B). Similarly, T values of the northwestern part are lower than 100 mm/year, while only in the upper of Huai River basin, the Pearl River basin, central China, the southwest and southeast parts of China are the values of T higher than 600 mm/year. Notably, in the south region of Yunnan province, the T values are larger than 800 mm/year ( Figure 9C). Furthermore, Es values lower than 100 mm/year were found in the arid and desert regions of the northwest part of China, while in some parts of southern China, the Es values are higher than 400 mm/year (see Figure 9D).  The soil evaporation-to-evapotranspiration ratios across different biomes are shown in Figure 8D. The average E s /ET values across biomes range from 0.11 (EBF biome) to 0.49 (OSH biome). The largest variation was noted for the GRA biome, having a greatest range of 0.27-0.39, followed by the Crop biome (range of 0.26-0.35). Our estimated E s /ET values are within the range of previously reported values.

Spatial Distributions of Mean Annual Evapotranspiration and Its Components
We applied the verified PT-SM algorithm to the whole of China for the period of 2003-2015 at a spatial resolution of 0.25 • using CMFD meteorological data, MODIS products, and integrated SM datasets as described in Sections 3.2-3.4. Figure 9 shows the maps of mean annual ET, E i , T, and E s over the period of 2003-2015. The mean annual ET distribution presented an intricate spatial structure with relatively high ET values in the southern China and low values located in the northwest. In most regions of the Xinjiang, Tibet, Qinghai, and Inner Mongolia provinces, the ET values are lower than 200 mm/year. In contrast, in the south of the Yangtze River and the Huai River, the ET values are higher than 800 mm/year. Generally, the mean annual ET increased from the northwest to the southeast ( Figure 9A). Almost all of the E i values in Xinjiang, Qinghai, and Inner Mongolia provinces are lower than 50 mm/year, especially in the northwestern part of China (barren lands) ( Figure 9B). Similarly, T values of the northwestern part are lower than 100 mm/year, while only in the upper of Huai River basin, the Pearl River basin, central China, the southwest and southeast parts of China are the values of T higher than 600 mm/year. Notably, in the south region of Yunnan province, the T values are larger than 800 mm/year ( Figure 9C). Furthermore, E s values lower than 100 mm/year were found in the arid and desert regions of the northwest part of China, while in some parts of southern China, the E s values are higher than 400 mm/year (see Figure 9D).

Soil Water Availability Constraint on Evaporation and Transpiration
Many previous studies have tried to estimate soil moisture by measurements, such as time domain reflectometry soil moisture probes [61], hydrologic models [62], remote sensing techniques [63] or statistical learning tools such as support vector machines [64]. Unfortunately, most of these studies only paid attention to the SM at the surface. For vegetation with deeper roots, the transpiration water requirements also need the SM at the root zone. However, although the estimation of SM at the root zone has drawn great attention [53,65], its retrieval is difficult for spatial heterogeneity. Consequently, in this study, the multiple linear regression model and trapezoidal method were employed to estimate the SM of both the surface and root zone. The simulated SM performs the observation well (see Figure 2), and it was noted that the observation-driven model SM displays superior strength with coupling precipitation, soil relative humidity, ∆TWS, and GLDAS soil moisture outputs. These results further confirmed that for most regions, SM is closely related to precipitation and terrestrial water storage change, which is consistent with the findings of previous studies [34,53,65].
The performances of the updated PT-SM algorithm with the SM constraint separately on E s and T are different (see Table 3). Overall, it was observed that the updated PT-SM algorithm with the SM constraint on T (PT_T) performs the observation better than the PT-SM algorithm with the SM constraint on E s (PT_ Es), and the PT-SM model showed increased R 2 and NSE and reduced RMSE and Bias, with the greatest improvements occurring in water-limited regions. With the aid of SMAP-derived surface soil moisture data, Purdy et al. (2018) indicated that the performance of the updated PT-SM algorithm with the soil moisture constraint on E s has improved performance on T across the 14 Ameriflux EC sites distributed across the US [24]. This result is different with our study, a difference that may have resulted from the difference in f REW and f TRM estimates. In particular, some studies also demonstrated that SMAP might underestimate SM for the low biases in the Global Modeling and Assimilation Office surface temperature [39]. While Purdy et al. (2018) estimated the ET components with SMAP-derived soil moisture at a single depth (5 cm), in our study, SM data of three layers (0-5, 5-20, and 20-100 cm) were employed to estimate the ET components, which is a possible explanation of the difference. Figure 10 is plotted to investigate the spatial difference of E s and T estimation by the PT-JPL algorithm and PT-SM algorithm. In contrast to the T estimated in July 2003 by the PT-JPL algorithm ( Figure 10A), the estimated T results of the PT-SM algorithm showed higher values, especially in the central, southern, and northeastern parts of China (see Figure 10B), regions showing higher SM values (see Figure 3) and LAIs [66]. For E s estimation, the large differences also occurred in the northeastern, southern, and central parts of China ( Figure 10C,D). The results further illustrate that estimation of ET components can only be improved by adequate soil moisture.

Contribution of Each ET Component to the Total ET
We also evaluated the contribution of each ET component to the total ET by different algorithms for the whole of China. Figure 11 shows that the PT-SM algorithm considering SM estimated a higher T/ET value than those of the original PT-JPL model (64.9% vs. 53.7%) and lower Ei/ET and Es/ET values (21.2% vs. 15.4%; 25.1% vs. 19.7%). Overall, the average estimates of the T/ET value in our study are somewhat lower than previous studies' results, although the T/ET value was improved by incorporating SM of multilayers into the original PT-JPL model. In essence, there are two possible reasons behind this phenomenon: (1) In terms of the site scale, Ei of vegetation has rarely been separately considered by previous studies, as usually, only Es and T are observed, and, therefore, the T/ET estimates would be higher; (2) different research scales, monitoring techniques or algorithms may also contribute to the variability-for instance, the isotope-based approach constrained by hydrologic decoupling always overestimates T/ET [60]. Additionally, the PT-SM Es/ET and Ei/ET values of our study are larger than the similarly reported fractions from the GLEAM model by 7-15% and 11-12%, respectively [23,57], but similar with Es/ET, T/ET, and Ei/ET at 23 ± 1.7%, 54 ± 1.6%, and 21 ± 0.8% of total ET, respectively, as reported in Purdy et al. (2018). It is worth noting that all of these studies focus on the global scale.

Contribution of Each ET Component to the Total ET
We also evaluated the contribution of each ET component to the total ET by different algorithms for the whole of China. Figure 11 shows that the PT-SM algorithm considering SM estimated a higher T/ET value than those of the original PT-JPL model (64.9% vs. 53.7%) and lower E i /ET and E s /ET values (21.2% vs. 15.4%; 25.1% vs. 19.7%). Overall, the average estimates of the T/ET value in our study are somewhat lower than previous studies' results, although the T/ET value was improved by incorporating SM of multilayers into the original PT-JPL model. In essence, there are two possible reasons behind this phenomenon: (1) In terms of the site scale, E i of vegetation has rarely been separately considered by previous studies, as usually, only E s and T are observed, and, therefore, the T/ET estimates would be higher; (2) different research scales, monitoring techniques or algorithms may also contribute to the variability-for instance, the isotope-based approach constrained by hydrologic decoupling always overestimates T/ET [60]. Additionally, the PT-SM E s /ET and E i /ET values of our study are larger than the similarly reported fractions from the GLEAM model by 7-15% and 11-12%, respectively [23,57], but similar with E s /ET, T/ET, and E i /ET at 23 ± 1.7%, 54 ± 1.6%, and 21 ± 0.8% of total ET, respectively, as reported in Purdy et al. (2018). It is worth noting that all of these studies focus on the global scale. Figure 11C reveals the spatial patterns of dominant ET components. It can be seen that the bare soil evaporation is dominated across almost the entire arid area of northwest China, especially in the arid and semi-arid regions, while in some warmer and wet regions, such as southeastern of Tibet, the interception loss is the main component. In contrast, transpiration typically takes up most of the remaining regions (especially the forested terrain), while in some areas in southern China (often rainforests), the total interception loss is generally the dominant component. These findings are consistent with the partitioning obtained from other datasets [6,24]. Nonetheless, this phenomenon is a comprehensive reflection of climate, water absorption capacity of vegetation, growth status, soil moisture, etc. Remote Sens. 2021, 13, x FOR PEER REVIEW 20 of 26  Figure 11C reveals the spatial patterns of dominant ET components. It can be seen that the bare soil evaporation is dominated across almost the entire arid area of northwest China, especially in the arid and semi-arid regions, while in some warmer and wet regions, such as southeastern of Tibet, the interception loss is the main component. In contrast, transpiration typically takes up most of the remaining regions (especially the forested terrain), while in some areas in southern China (often rainforests), the total interception loss is generally the dominant component. These findings are consistent with the partitioning obtained from other datasets [6,24]. Nonetheless, this phenomenon is a comprehensive reflection of climate, water absorption capacity of vegetation, growth status, soil moisture, etc.

Inter-Annual Variability of ET for 2003-2015
It can be seen that the whole of China is characterized by complex spatial patterns in ET and its change in components ( Figure 12). For the total ET, more than half of the entire

Inter-Annual Variability of ET for 2003-2015
It can be seen that the whole of China is characterized by complex spatial patterns in ET and its change in components ( Figure 12). For the total ET, more than half of the entire area of China show upward trends, while the other half show downward trends. Additionally, the largest upward trend of the annual ET series is recognized at the Qinghai-Tibet Plateau with a value of 3.35 mm/y, while the largest downward trend (i.e., −2.88 mm/y) can be found at the southeastern part of China ( Figure 12A). Moreover, based on a modified Priestley-Taylor algorithm, Yao et al. (2013) found that the LE decreased over large areas in central China, northwest China, and Inner Mongolia while increasing in the northeast, north, and south regions of China during the period of 2001-2010 by driving monthly MODIS products [29]. These two results are different partly due to the different study periods, and another reason could be the difference in data and methods used. By employing the ensembles of remote sensing-based physical models and machine learning algorithms, overall increasing trends (0.62 and 0.38 mm/y, respectively) in global terrestrial ET during the period of 1982-2011 were evaluated by Pan et al. (2020) [17]. A further investigation into the control factors of ET variations showed that anthropogenic Earth greening was the dominating factor, which has also been reported by previous studies [1,67]. By contrast, drought and reducing solar radiation caused by substantial increases in aerosol optical depth (AOD) were usually identified as the main possible factors of ET reduction [68,69].
large areas in central China, northwest China, and Inner Mongolia while increasing in the northeast, north, and south regions of China during the period of 2001-2010 by driving monthly MODIS products [29]. These two results are different partly due to the different study periods, and another reason could be the difference in data and methods used. By employing the ensembles of remote sensing-based physical models and machine learning algorithms, overall increasing trends (0.62 and 0.38 mm/y, respectively) in global terrestrial ET during the period of 1982-2011 were evaluated by Pan et al. (2020) [17]. A further investigation into the control factors of ET variations showed that anthropogenic Earth greening was the dominating factor, which has also been reported by previous studies [1,67]. By contrast, drought and reducing solar radiation caused by substantial increases in aerosol optical depth (AOD) were usually identified as the main possible factors of ET reduction [68,69]. The ET components are not only influenced by vegetation cover and growth period but also by environmental factors, such as climate conditions, topography, and soil characteristics [70]. However, the control degree of different influencing factors on individual ET components varies with energy and water conditions. Ei across northwest China The ET components are not only influenced by vegetation cover and growth period but also by environmental factors, such as climate conditions, topography, and soil characteristics [70]. However, the control degree of different influencing factors on individual ET components varies with energy and water conditions. E i across northwest China indicates upward trends, and this may be a result of the afforestation and grassland projection over the past two decades. For almost the whole southeastern part of China, E i shows downward trends with the largest downward trend (i.e., −2.57 mm/y) being found in central China ( Figure 12B), this result is relatively reasonable due to the development of modern agricultural technology, a large number of agricultural water-saving technologies have been implemented in the southeast region, the water use efficiency was increased, and the unnecessary evaporation dissipation was reduced [71]. As for the vegetation T, most areas presented upward trends, except small areas in the northeastern and southeastern parts of the country, and the larger upward trends of annual T series were recognized at the Tibetan Plateau, with a value over 2.00 mm/y ( Figure 12C). Unlike the general paradox regarding whether or not the growing-season precipitation amount or the precipitation pattern is the main reason for vegetation T change [6,72], increased NDVI and LAI as a result of large-scale afforestation and agricultural intensification have been widely identified to play the most important role in the increase in vegetation T [17,58,73]. E s across China (with the exception of the southeastern region) presented downward trends. Furthermore, the changing trends were mostly less than 3.00 mm/y ( Figure 12D). Similar to E i , E s also represents the ineffective water consumption that does not directly contribute to production. Although the slopes of the trends were different, all of the values of ET and its components have been found to show increasing trends with increasing precipitation by Gu et al. (2018) [58]. Meanwhile, the variation of f REW (relative extractable water) has the most significant impact on the E s simulation in arid areas with low vegetation [74]. Nevertheless, such E s reduction may be due to the decrease in solar radiation or drought. Additionally, many studies noted that the sensibility of ET components changes in relation to different climatic environments, biomasses, soil properties, plant types, and canopy heights [7,23,75]. Therefore, the reasons for the distinction of the performances of the updated PT-SM algorithm with a soil moisture constraint on E s and T separately are complex and warrant further investigation.

Conclusions
By incorporating soil moisture into the widely used PT-JPL ET algorithm, in this study, ET and its components series at a spatial scale of 0.25 × 0.25 • were estimated. For this purpose, SM data at three layers (0-5, 5-20, and 20-100 cm) were firstly reconstructed by the multiple linear regression model and trapezoidal method, integrating the datasets of humidity observations relative to farmland soil, observed precipitation, GRACE solutions, and SM datasets obtained from the LDAS-2 outputs and Yang et al. (2020). High SM values generally occurred in the southeastern part of China, while low values were distributed in the northwestern part of the country. Meanwhile, the two algorithms' parameters in the update PT-SM ET algorithm (f REW and f TRM ) were compared with the parameters in the original PT-JPL algorithm (f sm and f M ). The values of parameters in the updated PT-SM ET algorithm (i.e., f REW and f TRM ) were observed to be larger than those of the original PT-JPL algorithm (i.e., f sm and f M ), especially in the peak or valley values.
The updated algorithm (i.e., PT-SM) model shows increased R 2 and NSE and reduced RMSE and Bias, with the greatest improvements occurring in water-limited regions. SM incorporation into E s can improve ET estimates by increasing R 2 and NSE by 2% and 7%, respectively; RMSE and Bias were reduced by 20% and 3 mm, respectively; while SM incorporation into T improved ET estimates by increasing R 2 and NSE by 4% and 18%, respectively. RMSE and Bias were reduced by 34% and 7 mm, respectively. These updated PT-SM ET estimates distinctly show a reduced error and provide a rich dataset to evaluate land surface models, vegetation or anthropogenic perturbation response to changes in water availability. The mean annual ET distribution presented an intricate spatial structure with relatively high ET values in the southeast of China and low values in the northwestern part. As a whole, the soil moisture constraint resulted in a higher transpiration estimate and lower evaporation estimate. These results can provide more accurate water consumption data for regional water resource management and planning. In the future, a high spatial-temporal resolution ET is expected to be obtained. Therefore, the critical variables determining the model's ability to simulate ET components, vegetation, surface, and root-zone soil moisture at high spatial-temporal resolution, are urgently needed.

Institutional Review Board Statement:
Not applicable for studies not involving humans or animals.

Informed Consent Statement:
Not applicable for studies not involving humans.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

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