Estimation of Evapotranspiration in Sparse Vegetation Areas by Applying an Optimized Two-Source Model

: Evapotranspiration (ET) is an important part of the water, carbon, and energy cycles in ecosystems, especially in the drylands. However, due to the particularity of sparse vegetation, the estimation accuracy of ET has been relatively low in the drylands. Therefore, based on the dry climate and sparse vegetation distribution characteristics of the drylands, this study optimized the core algorithms (canopy boundary resistance, aerodynamic resistance, and sparse vegetation coverage) and explored an ET estimation method in the Shuttleworth–Wallace two-layer model (SW model). Then, the Beijing–Tianjin sandstorm source region (BTSSR) was used as the study area to evaluate the applicability of the improved model in the drylands. Results show that: (1) The R 2 value of the improved model results was increased by 1.4 and the RMSE was reduced by 1.9 mm, especially in extreme value regions of ET (maximum or minimum). (2) Regardless of the spatial distribution and seasonal changes of the ET (63–790 mm), the improved ET estimation model could accurately capture the differences. Furtherly, the different vegetation regions could stand for the different climate regions to a certain extent. The accuracy of the optimized model was higher in the semi-arid region (R 2 = 0.92 and 0.93), while the improved model had the best improvement effect in the arid region, with R 2 increasing by 0.12. (3) Precipitation was the decisive factor affecting vegetation transpiration and ET, with R 2 value for both exceeding 0.9. The effect of vegetation coverage (VC) was less. This method is expected to provide a more accurate and adaptable model for the estimation of ET in the drylands. vegetation transpiration, soil water evaporation, and ET in the four seasons in key areas.


Introduction
Evapotranspiration (ET) plays a very important role in the ecosystem, connecting the global water, carbon, and energy cycles [1,2]. Water-limited ecosystems (including arid, semi-arid and sub-humid regions) account for more than 40% of the Earth's land surface [3]. In these areas, 90% of the effective precipitation returns to the atmosphere through ET. This causes water resources to be the main limiting factor for economic development and ecological protection in arid areas [4,5]. Meanwhile, the drylands experience the frequent occurrence of climatic natural disasters, such as soil degradation and droughts. Early warnings of drought can be achieved through the dynamic monitoring of long-term ET [3,[6][7][8][9]. Therefore, carrying out research on the accurate estimation of ET in drylands has great practical significance. Through this data, we can analyze the status of water profit, loss, and use efficiency. This can allow us to better guide agricultural production and ecosystem management [10,11].
The ET process is complex and observation sites are scarce; thus, it is difficult to precisely measure ET data obtained at large scales. However, remote sensing provides a cost-effective method to estimate regional or even global ET [12]. There are three main types of remote sensing methods used for estimating ET: The statistical empirical model

Study Area
The BTSSR (36 • 52 -46 • 44 N, 105 • 10 -121 • E) extends from Wulatehou Banner in Inner Mongolia in the west, to Aluhorqin Banner in Inner Mongolia in the east, to Dingbian County in Shaanxi in the south, and Dongwuzhumuqin Banner in Inner Mongolia in the north, covering 138 counties (banners, cities and districts) of six provinces (autonomous regions and municipalities): Beijing, Tianjin, Hebei, Shanxi, Inner Mongolia and Shaanxi. It has a total area of 710,500 square kilometers and a desertification area of 226,900 square kilometers, accounting for 31.9% of the total area of the region. A schematic diagram of the study area is shown in Figure 1.
The tree-grass ecosystems (TGEs) are a common ecological structure in the drylands. Vicente et al. [49] assumed that grassland and broadleaved trees were the dominant vegetation structure and cover separately in the growing period and dry period in semi-arid. On this basis, the leaf area index (LAI) input data of the original model were optimized to improve the model accuracy. However, the particularity of the distribution and ET process of the areas with sparse vegetation has not yet been fully grasped. Therefore, there are still large gaps in estimating ET of the drylands from biological process and sparse vegetation characteristics, including the canopy boundary resistance, aerodynamic resistance, vegetation coverage (VC), and so on.
The Beijing-Tianjin sandstorm source region (BTSSR) covers a wide area, and its climate and land-type changes are very complicated. Optimization of the ecological environment in this area can greatly improve the ecological environment in the Beijing-Tianjin-Hebei region, effectively curb the impact of meteorological disasters such as sandstorms, and slow down (or even improve) land degradation in central and northern China. The main objectives are: (1) Optimizing the core algorithm of the original SW model, in terms of the canopy boundary resistance, aerodynamic resistance, and regional ET estimation models, in order to adapt to the characteristics of dry climate and sparse vegetation in the drylands; and (2) Evaluating the applicability of the improved model in the drylands by analyzing the characteristics of spatial distribution and seasonal changes in vegetation transpiration, soil water evaporation, and ET in the BTSSR. This research is expected to improve the accuracy of the ET estimation in sparse vegetation areas and evaluate the water and energy cycle process of the ecosystem in the drylands.

Study Area
The BTSSR (36°52′-46°44′N, 105°10′-121°E) extends from Wulatehou Banner in Inner Mongolia in the west, to Aluhorqin Banner in Inner Mongolia in the east, to Dingbian County in Shaanxi in the south, and Dongwuzhumuqin Banner in Inner Mongolia in the north, covering 138 counties (banners, cities and districts) of six provinces (autonomous regions and municipalities): Beijing, Tianjin, Hebei, Shanxi, Inner Mongolia and Shaanxi. It has a total area of 710,500 square kilometers and a desertification area of 226,900 square kilometers, accounting for 31.9% of the total area of the region. A schematic diagram of the study area is shown in Figure 1.  The BTSSR landform consists of three regions: Plains, mountains and plateau. The area includes the Mu Us Sandy Land, the Kubuqi Desert, the Otindag Sandy Land and the Horqin Sandy Land. The western arid areas are mainly deserts and the Gobi, with bare land and low vegetation cover. The BTSSR includes a sub-humid arid zone, a semi-arid region, and an arid region, from southeast to northwest. The average annual precipitation decreases from about 450 mm in the southeast to about 150 mm in the northwest. The vegetation in the area is distributed sequentially, from the southeast to the northwest, as forest grassland, typical grassland, desert grassland, steppe desert, and desert, with large areas of sand vegetation distributed on the sandy land. According to the results of China's vegetation zone division, the study area includes eight vegetation regions; namely, (1) A northern temperate meadow steppe subzone; (2) A temperate shrub, tree and semi-desert subzone; (3) A typical steppe subzone in the northern temperate zone; (4) A warm-temperate northern deciduous oak forest zone; (5) A temperate southern desert steppe subzone; (6) A typical steppe subzone in the southern temperate zone; (7) A temperate southern forest (meadow) grassland zone; and (8) A temperate shrub, semi-shrub and desert subzone. Their distribution is shown in Figure 2. The BTSSR landform consists of three regions: Plains, mountains and plateau. The area includes the Mu Us Sandy Land, the Kubuqi Desert, the Otindag Sandy Land and the Horqin Sandy Land. The western arid areas are mainly deserts and the Gobi, with bare land and low vegetation cover. The BTSSR includes a sub-humid arid zone, a semi-arid region, and an arid region, from southeast to northwest. The average annual precipitation decreases from about 450 mm in the southeast to about 150 mm in the northwest. The vegetation in the area is distributed sequentially, from the southeast to the northwest, as forest grassland, typical grassland, desert grassland, steppe desert, and desert, with large areas of sand vegetation distributed on the sandy land. According to the results of China's vegetation zone division, the study area includes eight vegetation regions; namely, (1) A northern temperate meadow steppe subzone; (2) A temperate shrub, tree and semi-desert subzone; (3) A typical steppe subzone in the northern temperate zone; (4) A warm-temperate northern deciduous oak forest zone; (5) A temperate southern desert steppe subzone; (6) A typical steppe subzone in the southern temperate zone; (7) A temperate southern forest (meadow) grassland zone; and (8) A temperate shrub, semi-shrub and desert subzone. Their distribution is shown in Figure 2.

Data
The remote sensing data used were Moderate Resolution Imaging Spectroradiometer (MODIS, http://earthdata.nasa.gov, accessed on 25 February 2021) images, including MOD13A2, MOD15A2, MOD11A1, MCD43A1, MCD43B3 and MCD12Q1C. Table 1 shows the specific data parameter information. Preprocessing was carried out using the Google Earth Engine (GEE, https://code.earthengine.google.com, accessed on 25 February 2021) platform, including cropping and stitching actions. The remote sensing data were mainly used as the parameters of the SW two-layer model for calculation.

Data
The remote sensing data used were Moderate Resolution Imaging Spectroradiometer (MODIS, http://earthdata.nasa.gov, accessed on 25 February 2021) images, including MOD13A2, MOD15A2, MOD11A1, MCD43A1, MCD43B3 and MCD12Q1C. Table 1 shows the specific data parameter information. Preprocessing was carried out using the Google Earth Engine (GEE, https://code.earthengine.google.com, accessed on 25 February 2021) platform, including cropping and stitching actions. The remote sensing data were mainly used as the parameters of the SW two-layer model for calculation. The meteorological data were obtained from the China Meteorological Data Network (http://data.cma.cn, accessed on 25 February 2021). They include the daily value data of meteorological stations from 365 days in 2019, including average temperature, relative humidity, precipitation, average wind speed, and soil moisture content. The station data were differentiated using the Kriging space difference method, in order to form surface data.
Digital elevation model (DEM) data (with a spatial resolution of 90 m) were obtained from the NASA Earth Observing System Data and Information System (EOSDIS, https: //earthdata.nasa.gov/eosdis, accessed on 25 February 2021), for topographic correction when estimating the net surface radiation. According to the vectorized 1:1 million Chinese soil data provided by the Harmonized World Soil Database (HWSD) version 1.1, the saturated water content, field capacity, and wilting point of different soil types were calculated for the estimation of soil surface resistance and canopy resistance.
The correlation analysis data came from the Global Land Data Assimilation System (GLDAS). Bai et al. evaluated the accuracy and adaptability of GLDAS ET products in various parts of China, and their results showed that the ET error in China was 0.44 mm/year [50]. As the GLDAS ET products have a high accuracy in China, they were used as reference data to analyze the authenticity of the ET results of this study.

Methods
In this study, we optimized the original SW model, based on the characteristics of sparse vegetation and environmental differences in the drylands, in order to improve its ET estimation accuracy. The optimization model was mainly used due to the limitations of the original SW model, in terms of canopy boundary resistance, aerodynamic resistance, vegetation coverage (VC), and leaf area index (LAI). Figure 3 shows the research flow chart. Vegetation transpiration and soil water evaporation were calculated separately, where both required canopy boundary resistance and aerodynamic resistance as input parameters. The VC and LAI directly affect the net surface radiation, canopy boundary resistance, and aerodynamic resistance.

SW Model
In order to simulate the energy change in water flowing in different media (e.g., from soil water to plant water) and different state transitions (e.g., gaseous to liquid), the SW model introduces the concept of a resistance system. There are five types of resistance: aerodynamic resistance from canopy to reference height ( ; s·m −1 ), aerodynamic resistance from the surface to the canopy height ( ; s·m −1 ), canopy boundary resistance ( ; s·m −1 ), soil surface resistance ( ; s·m −1 ), and canopy stomatal resistance ( ; s·m −1 ). These resistances are closely related to the physiological characteristics of the vegetation (e.g., stomatal conductance, photosynthesis, and so on), the atmospheric environment of the vegetation (e.g., solar radiation, atmospheric temperature, water vapor pressure, and so on), and the process of water and heat transfer between soil, vegetation, and the atmosphere.
The ET (λET) of the ecosystem can be decomposed into vegetation transpiration (λT) and soil water evaporation (λE) through the calculation of energy conversion by the resistance system: where is the latent heat flux of the canopy (W·m −2 ), is the latent heat flux of the soil covered by the canopy (W·m −2 ), and are the distribution coefficients, D is the saturated vapor pressure difference (KPa), Δ is the change rate of saturated water vapor pressure difference with temperature, γ is the psychrometer constant (0.067 kPa·K −1 ), R is the total effective energy (W·m −2 ), is the effective ground energy (W·m −2 ), ρ is the air density, and is the constant pressure specific heat. The calculation formulae of the distribution coefficients and are as follows:

SW Model
In order to simulate the energy change in water flowing in different media (e.g., from soil water to plant water) and different state transitions (e.g., gaseous to liquid), the SW model introduces the concept of a resistance system. There are five types of resistance: aerodynamic resistance from canopy to reference height (r a a ; s·m −1 ), aerodynamic resistance from the surface to the canopy height (r s a ; s·m −1 ), canopy boundary resistance (r c a ; s·m −1 ), soil surface resistance (r s s ; s·m −1 ), and canopy stomatal resistance (r c s ; s·m −1 ). These resistances are closely related to the physiological characteristics of the vegetation (e.g., stomatal conductance, photosynthesis, and so on), the atmospheric environment of the vegetation (e.g., solar radiation, atmospheric temperature, water vapor pressure, and so on), and the process of water and heat transfer between soil, vegetation, and the atmosphere.
The ET (λET) of the ecosystem can be decomposed into vegetation transpiration (λT) and soil water evaporation (λE) through the calculation of energy conversion by the resistance system: where PM c is the latent heat flux of the canopy (W·m −2 ), PM s is the latent heat flux of the soil covered by the canopy (W·m −2 ), C s and C c are the distribution coefficients, D is the saturated vapor pressure difference (KPa), ∆ is the change rate of saturated water vapor pressure difference with temperature, γ is the psychrometer constant (0.067 kPa·K −1 ), R is the total effective energy (W·m −2 ), R s is the effective ground energy (W·m −2 ), ρ is the air density, and C p is the constant pressure specific heat. The calculation formulae of the distribution coefficients C s and C c are as follows: where: The calculation formulae of the total effective energy R and ground effective energy R s are as follows: where R n and R s n are the net radiant flux (W·m −2 ) of the canopy and the ground, respectively; G is the soil heat flux (W·m −2 ); and LAI is the leaf area index. According to the Lambert-Beer law, R s n is closely related to the leaf area index and vegetation type, and C is the attenuation coefficient of energy in the vegetation canopy community. In this study, different attenuation coefficients are substituted for different land-cover types.

Adaptability Optimization of SW Model in the Drylands
(1) Optimization of Canopy Boundary Resistance r c a The calculation of canopy boundary resistance in the original model simply divides the entire pixel into two parts-vegetation and bare soil-and the canopy boundary curve is the boundary between these two parts. This is an extreme situation. All vegetation gathers together and affects each other, so moisture is not easily lost to the outside air, resulting in increased canopy boundary resistance. In the drylands, the vegetation grows sparsely, and have little influence on each other. Therefore, the resistance of the canopy boundary is lower. If the calculation is still carried out according to the original model, there will be deviations ( Figure 4). The left side of Figure 4 shows the accumulation of different vegetation within the pixel. The original models unified these, according to the ideal assumption on the right. Obviously, the two situations on the left differ, in terms of their canopy boundary resistance.
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 20 where: The calculation formulae of the total effective energy R and ground effective energy are as follows: where and are the net radiant flux (W·m −2 ) of the canopy and the ground, respectively; G is the soil heat flux (W·m −2 ); and LAI is the leaf area index. According to the Lambert-Beer law, is closely related to the leaf area index and vegetation type, and C is the attenuation coefficient of energy in the vegetation canopy community. In this study, different attenuation coefficients are substituted for different land-cover types.

Adaptability Optimization of SW Model in the Drylands (1) Optimization of Canopy Boundary Resistance
The calculation of canopy boundary resistance in the original model simply divides the entire pixel into two parts-vegetation and bare soil-and the canopy boundary curve is the boundary between these two parts. This is an extreme situation. All vegetation gathers together and affects each other, so moisture is not easily lost to the outside air, resulting in increased canopy boundary resistance. In the drylands, the vegetation grows sparsely, and have little influence on each other. Therefore, the resistance of the canopy boundary is lower. If the calculation is still carried out according to the original model, there will be deviations ( Figure 4). The left side of Figure 4 shows the accumulation of different vegetation within the pixel. The original models unified these, according to the ideal assumption on the right. Obviously, the two situations on the left differ, in terms of their canopy boundary resistance.  Even under the same VC and LAI, due to the difference in the degree of leaf aggregation within the pixel, the canopy boundary resistance will still vary greatly. Therefore, we creatively established the escaping index (EI), in order to revise the canopy boundary resistance in this situation. The biggest influencing factor of the EI is the canopy structure and vegetation community structure. Under the same VC and LAI, the denser the canopy structure and vegetation community structure are, the lower the EI is; meanwhile, the sparser the canopy structure and vegetation community structure are, the higher the EI is. Thus, the EI can be used to optimize the canopy boundary resistance of the original model in order to improve the calculation accuracy. Due to the differing degrees of vegetation accumulation, the shadow occlusion in different directions is different, which leads to different reflectivity in different directions. Therefore, the bidirectional reflectance distribution function (BRDF) contains vegetation aggregation and dispersion information. The large-area EI can be obtained based on remote sensing data and through the normalized difference hot spots and dark spots (NDHD) [51,52]: where A and B are simulation coefficients, which vary with the solar zenith angle (SZA) (θ s ) and crown shape, and the definition of the NDHD is: where ρ h and ρ d represent the reflectivities of the hot and dark spots in the main plane, respectively. A hot spot is the peak of backscattering when the solar radiation is in the same direction as the viewing direction and there is no shadow. A dark spot is in the forward scattering direction away from a hot spot, where the largest shadow is observed, and corresponds to the smallest reflectivity in the main plane.
The ρ h and ρ d of each pixel are calculated from the BRDF product (MCD43) of the MODIS data. This product uses the RossThick-LiSparse-Reciprocal (RTLSR) model to calculate the anisotropic surface reflectance of each pixel [53,54]. The core is represented by three kernel functions: isotropic scattering, volume scattering and geometric optical surface scattering. The values of ρ h and ρ d are related to the vegetation coverage. However, there are many low-coverage vegetation areas in the study area and, so, it is difficult to obtain hot spots without background interference. Therefore, for areas with low vegetation coverage, the relative observation angle needs to be increased. In this study, when the vegetation coverage is greater than 25%, the relative observation angle is the image observation angle; when the vegetation coverage is less than 25%, the relative observation angle is 60 • [49].
The calculation formulae for the aerodynamic resistance inside the vegetation from leaf resistance (r b ) to canopy resistance (r c a ) are as follows [55][56][57]: where L is the LAI, w is the leaf width, and u h is the wind speed at the top of the canopy. In summary, the optimized formula for calculating the canopy boundary resistance is: (2) Optimization of Aerodynamic Resistance r a a and r s The aerodynamic resistance can be mainly divided into two parts: Closed vegetation and bare soil. The calculation formulae for the aerodynamic resistance of the closed vegetation area are shown in Equations (17) and (18): r a a (a) = where h r is the reference height, z 0 is the roughness, d is the zero-plane displacement, L is the LAI, h is the canopy height, u is the wind speed, and k is the Kaman constant. The calculation formulae for the aerodynamic resistance of bare soil are shown in Equations (19) and (20): The aerodynamic resistance revision in the original SW model is based on the LAI [37,54]. When LAI > 4, it is assumed that the vegetation is completely closed, while when LAI < 4, it is assumed that there is a linear correlation between the aerodynamic resistance and the LAI. However, there are significant differences among different types of vegetation, regarding leaf type, crown structure, and tree shape. Therefore, there will be great errors in dividing closed vegetation or not when using a single LAI. In response to this problem, in this paper, we used the porosity to quantify the canopy closure of vegetation. Each pixel was divided completely into closed vegetation and bare soil, based on the porosity. Then, the aerodynamic resistance was calculated for the closed vegetation and bare soil, respectively. The formulae are as follows: r a a = (1 − P)·r a a (a) + P·r a a (0).
The calculation formula for porosity is based on the Lambert-Beer formula: where L is the LAI of the pixel, G is the leaf correlation coefficient, and Ω is the aggregation index. (

3) Optimization of the ET Based on Poisson Distribution
The original SW model uses a single VC in each pixel to calculate the vegetation transpiration and soil water evaporation. However, in areas with sparse vegetation, the VC changes dramatically and different VCs may coexist within a 500 m pixel. For areas with sparse vegetation, vegetation growth can be regarded as random, and the value of each sub-pixel is a binary value of vegetation or bare soil. Therefore, from the perspective of the entire pixel, it can approximate a Poisson distribution. This can introduce the influence of different vegetation coverage into the model, thereby better characterizing areas with sparse vegetation by using a Poisson distribution to optimize the model. In this study, we divide the entire pixel into nine parts, each with a probability of 1-9 (VC × 10), where the VC × 10 of the entire pixel is used as the mean value of the Poisson distribution. The vegetation transpiration and soil water evaporation of a pixel are the sum of the vegetation transpiration and soil water evaporation of all sub-pixels. In this way, different vegetation coverage information can be fully utilized, and the characteristics of vegetation transpiration and soil water evaporation under different levels of vegetation coverage can also be expressed. The calculation formulae are as follows: where X i is the VC × 10 of sub-pixel i (with i ranging from 1 to 9), E c is the vegetation transpiration of the pixel, E c,i is the vegetation transpiration of sub-pixel i, E s is the soil water evaporation of the pixel, and E s,i is the soil water evaporation of sub-pixel i. Similarly, the LAI can also be approximated as a Poisson distribution. The method is the same as that of VC; the formulae are as follows: where X i is the LAI of sub-pixel i (with i ranging from 1 to 9), E c is the vegetation transpiration of the pixel, E c,i is the vegetation transpiration of sub-pixel i, E s is the soil water evaporation of the pixel, and E s,i is the soil water evaporation of sub-pixel i.

Relative Validation of Results
As it is difficult to determine the relative true value of the ET in a large region, we adopted the cross-analysis method, which is a common method used in the analysis of ET estimation results [58][59][60]. This method uses a reliable surface ET product, which has been tested against a large number of absolute true values [61]. In addition, the cross-analysis method can also obtain the correlation results on a surface scale. Therefore, it can show the spatial distribution, which is important in order for us to grasp the improved model effect in different regions. The GLDAS ET products were used as the reference data in this study, and the root mean square error (RMSE) and R-squared (R 2 ) were used to analyze the results. The calculation formulae are as follows: whereŷ i and y i are the reference data value and model-estimated value, respectively; y i is the average value of the reference data; and n is the total number of samples.

Changes Analysis of ET in the Vegetation Regions
The different vegetation regions could stand for the different climate regions to a certain extent. In order to better study the temporal and spatial characteristics of the ET and to explore the heterogeneity of ET in different vegetation regions, five statistics were calculated including minimum (min), maximum (max), mean (mean), median (med), and standard deviation (SD) of the eight vegetation regions and the entire study area. The minimum, maximum, and mean values were used to represent the numerical distribution of the vegetation regions, while the mean-median difference and standard deviation were used to characterize the heterogeneity of ET in this vegetation zone and to characterize the deviation when taking the mean value as the ET in the whole area.

Validation Results of the Model and the Reference Data
According to the analysis method introduced in Section 2.3.3, the effect of the optimization model is shown in Figure 5. Compared with the estimation results of the original model, different parts of the optimization model had different effects on the correlation coefficient. Among them, the R 2 of the original model was 0.68 and the RMSE was 8.2 mm. After the three-part optimization algorithm was carried out at the same time, the R 2 increased by 0.14 and the RMSE decreased by 1.9 mm. This indicates that the optimized model can effectively improve the estimation correlation of the original model. From the change in the scattered point distribution, it can be seen that the shape of the optimized model in Figure 5 was slenderer, which indicates that the optimized results were more stable and the results were more accurate. In particular, the optimization model was significantly optimized for the high-value area underestimated in the original model, while the low-value area was also significantly improved. Therefore, this indicates that the model has a significant improvement in extreme value regions (maximum or minimum) and has strong adaptability to the drylands.

Temporal and Spatial Distribution of ET in the BTSSR
According to the optimized SW model, the annual soil water evaporation, vegetation transpiration, and ET in the BTSSR in 2019 were obtained, as shown in Figure 6.
Because of the discontinuity and heterogeneity of vegetation in the drylands, the spatial distribution of the soil water evaporation, vegetation transpiration, and ET showed obvious differences. The improved model could accurately capture the differences. In the BTSSR, the highest value (261 mm) of soil water evaporation was observed in the northwest of the study area. The main land-cover type is Gobi in this area, with a low VC. The limited precipitation mainly dissipates into the atmosphere, in the form of soil water evaporation. On the other hand, the lowest value (19 mm) of soil water evaporation was in the southeast of the study area. The main land-cover type is forest in this area, with relatively closed VC. The precipitation is mainly dissipated in the form of vegetation transpiration. The regional distribution of vegetation transpiration showed an opposite trend to the regional distribution of soil water evaporation, being lower (>11 mm) in the northwest and higher (<573 mm) in the southeast. The highest values (790 mm) for ET were in the central and southern parts of the study area, while the lowest values (63 mm) were in the northwestern parts. It can be seen, from Figure 6, that no matter whether soil water evaporation, vegetation transpiration, or ET was considered, there were obvious spatial differences and gradual changes. Additionally, the higher the VC, the higher the vegetation transpiration, the higher the ET, and the lower the soil water evaporation. On the time scale, we analyzed the spatial distribution of soil water evaporation, vegetation transpiration, and ET over the four seasons (January-March, April-June, July-September and October-December) in key areas of the study area. The selection factor for the key areas was to have as many vegetation regions as possible. The changes in precipitation and vegetation distribution were also complicated, and the types of land-cover were also as diverse as possible. The results are shown in Figure 7.
From the seasonal changes in Figure 7, we can see that the vegetation transpiration, soil water evaporation, and ET were significantly higher in the second and third quarters than in the first and fourth quarters, with a magnitude difference of about 10 times. In addition, vegetation transpiration accounted for a large proportion of the ET-exceeding 80% in most periods and regions-and the proportions in the second and third quarters were significantly higher than those in the first and fourth quarters. The seasonal spatial distribution of the soil water evaporation had little change, and there was no obvious change in the high-and low-value areas. Due to the influence of perennial vegetation, the high-value areas of vegetation transpiration in the first and fourth quarters disappeared in the second and third quarters, and were replaced by areas with better precipitation conditions and a higher vegetation coverage. As vegetation transpiration accounts for a relatively high proportion of the ET, the seasonal spatial variation in ET was similar to that of vegetation transpiration. This demonstrates that the precipitation conditions and annual maximum VC have a greater impact on the regional vegetation transpiration and On the time scale, we analyzed the spatial distribution of soil water evaporation, vegetation transpiration, and ET over the four seasons (January-March, April-June, July-September and October-December) in key areas of the study area. The selection factor for the key areas was to have as many vegetation regions as possible. The changes in precipitation and vegetation distribution were also complicated, and the types of landcover were also as diverse as possible. The results are shown in Figure 7.
From the seasonal changes in Figure 7, we can see that the vegetation transpiration, soil water evaporation, and ET were significantly higher in the second and third quarters than in the first and fourth quarters, with a magnitude difference of about 10 times. In addition, vegetation transpiration accounted for a large proportion of the ET-exceeding 80% in most periods and regions-and the proportions in the second and third quarters were significantly higher than those in the first and fourth quarters. The seasonal spatial distribution of the soil water evaporation had little change, and there was no obvious change in the high-and low-value areas. Due to the influence of perennial vegetation, the high-value areas of vegetation transpiration in the first and fourth quarters disappeared in the second and third quarters, and were replaced by areas with better precipitation conditions and a higher vegetation coverage. As vegetation transpiration accounts for a relatively high proportion of the ET, the seasonal spatial variation in ET was similar to that of vegetation transpiration. This demonstrates that the precipitation conditions and annual maximum VC have a greater impact on the regional vegetation transpiration and ET. This is also the main reason why precipitation and VC were selected as influencing factors for the following analysis. In general, regardless of spatial or seasonal distribution, the improved model could accurately capture the differences in the drylands, which reflects the robustness of the improved model.

Temporal and Spatial Characteristics of ET in Different Vegetation Regions
The results of the five statistics are shown in Figure 8, in which the numbers 1-8 correspond to the vegetation regions in Figure 2 and the category "All" indicates the statistics referring to the entire study area.
Among the results, although the annual maximum occurred in the typical steppe subzone in the southern temperate zone, the mean and median of this area were lower than in the warm-temperate northern deciduous oak forest zone, while the mean-median differences and standard deviation were higher than those in this zone. Therefore, these results show that the overall ET was relatively high and the heterogeneity was low in the warm-temperate northern deciduous oak forest zone. The regional heterogeneity was relatively high and the distribution of ET was uneven in the typical steppe subzone in the southern temperate zone. On the other hand, the highest heterogeneity and the most uneven distribution of ET in this area was in the typical steppe subzone in the northern temperate zone. The smallest of the five statistics were observed in the temperate shrub, semishrub, and desert subzone, where the mean-median difference was also the lowest. This

Temporal and Spatial Characteristics of ET in Different Vegetation Regions
The results of the five statistics are shown in Figure 8, in which the numbers 1-8 correspond to the vegetation regions in Figure 2 and the category "All" indicates the statistics referring to the entire study area.
Among the results, although the annual maximum occurred in the typical steppe subzone in the southern temperate zone, the mean and median of this area were lower than in the warm-temperate northern deciduous oak forest zone, while the mean-median differences and standard deviation were higher than those in this zone. Therefore, these results show that the overall ET was relatively high and the heterogeneity was low in the warm-temperate northern deciduous oak forest zone. The regional heterogeneity was relatively high and the distribution of ET was uneven in the typical steppe subzone in the southern temperate zone. On the other hand, the highest heterogeneity and the most uneven distribution of ET in this area was in the typical steppe subzone in the northern temperate zone. The smallest of the five statistics were observed in the temperate shrub, semi-shrub, and desert subzone, where the mean-median difference was also the lowest. This indicates that the annual ET in this area was not only scarce, but also very evenly distributed.
Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 20 indicates that the annual ET in this area was not only scarce, but also very evenly distributed. On the time scale, Figure 9 shows the daily average changes in vegetation transpiration and soil water evaporation in the eight vegetation regions in 2019, and the ET could be calculated by adding the two. The northern temperate meadow steppe subzone reached 5.31 mm on 21 July, which was the highest regional daily average value; however, the seasonal variation in ET in this area was large. This is mainly because this area is a typical grassland area, where there is more annual vegetation. This led to the ET from July to August being higher and being lower in other months. The opposite was true in the warm-temperate northern deciduous oak forest zone, where there are many arbor forests and the seasonal changes in ET are not very dramatic. By analyzing the daily average changes in the temperate shrub, semi-shrub, and desert subzone, it can be seen that, although the annual ET in this area was relatively low, there were still high values of daily averages (e.g., exceeding 4 mm) on certain dates, due to the influence of precipitation. The spatial distribution of the annual ET in the northern temperate meadow steppe subzone was relatively uniform, but with large seasonal variation, indicating that the spatial and temporal characteristics of the region were not consistent.
The different vegetation regions could stand for different climate regions to a certain extent. The comparative analysis of the eight vegetation regions on the daily scale is shown in Figure 10. The red dots and line mean the comparative result between the original model with reference data, while the green ones mean the comparative result between the improved model with reference data. Among them, the northern temperate meadow steppe subzone and the warm-temperate northern deciduous oak forest zone are in the semi-arid region, while the temperate shrub, semi-shrub, and desert subzone is in the arid region, and other vegetation regions belong to the transitional areas between the arid and semi-arid regions. By analysis, the accuracy of the optimized model was higher in semiarid regions than in other regions, with R 2 values exceeding 0.9. In the arid region, the model had the best improvement effect, with R 2 increasing by 0.12. On the daily scale, the On the time scale, Figure 9 shows the daily average changes in vegetation transpiration and soil water evaporation in the eight vegetation regions in 2019, and the ET could be calculated by adding the two. The northern temperate meadow steppe subzone reached 5.31 mm on 21 July, which was the highest regional daily average value; however, the seasonal variation in ET in this area was large. This is mainly because this area is a typical grassland area, where there is more annual vegetation. This led to the ET from July to August being higher and being lower in other months. The opposite was true in the warm-temperate northern deciduous oak forest zone, where there are many arbor forests and the seasonal changes in ET are not very dramatic. By analyzing the daily average changes in the temperate shrub, semi-shrub, and desert subzone, it can be seen that, although the annual ET in this area was relatively low, there were still high values of daily averages (e.g., exceeding 4 mm) on certain dates, due to the influence of precipitation. The spatial distribution of the annual ET in the northern temperate meadow steppe subzone was relatively uniform, but with large seasonal variation, indicating that the spatial and temporal characteristics of the region were not consistent.
The different vegetation regions could stand for different climate regions to a certain extent. The comparative analysis of the eight vegetation regions on the daily scale is shown in Figure 10. The red dots and line mean the comparative result between the original model with reference data, while the green ones mean the comparative result between the improved model with reference data. Among them, the northern temperate meadow steppe subzone and the warm-temperate northern deciduous oak forest zone are in the semi-arid region, while the temperate shrub, semi-shrub, and desert subzone is in the arid region, and other vegetation regions belong to the transitional areas between the arid and semi-arid regions. By analysis, the accuracy of the optimized model was higher in semi-arid regions than in other regions, with R 2 values exceeding 0.9. In the arid region, the model had the best improvement effect, with R 2 increasing by 0.12. On the daily scale, the regions with the low ET (in the drylands) were the most difficult to estimate ET; However, the optimized model had the best improvement effect.

Influencing Factors on ET
From the input parameters of the SW model, it can be seen that there were three main types of factors affecting ET: The vegetation itself, soil topographic factors, and climate conditions. The relevant analysis diagram of influencing factors is shown in Figure 11. After the analysis, it can be seen that the annual precipitation had a strong influence on the vegetation transpiration and ET, with R 2 exceeding 0.9. This indicates that, in drylands, annual precipitation is the biggest influencing factor for vegetation transpiration and ET. The annual maximum NDVI (Normalized Difference Vegetation Index) also had a certain impact on ET, with an R 2 of 0.52. This was because areas with a higher VC generally have higher rainfall and ET; however, these correlations are indirect effects. The annual maximum NDVI had a greater impact on vegetation transpiration than ET. This was because the VC affects the vegetation transpiration to a certain extent, and vegetation transpiration accounts for a far greater proportion of ET than soil water evaporation. However, in drylands, the vegetation is more drought-tolerant and the ET may be quite different under the same VC. Therefore, in general, the annual precipitation has a strong effect on the ET and the effect of VC is less in this area.

Conclusions
Taking into account the complex characteristics of sparse vegetation in the drylands, we optimized the canopy boundary resistance, aerodynamic resistance, VC, and LAI of the SW model. Then, we comprehensively analyzed the spatial distribution and temporal changes of the vegetation transpiration, soil water evaporation, and ET in the BTSSR in 2019. Results demonstrated that: (1) The R 2 value of the improved model results was increased by 1.4 and the RMSE was reduced by 1.9 mm, especially in extreme value regions of ET (maximum or minimum). This indicates that the optimized model can effectively improve the estimation accuracy. (2) In general, the higher the VC, the higher the vegetation transpiration, the higher the ET, and the lower the soil water evaporation. Additionally, the seasonal spatial variation in soil water evaporation did not change much, while the vegetation transpiration and ET changed along with changes in precipitation and VC.
Regardless of the spatial distribution and seasonal changes of the ET (63-790 mm), the improved ET estimation model could accurately capture the differences in the drylands. Furtherly, the different vegetation regions could stand for different climate regions to a certain extent. By analysis, the accuracy of the optimized model was higher in semi-arid regions (R 2 = 0.92 and 0.93). In the arid region, the model had the best improvement effect, with R 2 increasing by 0.12. (3) In terms of influencing factors, precipitation was the decisive factor affecting vegetation transpiration and ET, with R 2 value for both exceeding 0.9. The effect of VC on the vegetation transpiration and ET in the region was less. The results show that the optimized SW model has high accuracy and applicability for the estimation of vegetation transpiration, soil water evaporation, and ET in drylands.

Discussion
Based on the climate and sparse vegetation characteristics of drylands, this study used the BTSSR as the study area and optimized the core algorithms, such as canopy boundary resistance, aerodynamic resistance, and sparse vegetation coverage, in the SW model. This method achieved the aim of further improving the estimation ability for regional ET in the drylands.
(1) The estimation results of the improved model are based on the BTSSR, and it is very site-specific. It may have the applicability in other drylands of the world which have the similar characteristics of climate and vegetation, but this needs to be verified.
Although the field data of regional ET is lacking to estimate the improved model, the spatial and temporal trends of ET in many previous studies were consistent with the results of this paper [62][63][64]. (2) When calculating the canopy boundary resistance, the improved model assumes that the vegetation source and the soil source are on the same plane, and ignores the water and energy exchange in the vertical direction. This may lead to underestimating the results. In addition, the improved model assumes that there is indeed an area of bare soil within the pixel, and makes it more accurate in the sparse vegetation. As for the areas with continuous vegetation, it may be a little difficult to parameterize the model, and there is need to judge in advance whether the pixel is a continuous vegetation pixel. (3) Considering the other factors simplified in this paper, the SW model may further improve the estimation accuracy of ET. First, the precipitation interception was not calculated, which may have caused the ET to seem higher. Second, we ignored the impact of snowmelt. There are freezing periods, ranging from the north to south in the study area, based on the geographical location. Although there is less precipitation in winter, the area is still affected by snowmelt. Zhao et al. found that snowmelt can cause the SW model to estimate vegetation transpiration to be too high and soil water evaporation to be too low [65]. Third, when discussing the analysis of factors affecting regional ET, we only considered the two most important factors-precipitation and VC-in this study. The next step should further analyze the vegetation itself, as well as atmospheric and surface environmental factors, in order to make the results more accurate.