Improved soil evaporation remote sensing retrieval algorithms and associated uncertainty analysis on the Tibetan Plateau

Actual evapotranspiration (ET) is the key link between water and energy cycles. However, accurate evaporation estimation in alpine barren areas remains understudied. In this study, we aimed to improve the satellite-driven Process-based Land Surface ET/Heat fluxes algorithm (P-LSH) for better satellite retrieval of ET on the Tibetan Plateau by introducing two effective soil moisture constraint schemes in which normalized surface soil moisture and the ratio of cumulative antecedent precipitation to cumulative antecedent equilibrium evaporation are used to represent soil water stress, respectively, based on the intercomparison and knowledge-learning of the existing schemes. We first conducted intercomparison of six existing soil evaporation algorithms and sorted out the two most effective soil moisture constraint schemes. We then introduced the modified versions of the two constraint schemes into the P-LSH algorithm and further optimized the parameters using the differential evolution method. As a result, it formed two improved P-LSH algorithms. We systematically assessed the performances of the two improved P-LSH algorithms and six existing remote sensing ET retrieval algorithms on two barren-dominated basins of the Tibetan Plateau using reconstructed ET estimates derived from the terrestrial water balance method as a benchmark. The two moisture constraint schemes largely improved the performance of the P-LSH algorithm and showed better performance in both basins (root mean square error (RMSE) = 7.36 and 7.76 mm per month; R²=0.86 and 0.87), resulting in a higher simulation accuracy than all six existing algorithms. We used five soil moisture datasets and five precipitation datasets to further investigate the impact of moisture constraint uncertainty on the improved P-LSH algorithm. The ET estimates of the improved P-LSH algorithm, driven by the GLDAS_Noah soil moisture, performed best compared with those driven by other soil moisture and precipitation datasets, while ET estimates driven by various precipitation datasets generally showed a high and stable accuracy. These results suggest that high-quality soil moisture can optimally express moisture supply to ET, and that more accessible precipitation data can serve as a substitute for soil moisture as an indicator of moisture status for its robust performance in barren evaporation.

1 Introduction

As a key link between the water and energy cycles, actual evapotranspiration (ET) is critical for assessing regional water and energy balances (Zhang et al., 2011). Oki and Kanae (2006) reported that approximately 60 % of precipitation returns to the atmosphere in the form of ET, whereas the proportion can reach more than 90 % in arid and semi-arid regions (Glenn et al., 2007; Morillas et al., 2013). Hence, accurate ET estimation is extremely important for irrigation planning, watershed management, and meteorology and climate change studies in arid and semi-arid regions.

Satellite remote sensing is an important means of estimating regional and global ET. A series of ET estimation algorithms have been developed over the past decade, including remote-sensing-based physical models, process-based land surface models, and vegetation-index-based empirical algorithms. In remote-sensing-based physical models, the Penman–Monteith (PM) method (Monteith, 1965; Cleugh et al., 2007; Mu et al., 2011; K. Zhang et al., 2010) and Priestley–Taylor (PT) method (Fisher et al., 2008; Martens et al., 2017; Priestley and Taylor, 1972; Yao et al., 2013) are the main representative methods for estimating ET. Several studies have combined these two methods to calculate canopy transpiration and soil evaporation (Leuning et al., 2008; Wang et al., 2018; Y. Zhang et al., 2019). The PT equation simplifies the PM equation and avoids the difficulty of quantifying aerodynamic and surface conductance. However, the PT equation simplifies the physical process, leading to a weaker physical basis than that of the PM equation. Land surface models reflect interactions and feedback between physical, biological, and biogeochemical processes in a predictive manner (Jiménez et al., 2011). These methods do not require remote sensing data; however, different parameterization schemes in land surface models for various physiological processes lead to considerable uncertainty in ET estimation (Famiglietti and Wood, 1991; Pan et al., 2020; Schwalm et al., 2013). In addition, ET has a close relationship with the ecophysiological processes that can be represented by satellite spectral products such as the normalized difference vegetation index (NDVI), leaf area index (LAI), and land surface temperature (LST). As a result, a number of vegetation-index-based empirical algorithms have been developed (Wang et al., 2006; Glenn et al., 2010). Subsequent developments in machine learning have attracted further attention in ET estimation because of their advantages in capturing the complex and nonlinear relationship between ET and its controlling environmental factors (Abdullah et al., 2015; Bai et al., 2021; Jung et al., 2010).

Although considerable effort has been made to estimate ET using the above methods, there are still significant uncertainties in quantifying the temporal and spatial characteristics and components of regional ET, especially in arid and semi-arid regions (Miralles et al., 2016; Pan et al., 2020). ET in these regions is dominated by water supply and climatic water deficits, whereas in humid regions it is dominated by available energy (Vinukollu et al., 2011; Zhang et al., 2016; Ma and Zhang, 2022). It is worth studying how to accurately reflect the influences of water supply and climatic water deficits. In remote-sensing-based physical models, both the PM and PT equations use the moisture constraint f to downscale the equilibrium (i.e., potential) evaporation at the soil surface to actual soil evaporation. Based on the hypothesis that surface moisture status is related to the adjacent atmospheric humidity (Bouchet, 1963), Fisher et al. (2008) used relative humidity (RH) and vapor pressure deficit (VPD) to reflect soil moisture supply and atmospheric water deficit and applied this method to a wide variety of ecosystems, vegetation types, footprints, and climatic regimes. Y. Zhang et al. (2019) selected the cumulative precipitation and cumulative equilibrium evaporation rates over the past 32 d to estimate f, based on which a continuous ET dataset was generated. Morillas et al. (2013) improved the method proposed by Y. Zhang et al. (2010) by adding a soil drying simulation factor after rainfall events and compared the uncertainties between three different methods in semi-arid and sub-humid flux towers in the Mediterranean. Miralles et al. (2011) also identified environmental factors that constrain potential evaporation by the moisture constraint f, parameterized for tall canopies, short vegetation, and barren areas. For barrens with sparse vegetation, the f estimates are based only on surface soil moisture (θ) conditions (Miralles et al., 2011; Martens et al., 2017), and soil moisture is normalized by the wilting point and critical moisture level, with an exponential (subsequently simplified to linear) form to estimate f. However, this method relies heavily on soil properties. Yao et al. (2013) incorporated diurnal temperature changes into apparent thermal inertia (ATI) estimation to calculate the moisture constraint f. This method was then compared with the relative extractable water (REW) of 16 flux towers in China and showed good agreement. García et al. (2013) also expressed the moisture constraint f using ATI, which was calculated using LST and albedo from the Meteosate Second Generation-Spinning Enhanced Visible and InfraRed Imager (MSG-SEVIRI) satellite. Their results showed that ET estimates derived from both towers and satellites performed better than the two-source model or the Penman–Monteith–Leuning model in the African Sahelian savanna and Mediterranean grasslands. However, this ATI-based method requires fine spatial and temporal resolutions of LST. Brust et al. (2021) calculated REW as moisture control directly, using soil moisture data from the NASA Soil Moisture Active Passive (SMAP) mission. Their results showed that the accuracy of the method with soil moisture control was better than that of the baseline MOD16. In summary, the f estimations proposed above performed well in their respective studies, but their applicability has not been sufficiently tested on barrens with sparse vegetation in arid or semi-arid basins, such as those found on the Tibetan Plateau.

Known as the “Asian Water Tower”, the Tibetan Plateau (TP) is crucial to the development of the Asian monsoon and water and energy cycles (Yao et al., 2012). Although great efforts have been made to evaluate ET in the sub-basin of the TP over the past few years (Xue et al., 2013; Hu et al., 2018; Wang et al., 2018; Li et al., 2019; Xu et al., 2018), most studies have focused on the headwaters of rivers in eastern or southern TP and have ignored the central and western inland arid and semi-arid regions. Ma et al. (2020) provided some hourly land-atmosphere interaction observations of inner regions with sparse vegetation; however, accurate soil evaporation estimates involving barrens remain a challenge. Li et al. (2014) reconstructed monthly ET estimates using the water balance method to evaluate five existing global ET products. They found that existing ET products were still not satisfactory for the Qaidam Basin and Qiangtang Plateau, two barrens-dominated sub-basins on the TP. In brief, the surface energy balance and land–atmosphere interaction mechanisms in alpine barren areas have not been explicitly revealed.

Therefore, we aim to discover the best mathematical representation of water supply, namely soil moisture, constraint to soil evaporation by learning from the existing schemes and further improve the satellite retrieval of soil evaporation in arid and semi-arid regions, especially in these understudied barren areas of the Tibetan Plateau. The specific objectives are: (1) to investigate the differences between the six existing soil evaporation algorithms and their applicability to alpine barren areas; (2) to improve the P-LSH algorithm by introducing two schemes for quantifying moisture constraints to ET in terms of surface soil moisture and precipitation, respectively; and (3) to test the applicability of satellite soil moisture and precipitation data for improving ET retrieval and analyze the influence of soil moisture and precipitation uncertainties on ET estimation on alpine barren areas.

2 Materials and study area

2.1 Study area

The Qaidam Basin is located in the northeastern TP (35^∘55^′–39^∘10^′ N, 90^∘00^′–98^∘20^′ E) and occupies an area of 257 768 km². The elevation of the Qaidam Basin is between 2676 and 6860 m, and the annual average temperature ranges from −6.4 to 14.5 ^∘C. Saline lakes and deserts cover approximately one-quarter and one-third of the Qaidam Basin, respectively. The Qaidam Basin has a typical continental climate with an average annual precipitation ranging from 29 to 387 mm, with approximately 80 % of the precipitation occurring in summer. Its drought conditions, high salinity, large diurnal and seasonal temperature ranges, and high ultraviolet radiation make the basin unsuitable for living. According to the MODIS IGBP classification (Friedl et al., 2010), 79.1 % of the Qaidam Basin is barren, 20.2 % is grassland, and other land uses/land cover types represent less than 1 %. Grassland is concentrated at the edge of the eastern and southern basins, whereas barren land is widely distributed across the remaining basins (Fig. 1a).

Figure 1 Locations and land cover/land use of (a) the Qaidam Basin and (b) the Qiangtang Plateau within China (OSH: open shrublands; GRA: grasslands; WET: wetlands; CRO: croplands; URB: urban and built-up lands; SNO: snow and ice; BAR: barren; WAT: water bodies).

The Qiangtang Plateau is located in the central hinterland of the TP, close to the Qaidam Basin. It forms the main feature of the TP with an area of 700 000 km². The average annual precipitation on the Qiangtang Plateau ranges from 50 to 300 mm in solid forms, such as snow, graupel, and hail, with precipitation being concentrated in the summer. The high altitude and inland surrounding high mountains make the Qiangtang Plateau a uniquely cold and arid region with widely distributed permafrost. Similar to the Qaidam Basin, barrens account for the largest proportion of the Qiangtang Plateau, reaching 55.7 %, whereas grassland and open water account for the second and third proportions, with values reaching 39.7 % and 3.0 %, respectively (Fig. 1b). The lakes on the Qiangtang Plateau cover an area of 21 400 km², accounting for approximately a quarter of all lake areas in China. The unique geographical structure makes the Qiangtang Plateau an endorheic area, which is also true of the Qaidam Basin, where water is retained and no outflow to other external rivers or oceans occurs. In an endorheic basin, drainage converges into inner lakes or swamps and equilibrates through evaporation.

2.2 Satellite and meteorological inputs

Table 1 summarizes the datasets used in this study. All input datasets were resampled from the original spatial resolution to a common $\mathrm{1}/\mathrm{12}$^∘ grid with a temporal resolution on a daily scale using the bilinear interpolation method. The daily meteorological inputs required by remote sensing algorithms are derived from the China Meteorological Forcing Dataset (CMFD) (He et al., 2020), including air temperature (T), specific humidity (q), air pressure (P_air), wind speed (u_m), and precipitation (P). The dataset incorporates existing reanalysis datasets and in situ observations, and shows better accuracy than existing reanalysis datasets (Yang et al., 2010; He et al., 2020). Radiation inputs come from the Clouds and the Earth's Radiant Energy System (CERES) SYN1deg radiative fluxes (Wielicki et al., 1996), which have provided continuous products since March 2000 with a resolution of 1^∘ globally. In this study, we used all-sky incoming shortwave radiation and net radiation. The NDVI product used in this study is from the MODIS MOD13Q1 Version 6 (https://lpdaac.usgs.gov, last access: 7 January 2023).

Table 1 List of the forcing datasets used in this study with their original resolutions and references.

In our algorithm, the surface soil moisture and precipitation were used to restrain soil evaporation. We selected various surface soil moisture and precipitation datasets from satellites, microwave assimilation, machine-learning methods, and reanalysis. The surface soil moisture comes from five datasets including the following: (i) the soil moisture dataset of China based on microwave data assimilation (Yang et al., 2020) (denoted as θ_Yang in this study); (ii) the land surface soil moisture dataset of SMAP time-expanded daily $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ over the Qinghai–Tibet Plateau Area (Qu et al., 2019) (denoted as θ_Qu); (iii) the combined product from the European Space Agency's Climate Change Initiative (ESA CCI) Soil Moisture Version 06.1 (Gruber et al., 2019) (denoted as θ_ESA CCI); (iv) Global Land Data Assimilation System (GLDAS) Noah Land Surface version 2.1 (Rodell et al., 2004) (denoted as θ_GLDAS Noah); and (v) the second Modern-Era Retrospective Analysis for Research and Applications (MERRA) version 2 (Molod et al., 2015) (denoted as θ_MERRA). The precipitation comes from five datasets including the following: (i) CMFD (denoted as P_CMFD); (ii) Global Precipitation Measurement (GPM) IMERG Final Precipitation L3 Version 06 (Hou et al., 2014) (denoted as P_GPM); (iii) Multi-Source Weighted-Ensemble Precipitation (MSWEP) version 2.8 (Beck et al., 2019) (denoted as P_MSWEP); (iv) GLDAS Noah (denoted as P_GLDAS Noah); and (v) MERRA (denoted as P_MERRA). All of the above soil moisture and precipitation sequences were resampled to $\mathrm{1}/\mathrm{12}$^∘.

Our algorithm adopts different parameterization schemes according to pixelated land cover, which comes from the MODIS Land Cover Type Yearly L3 Global 500 m SIN Grid (MCD12Q1) (Friedl et al., 2010). The MCD12Q1 product provides land cover properties, which come from observations spanning 1 year from the Terra and Aqua satellites. Here, we used data from 2003 and regarded them as static values. We calculated the percentage of various land covers for each pixel ($\mathrm{1}/\mathrm{12}$^∘), estimated the ET of various land covers, and then weighted each pixel by the percentage. Soil properties, including residual soil moisture and saturated water content, were obtained from the China Dataset of Soil Hydraulic Parameters Pedotransfer Functions for Land Surface Modeling (Dai et al., 2013). We aggregated the dataset from the original 30^′′ resolution to $\mathrm{1}/\mathrm{12}$^∘ using the arithmetic averaging method.

To evaluate the robustness and uncertainty of various remote-sensing algorithms, this study used reconstructed ET estimates derived from the terrestrial water balance method (ET_recon) as a benchmark. For endorheic basins, river discharge is zero, and ET is equal to the residue between precipitation and change in terrestrial water storage (ΔS). Based on this method, Li et al. (2014) established a monthly ET_recon for the Qaidam Basin and Qiangtang Plateau from 2003 to 2012. The gridded precipitation data for this study were obtained from the National Meteorological Information Center of the China Meteorological Administration (CMA), and ΔS was obtained from Gravity Recovery and Climate Experiment (GRACE) land data.

3 Methodology

3.1 Description of the baseline algorithm: P-LSH

The Process-based Land Surface Evapotranspiration/Heat Fluxes (P-LSH) algorithm (K. Zhang et al., 2009, 2010) is an ET algorithm evolved from the PM equation, in which canopy conductance comes from the Jarvis–Stewart formula (Jarvis, 1976; Stewart, 1988) and an empirical g₀-NDVI equation (K. Zhang et al., 2009, 2010). The P-LSH algorithm distinguishes between open water and vegetation pixels using land cover classification. Vegetation pixels include canopy transpiration and soil evaporation, whereas open water pixels only contain water evaporation.

(1) Canopy transpiration

The P-LSH algorithm calculates canopy transpiration (E_c: mm) by a modified PM equation,

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

where λ (J kg⁻¹) is the latent heat of vaporization, Δ (Pa K⁻¹) is the slope of the curve relating saturated water vapor pressure to air temperature, VPD (Pa) is the vapor pressure deficit, ρ (kg m⁻³) is the air density, C_p (J kg⁻¹ K⁻¹) is the specific heat capacity of air, γ (–) is the psychrometric constant, A_c (W m⁻²) is the available energy component allocated to the canopy based on fractional vegetation cover, and g_{a_c} (m s⁻¹) is the aerodynamic conductance of the canopy. Based on various vegetation types, K. Zhang et al. (2010) established an empirical relationship between the maximum canopy conductance (g₀: m s⁻¹) and NDVI based on observations from flux towers and reduced conductance from the maximum (g₀: m s⁻¹) to the actual value (g_c: m s⁻¹) through restraints from T (^∘C, VPD (Pa), and CO₂ (ppm). Feng et al. (2022) added incoming shortwave radiation and surface soil moisture to strengthen restraints on g_c over three TP grasslands. More details regarding canopy transpiration are available in Feng et al. (2022) and Zhang et al. (2015).

(2) Soil evaporation

The P-LSH algorithm combines the modified PM equation and complementary relationship hypothesis to quantify soil evaporation (E_s: mm) (Bouchet, 1963; Fisher et al., 2008), which can be expressed as

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{s}}=f{E}_{\mathrm{eqs}},\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{eqs}}={\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_{\mathrm{p}}\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}/g_{\mathrm{totc}}}},\text{(4)}& {\displaystyle }& {\displaystyle }f={\mathrm{RH}}^{\frac{\mathrm{VPD}}{k}},\end{array}$$

where f (–) is the moisture constraint, RH (–) is the relative humidity, k (Pa) is a parameter to fit the complementary relationship, E_eqs (mm) is the equilibrium (i.e., potential) evaporation, A_s (W m⁻²) is the available energy component allocated to the soil surface, and g_totc (m s⁻¹) is the corrected value of g_tot (m s⁻¹) based on the standard temperature and pressure. In this study, the g_tot term was expressed in the form of resistance r_tot ($r_{\mathrm{tot}}=\mathrm{1}/g_{\mathrm{tot}}$: s m⁻¹) and g_{a_s} (m s⁻¹) is the aerodynamic conductance of the soil surface. More details regarding soil evaporation are available in Mu et al. (2007) and K. Zhang et al. (2010).

(3) Open water

For open water pixels, the P-LSH algorithm uses the Penman equation rewritten by Shuttleworth (1993) to quantify the effects of the surface wind speed on open water evaporation (E_w: mm). The surface resistance r_s (s m⁻¹) is assumed to be zero on the open water surface; therefore, the PM equation is revised as

$$\begin{array}{}\text{(5)}& \mathit{\lambda }{E}_{\mathrm{w}}={\displaystyle \frac{\mathrm{\Delta }A+\mathit{\rho }{C}_{\mathrm{p}}\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{w}}}{\mathrm{\Delta }+\mathit{\gamma }}},\end{array}$$

where A (W m⁻²) is the available energy component for open water, following K. Zhang et al. (2010). The g_{a_w} (m s⁻¹) term is the aerodynamic conductance of the open water and is estimated by the wind speed,

$$\begin{array}{}\text{(6)}& g_{\mathrm{a}\mathrm{\_}\mathrm{w}}={\displaystyle \frac{\mathrm{1}+\mathrm{0.536}{U}_{\mathrm{2}}}{\mathrm{4.72}{\left[\ln \left(z_{\mathrm{m}}/z_{\mathrm{0}}\right)\right]}^{\mathrm{2}}}},\end{array}$$

where U₂ (m s⁻¹) is the 2 m height wind speed calculated from the reanalysis data and the vertical wind speed function, z_m (m) is the wind measurement height, and z₀ (m) is the aerodynamic roughness of the water surface, which is set to 0.00137.

3.2 Five existing soil evaporation algorithms

In this study, we further selected the soil evaporation schemes from five existing ET algorithms, including the Penman–Monteith–Leuning (PML) algorithm (Y. Zhang et al., 2010, 2019), Global Land Evaporation Amsterdam Model (GLEAM) algorithm (Martens et al., 2017), the Priestley Taylor-Jet Propulsion Laboratory (PT-JPL) algorithm (Fisher et al., 2008), the Priestley Taylor-Yao (PT-Yao) algorithm (Yao et al., 2013), and the Penman–Monteith–Brust (PM-Brust) algorithm (Brust et al., 2021).

(1) PML soil evaporation algorithm

The PML algorithm quantifies soil evaporation using the modified PT equation, which avoids the difficulty of parameterizing the resistances in the PM equation (Y. Zhang et al., 2010, 2019),

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{s}}=f{E}_{\mathrm{eqs},n},\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{eqs},n}={\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }}},\end{array}$$

where A_s (W m⁻²), Δ (Pa K⁻¹), and γ (–) represent the same physical meanings as in Eq. (3). The moisture constraint f (–) is estimated by the cumulative precipitation and equilibrium evaporation in the previous periods, without any observation of soil moisture as input,

$$\begin{array}{}\text{(9)}& f=\min \left({\displaystyle \frac{\sum _{n=\mathrm{1}}^N{P}_n}{\sum _{n=\mathrm{1}}^N{E}_{\mathrm{eqs},n}}},\mathrm{1}\right),\end{array}$$

where P_n (mm) and E_eqs,n (mm) are the cumulative precipitation and equilibrium evaporation of the surface in the nth period, respectively, and N is the number of periods.

(2) GLEAM soil evaporation algorithm

Similar to the PML algorithm, GLEAM takes the PT equation as the equilibrium soil evaporation and reduces it to actual soil evaporation through the moisture constraint f (Martens et al., 2017). The difference is that the GLEAM algorithm estimates f for tall canopies, short vegetation, and barren areas. For barren areas with sparse vegetation, the surface soil moisture is linearized by the critical moisture level and residual soil moisture, and is then used to estimate soil evaporation, which is expressed as

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{s}}=f{E}_{\mathrm{p}},\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{p}}=\mathit{\alpha }{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{A}_{\mathrm{s}},\text{(12)}& {\displaystyle }& {\displaystyle }f=\mathrm{1}-{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{c}}-\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}-{\mathit{\theta }}_{\mathrm{r}}}},\end{array}$$

where f (–) is the same as that in Eq. (9) to explain the restraints of the suboptimal environment on soil evaporation; E_p (mm) is the potential soil evaporation; α (–) is the PT dimensionless coefficient, and 1.26 for barrens; θ (cm³ cm⁻³) is the actual surface soil moisture; and θ_c (cm³ cm⁻³) is the critical moisture level and is set as θ_c=0.75θ_s following Zhu et al. (2013), where θ_s (cm³ cm⁻³) is the saturated water content and θ_r (cm³ cm⁻³) is the residual soil moisture.

(3) PT-JPL soil evaporation algorithm

The PT-JPL algorithm uses the same equilibrium soil evaporation as the GLEAM algorithm, with the difference being in the f estimation (Fisher et al., 2008). In the PT-JPL algorithm, f is constituted by f_SM (–) and f_wet (–), where f_SM comes from RH and VPD (the same as in the P-LSH algorithm), whereas f_wet is only determined by RH as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{SM}}={\mathrm{RH}}^{\frac{\mathrm{VPD}}{k}},\text{(14)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{wet}}={\mathrm{RH}}^{\mathrm{4}},\text{(15)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{s}}=\left[f_{\mathrm{wet}}+f_{\mathrm{SM}}\left(\mathrm{1}-f_{\mathrm{wet}}\right)\right]E_{\mathrm{p}},\end{array}$$

where k (Pa) is a parameter to fit the complementary relationship, which is calibrated using the bisection method in this study since it is the only parameter in the algorithm and lacks prior values for barren. The E_p (mm) is the equilibrium soil evaporation calculated using Eq. (11).

(4) PT-Yao soil evaporation algorithm

Yao et al. (2013) used the diurnal land surface temperature range (DTsR: ^∘C) and air temperature range (DTaR: ^∘C) to simplify the calculation of the apparent thermal inertia (ATI: ^∘C⁻¹) for f_SM (–) estimation with equilibrium soil evaporation using the PT equation, same as the GLEAM and PT-JPL algorithms:

$$\begin{array}{}\text{(16)}& f_{\mathrm{SM}}={\mathrm{ATI}}^k={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{DT}}}\right)}^{\mathrm{DT}/{\mathrm{DT}}_{\max }},\end{array}$$

where DT_max (^∘C) is defined as the maximum daily temperature range (DT: ^∘C), which reflects the relative sensitivity to changes in the daily temperature range and is set as a constant (DTaR_max=40 ^∘C, DTsR_max=60 ^∘C). Yao et al. (2013) showed that the performances of soil evaporation from DTaR and DTsR are similar; therefore, in this study, we only used DTaR for f_SM estimation.

(5) PM-Brust soil evaporation algorithm

The PM-Brust algorithm (Brust et al., 2021) originated from the MOD16 algorithm that is based on the PM equation (Mu et al., 2011). The equilibrium soil evaporation in the PM-Brust algorithm is similar to Eq. (3), with the resistance estimations slightly different from those of the P-LSH algorithm. The PM-Brust algorithm assumes that the boundary layer resistance is equal to the total aerodynamic resistance (r_tot: s m⁻¹), which is determined by VPD and four biome-specific constants, including maximum resistance (rbl_max: s m⁻¹), minimum resistance (rbl_min: s m⁻¹), VPD at which canopy stomata are completely open (VPD_open: Pa), and VPD at which canopy stomata are completely closed (VPD_close: Pa). In contrast, the P-LSH algorithm assumes that the boundary layer resistance and total aerodynamic resistance are biome-specific constants. Brust et al. (2021) estimated f_SM with a more direct soil moisture control outline (i.e., REW):

$$\begin{array}{}\text{(17)}& {\displaystyle }{f}_{\mathrm{SM}}=\mathrm{REW}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\min }}{{\mathit{\theta }}_{\max }-{\mathit{\theta }}_{\min }}},\text{(18)}& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{eqs}}={\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_{\mathrm{p}}\left(\mathrm{1}-f_{\mathrm{c}}\right)\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }\times g_{\mathrm{a}\mathrm{\_}\mathrm{s}}/g_{\mathrm{totc}}}},\text{(19)}& {\displaystyle }{g}_{\mathrm{totc}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{tot}}\times r_{\mathrm{corr}}}},\text{(20)}& {\displaystyle }\begin{array}{rl}& r_{\mathrm{tot}}=\\ & \left\{\begin{array}{ll}{\mathrm{rbl}}_{\max }& \mathrm{VPD}\le {\mathrm{VPD}}_{\mathrm{open}}\\ {\mathrm{rbl}}_{\max }-\frac{\left({\mathrm{rbl}}_{\max }-{\mathrm{rbl}}_{\min }\right)\times \left({\mathrm{VPD}}_{\mathrm{close}}-\mathrm{VPD}\right)}{{\mathrm{VPD}}_{\mathrm{close}}-{\mathrm{VPD}}_{\mathrm{open}}}& \\ {\mathrm{rbl}}_{\min }& \mathrm{VPD}\ge {\mathrm{VPD}}_{\mathrm{close}}\end{array}\right.\\ & {\mathrm{VPD}}_{\mathrm{open}}<\mathrm{VPD}<{\mathrm{VPD}}_{\mathrm{close}},\end{array}\text{(21)}& {\displaystyle }{E}_{\mathrm{s}}=\left[f_{\mathrm{wet}}+f_{\mathrm{SM}}\left(\mathrm{1}-f_{\mathrm{wet}}\right)\right]E_{\mathrm{eqs}},\end{array}$$

where REW (–) is the relative extractable water, θ (cm³ cm⁻³) is the surface soil moisture, θ_min (cm³ cm⁻³) and θ_max (cm³ cm⁻³) are the minimum and maximum values of θ for the period of record, respectively, and f_c (–) is the vegetation cover fraction. The g_{a_s} (m s⁻¹) is the aerodynamic conductance of the soil surface and is the sum of the conductance to radiative heat transfer, which is calculated using the same method proposed by Choudhury and DiGirolamo (1998), and the conductance to convective heat transfer, which is set equal to r_tot following Mu et al. (2011). The r_corr (–) is the correction coefficient of r_tot following Mu et al. (2011). In this study, the parameters rbl_max and rbl_min are calibrated using the Differential Evolution method (Storn and Price, 1997), while the parameters VPD_open and VPD_close are set to 650 and 4200 Pa following Mu et al. (2011) because they are relatively insensitive parameters (K. Zhang et al., 2019; Feng et al., 2022).

3.3 Improvements in the P-LSH soil evaporation algorithm

We attempted two strategies to improve soil evaporation in the P-LSH algorithm. One strategy was to directly control f through the surface soil moisture as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathit{\lambda }{E}_{\mathrm{s}}=& \left[f_{\mathrm{wet}}+{\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\min }}{{\mathit{\theta }}_{\max }-{\mathit{\theta }}_{\min }}}\left(\mathrm{1}-f_{\mathrm{wet}}\right)\right]\\ & {\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_{\mathrm{p}}\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}/g_{\mathrm{totc}}}},\end{array}\end{array}$$

where each item has the same meaning as that in Eqs. (3) and (17). The only parameter r_tot is calibrated using the bisection method.

The second strategy was to use the ratio of cumulative precipitation to equilibrium evaporation in the previous periods to quantify moisture constraint, with equilibrium evaporation estimated by the modified PM equation as follows:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\mathit{\lambda }{E}_{\mathrm{s}}=& \min \left({\displaystyle \frac{\sum _{n=\mathrm{1}}^N{P}_n}{\sum _{n=\mathrm{1}}^N\frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_{\mathrm{p}}\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}/g_{\mathrm{totc}}}}},\mathrm{1}\right)\\ & {\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_{\mathrm{p}}\mathrm{VPD}{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}}{\mathrm{\Delta }+\mathit{\gamma }{g}_{\mathrm{a}\mathrm{\_}\mathrm{s}}/g_{\mathrm{totc}}}},\end{array}\end{array}$$

where each item has the same meaning as that in Eqs. (3) and (9).

We combined the vegetation evapotranspiration (i.e., vegetation transpiration and vegetation surface evaporation) and open water evaporation components with the new soil evaporation component based on the first strategy to form an improved P-LSH algorithm, which is called P-LSH_θ. Similarly, we built the second improved P-LSH algorithm based on the second strategy (hereafter called P-LSH_P). By contrast, the original P-LSH soil evaporation algorithm is called P-LSH_ori in this study.

3.4 Evaluation of algorithm performance

Because we do not have direct observation of soil evaporation, we have to rely on the ET_recon as the benchmark to assess our improved soil evaporation algorithms and their associated ET retrieval algorithms. Therefore, we need to assemble the soil evaporation algorithm with the vegetation evapotranspiration and water evaporation algorithms to form a complete ET retrieval algorithm to estimate ET. To this end, we coupled the vegetation evapotranspiration scheme and water evaporation scheme of the P-LSH algorithm with the six existing soil evaporation algorithms (namely, the soil evaporation algorithms of the PML, GLEAM, PT-JPL, PT-Yao, PM-Brust, and P-LSH_ori) to produce six ET retrieval algorithms (i.e., A1–A6 of Table 2 and Fig. 2). Therefore, A1–A6 are comparable to P-LSH_θ and P-LSH_P because the only difference between these algorithms is their soil evaporation component. We then compared the performances of A1–A6, P-LSH_θ, and P-LSH_P for barren areas from January 2003 to August 2011 using ET_recon as the benchmark.

Table 2 Combinations of the six existing soil evaporation algorithms with the P-LSH vegetation evapotranspiration and water evaporation schemes.

Figure 2 Flowchart of the gridded evapotranspiration estimation for a basin used in this study.

The total ET in a pixel is expressed as

$$\begin{array}{}\text{(24)}& E=\sum _i{E}_i{a}_i,\end{array}$$

where i represents the ith land cover in the basin. We ignored land cover that accounted for less than 1 %, so there were grasslands, barrens, and open water for the Qaidam Basin and open shrublands, grasslands, barrens, and open water for the Qiangtang Plateau. The E_i (mm) is the evapotranspiration estimated from the ith land cover, and a_i (–) is the proportion of the ith land cover in a pixel. The open shrubland and open water pixels take the vegetation evapotranspiration scheme and water evaporation scheme from the P-LSH algorithm following K. Zhang et al. (2010, 2015), and the grassland pixels take the vegetation evapotranspiration scheme from the revised P-LSH algorithm following Feng et al. (2022). For barrens, we took the assumption that there was no canopy transpiration, and the performance of the six existing and two improved soil evaporation algorithms were compared. A flowchart of the total ET estimation for the basin is shown in Fig. 2.

We chose the root mean square error (RMSE), coefficient of determination (R²), deviation (bias), and Nash–Sutcliffe efficiency coefficient (NSE) to quantify the performances of remote sensing algorithms compared with the ET_recon as follows:

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{T}}\sum _{i=\mathrm{1}}^T{\left(O_i-S_i\right)}^{\mathrm{2}}},\text{(26)}& {\displaystyle }& {\displaystyle }{R}^{\mathrm{2}}={\left({\displaystyle \frac{\sum _{i=\mathrm{1}}^T\left(O_i-\stackrel{\mathrm{‾}}{O}\right)\left(S_i-\stackrel{\mathrm{‾}}{S}\right)}{\sqrt{\sum _{i=\mathrm{1}}^T{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}\sqrt{\sum _{i=\mathrm{1}}^T{\left(S_i-\stackrel{\mathrm{‾}}{S}\right)}^{\mathrm{2}}}}}\right)}^{\mathrm{2}},\text{(27)}& {\displaystyle }& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{i=\mathrm{1}}^T\left(O_i-S_i\right),\text{(28)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^T{\left(O_i-S_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^T{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}},\end{array}$$

where T is the number of months in the period of record, O is the reconstructed ET, S is the simulated ET, $\stackrel{\mathrm{‾}}{O}$ is the average of all reconstructed values O_i, and $\stackrel{\mathrm{‾}}{S}$ is the average of all simulated S_i.

4 Results

4.1 Performance of existing soil evaporation algorithms

We estimated the daily and $\mathrm{1}/\mathrm{12}$^∘ pixel ET in the Qaidam Basin and Qiangtang Plateau from January 2003 to August 2011 using the six coupling algorithms listed in Table 2. All daily and gridded estimates were aggregated to monthly and basin scales to match ET_recon. Generally, the ET estimates derived from the six coupling algorithms showed large differences. In the Qaidam Basin, the ET estimates of the A1 algorithm (ET_A1) and the A5 algorithm (ET_A5) demonstrated good consistency with the ET_recon, while the ET estimates of the A3 algorithm (ET_A3) and the A6 algorithm (ET_A6) matched the worst. The ET_A1 estimates performed best among all the existing algorithms (Fig. 3a), with an RMSE of 4.06 mm per month, an NSE of 0.88, and an R² of 0.92. The ET estimates of the A2 algorithm (ET_A2) with a linear formula for f were well-simulated for low intervals and were always underestimated for the middle and high intervals (Figs. 3b and 4a). Parameter k in the PT-JPL algorithm was a biome-specific constant and took the same value for all barren pixels, set to 926 Pa, which was calibrated by the ET_recon. Although the parameter k has been calibrated, ET_A3 still could not accurately describe the seasonal variability of ET (Figs. 3c and 4a), mainly because of errors involving f estimates derived by RH and VPD. The medium ET estimates of the A4 algorithm (ET_A4) were always overestimated for the Qaidam Basin (Fig. 3d), which specifically occurred in spring (Fig. 4a). In the PM-Brust method, the biome-specific constants rbl_max and rbl_min for the r_tot estimation were calibrated as 500 and 200 s m⁻¹, respectively, for the Qaidam Basin. The ET_A5 presented good performance (Fig. 3e), with an RMSE of 4.36 mm per month, an NSE of 0.87, and an R² of 0.88. The ET_A6 estimates used RH and VPD to estimate f, with parameter k of 359.1 Pa and r_tot of 462.4 s m⁻¹ following Feng et al. (2022). However, ET_A6 could not adequately describe seasonal variability (Figs. 3f and 4a) in the Qaidam Basin, and seasonal mean values also varied by a large margin compared with ET_recon (Fig. 4a).

Figure 3 Simulated ET derived from the six existing coupling algorithms and the two improved algorithms versus the reconstructed ET (ET_recon). The A1–A6, P-LSH_θ, and P-LSH_P are the coupling algorithms listed in Table 2. The blue dots represent results for Qaidam Basin and the red for Qiangtang Plateau.

Figure 4 Seasonal average reconstructed ET (ET_recon) and ET estimates derived from six existing coupling algorithms and the two improved algorithms for (a) the Qaidam Basin and (b) the Qiangtang Plateau. The A1– A6, P-LSH_θ, and P-LSH_P are the coupling algorithms listed in Table 2.

On the Qiangtang Plateau, almost all algorithms overestimated ET for barren areas in spring and winter and underestimated ET in summer and autumn (Fig. 4b). The multi-year average ET_recon in spring and winter was 6.3 mm per month, while the multi-year average ET derived from six coupling remote sensing algorithms was 14.4±6.8 mm per month. The multi-year average ET_recon in summer and autumn was 38.3 mm per month, and it was 28.2±12.1 mm per month from six remote sensing algorithms. In the comparison of the six algorithms, the ET_A1 estimates still performed best among all algorithms, with an RMSE of 11.14 mm per month, and ET_A2 estimates performed the worst, with an RMSE of 14.46 mm per month. The biome-specific constant, k, in the PT-JPL algorithm was recalibrated to 566 Pa using ET_recon for the Qiangtang Plateau. In the PT-JPL and P-LSH_ori algorithms, unreasonable f estimates also led to the homogenization of strong seasonal variability (Figs. 3c, f, and 4b). Similar to the Qaidam Basin, the ET_A4 estimates showed moderate performance for the Qiangtang Plateau (Fig. 3d), and the ET_A5 estimates showed good performance next to ET_A1, with an RMSE of 11.42 mm per month, an NSE of 0.72, and an R² of 0.85.

Figure 5 The monthly average f derived from the six existing soil evaporation algorithms for the Qaidam Basin and the Qiangtang Plateau.

We calculated the moisture constraint f in the soil evaporation of each pixel and used the same method as ET to aggregate f at the monthly and basin scales. The f estimates derived from various algorithms are shown in Fig. 5. The f estimates of the PML algorithm (f_PML) were high in summer and low in winter, with distinct seasonality in both basins, with small peaks occasionally appearing in winter. The f estimates of the GLEAM algorithm (f_GLEAM) hardly showed seasonality and were always low in both basins, which was the main reason for the poor performance of ET_A2. Compared with f_PML and f_GLEAM estimates, the f estimates of the PT-Yao algorithm (f_PT-Yao) were overestimated in spring and winter, partly causing the overestimation of ET_A4, and this overestimation was larger than that of f_PML (Figs. 3d and 4). Considering the positive relationship between precipitation and soil moisture, the f estimates of the PM-Brust algorithm (f_PM-Brust) from soil moisture and the f_PML estimates from precipitation showed great consistency, with correlation coefficients of 0.86 and 0.85 (p<0.001) for the Qaidam Basin and Qiangtang Plateau, respectively. However, the f_PM-Brust estimates were higher overall than f_PML in spring and winter, and hardly ever close to zero, indicating that the soil moisture sequences over basins seldom reached their minimum at the same time. In addition, compared with f_PML, the overestimation of f_PM-Brust was also a reason for the overestimation of ET_A5 in spring and winter (Figs. 4a and 5). The PT-JPL and P-LSH_ori algorithms shared a similar f estimation and had the same temporal characteristics, with high values in winter and low values in summer, which showed the opposite seasonal variability to soil moisture (expressed in the form of f_PM-Brust). Therefore, the performances of ET_A3 and ET_A6 were unsatisfactory. This is because the VPD sequence for both basins on the TP had stronger seasonality (high in summer and low in winter) compared with the milder RH. Although ET estimates derived from the PT-JPL and P-LSH_ori algorithms have been well validated in some flux towers (Fisher et al., 2008; K. Zhang et al., 2010; Mu et al., 2011), this method is no longer applicable because of the unique meteorology of the TP (mainly manifested in the seasonality of RH and VPD) and the possible decoupling of VPD and soil moisture on a daily scale.

4.2 Performance of the two improved P-LSH algorithms

Because of the good performance of surface soil moisture and precipitation in moisture constraints of the land surface, both were used to improve the P-LSH algorithm, called P-LSH_θ and P-LSH_P. The soil moisture sequence was obtained from the assimilation-based θ_Yang, and the precipitation sequence was obtained from the satellite-based P_GPM. In the original P-LSH algorithm, r_tot was a biome-specific constant sensitive to soil evaporation (Feng et al., 2022). Therefore, we separately calibrated r_tot for both basins in the P-LSH_θ algorithm using the bisection method, with the RMSE as the objective function. The calibrated r_tot values were 575 and 290 s m⁻¹ for the Qaidam Basin and the Qiangtang Plateau, respectively. The ET estimates derived from P-LSH_θ (ET_{P-LSH_θ}) matched well with the ET_recon and captured the strong seasonality of both basins (Fig. 6). The P-LSH_θ algorithm had advantages in normalized standard deviation and centered RMSE, with values of 0.80 and 0.40, while they were 0.61±0.08 and 0.55±0.08 of existing coupling algorithms in Sect. 4.1 (Fig. 7). The r_tot value in the P-LSH_P algorithm for each basin was set the same as that in the P-LSH_θ algorithm. The ET estimates derived from the P-LSH_P (ET_P-LSH¶) were similar to ET_P-LSHtheta and showed a better simulation of the Qaidam Basin, especially the simulations of low values in spring and winter (Fig. 6). However, the ET_{P-LSH_P} estimates were always underestimated on the Qiangtang Plateau, much lower than the ET_recon in summer, which may have been caused by the error of the GPM satellite precipitation on the Qiangtang Plateau (see Sect. 4.3) (Li et al., 2020). Nevertheless, ET_{P-LSH_P} estimates still performed well, second only to ET_{P-LSH_θ} among all algorithms for both basins (Fig. 7). In summary, the P-LSH_θ and P-LSH_P algorithms for both basins showed better performance than the existing algorithms in Sect. 4.1.

Figure 6 Comparisons of the monthly regional average ET estimates derived from two improved retrieval algorithms (P-LSH_θ and P-LSH_P) with the ET_recon for (a) the Qaidam Basin and (b) the Qiangtang Plateau.

Figure 7 Taylor diagram comparing the retrieved ET by the six existing coupling algorithms and the two improved retrieval algorithms (P-LSH_θ and P-LSH_P) in the two basins. The green dashed line represents the centered root mean square error.

The multi-year average annual ET_{P-LSH_θ} and ET_{P-LSH_P} estimates for both basins are shown in Fig. 8. The estimations of the two algorithms shared a similar spatial pattern, with a decreasing trend from the southeastern to northwestern basins. From the perspective of the regional average, ET_{P-LSH_θ} and ET_{P-LSH_P} were 177 and 148 mm for the Qaidam Basin, respectively, and 300 and 232 mm for the Qiangtang Plateau, respectively. However, in the central Qaidam Basin and northwest of the Qiangtang Plateau, the ET_{P-LSH_P} estimates were generally lower than those of ET_{P-LSH_θ}, and these underestimations existed in almost all seasons (Fig. 9). This underestimation was little in winter because both precipitation and soil moisture in winter were low, and the spatial differences between ET_{P-LSH_P} and ET_{P-LSH_θ} almost disappeared.

Figure 8 The spatial distributions of multi-year (January 2003–August 2011) average annual ET derived from (a, b) the P-LSH_θ and (c, d) P-LSH_P for (a, c) the Qaidam Basin and (b, d) the Qiangtang Plateau.

Figure 9 The spatial distributions of multi-year (January 2003–August 2011) seasonal ET derived from (a–d, i–l) P-LSH_θ and (e–h, m–p) P-LSH_P for (a–h) the Qaidam Basin and (i–p) the Qiangtang Plateau.

Figure 9 shows the multi-year spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February) ET_{P-LSH_θ} and ET_{P-LSH_P} in both basins. The pattern of seasonal estimates was similar to that of the annual values. Generally, the ET in autumn was higher than that in spring, with 71 % of ET_{P-LSH_θ} and 97 % of ET_{P-LSH_P} for the Qaidam Basin and 72 % of ET_{P-LSH_θ} and 85 % of ET_{P-LSH_P} for the Qiangtang Plateau (percentage represents the number of pixels accounting for the basin). The multi-year seasonal ET_{P-LSH_θ} and ET_{P-LSH_P} averaged over the Qaidam Basin were 36, 88, 40, and 13 mm, and 20, 87, 33, and 8 mm, for spring, summer, autumn, and winter, respectively. The multi-year seasonal ET_{P-LSH_θ} and ET_{P-LSH_P} averaged over the Qiangtang Plateau were 61, 142, 68, and 29 mm, and 41, 114, 55, and 22 mm, for spring, summer, autumn, and winter, respectively.

Figure 10 Maps of the coefficients of variation of monthly average soil evaporation from 2003 to 2011 for barrens estimated by (a, b) the P-LSH_θ algorithm driven by five soil moisture datasets and (c, d) the P-LSH_P algorithm driven by five precipitation datasets in (a, c) the Qaidam Basin and (b, d) the Qiangtang Plateau. The gray shading indicates the non-barren areas within the basin.

4.3 Uncertainty of soil moisture and precipitation influence on soil evaporation

Two soil moisture constraint schemes based on the respective surface soil moisture and precipitation data were used to improve the P-LSH algorithm (Sect. 4.2), leading to two improved P-LSH algorithms, namely P-LSH_θ and P-LSH_P. Our results show that both of the improved algorithms generally performed well in the two study basins. However, the two improved algorithms are highly dependent on high-quality gridded data, so it is necessary to quantify the influence of uncertainty in the soil moisture and precipitation data on ET estimation. To this end, we selected five surface soil moisture and five precipitation datasets to investigate the impact of uncertainty in the moisture constraint quantification. To be specific, we investigated the difference between ET estimates derived from five surface soil moisture/precipitation datasets, together with the difference among the five soil moisture/precipitation datasets. The daily and $\mathrm{1}/\mathrm{12}$^∘ pixel soil evaporation estimates for both basins were estimated and aggregated to monthly and basin scales. We calculated the coefficient of variation (Cv, defined as the ratio of the standard deviation to the mean) between five barren evaporation estimates from the P-LSH_θ and five barren evaporation estimates from the P-LSH_P algorithms, where the non-barren estimate was masked (hereafter E_{s_P-LSH_θ} and E_{s_P-LSH_P}, where the subscript s denotes soil evaporation for barrens). In the following part, we further discuss the impacts of temporal and spatial uncertainties in soil moisture and precipitation on soil evaporation estimates.

To quantify the temporal uncertainties in soil moisture and precipitation, and their resultant ET uncertainties, we calculated the multi-monthly (i.e., monthly average from 2003 to 2011) values for every soil moisture and precipitation dataset and its associated ET estimate on a grid-cell basis. There was a clear spatial pattern of the Cv among multi-monthly (i.e., monthly average from 2003 to 2011) average E_{s_P-LSH_θ} and E_{s_P-LSH_P} from various datasets. The Cv of E_{s_P-LSH_θ} showed little variation in both basins (Fig. 10a and b). The Cv of E_{s_P-LSH_θ} on the Qaidam Basin ranged from 0.05 to 0.65 with a mean of 0.29, and ranged from 0.03 to 0.71 with a mean value of 0.29 on the Qiangtang Plateau. In contrast, the Cv of E_{s_P-LSH_P} was not consistent in both basins. In the Qaidam Basin, the Cv of E_{s_P-LSH_P} was at a lower level (Fig. 10c), with an average of 0.29, which was comparable to that of E_{s_P-LSH_θ}. On the Qiangtang Plateau, the Cv of E_{s_P-LSH_P} increased from the southeast to the northwest of the basin (Fig. 10d), with an average of 0.46, which was higher than that of E_{s_P-LSH_θ}.

To further distinguish the impact of the datasets and algorithm structure on barren evaporation estimates, we compared the variation for various surface soil moisture and precipitation datasets, as shown in Fig. 11. The surface soil moisture had high uncertainty in the central and northern Qaidam Basin and western Qiangtang Plateau (Fig. 11a and b), but these uncertainties were not reflected in E_{s_P-LSH_θ}, indicating that the moisture constraint calculated by Eq. (20) reduced the uncertainty of soil moisture and, instead, focused more on the relative changes in the soil moisture of each dataset. The Cv of the precipitation showed a similar spatial distribution to that of E_{s_P-LSH_P} in both basins, and their correlation coefficients were 0.93 and 0.98 (p<0.001) for the Qaidam Basin and Qiangtang Plateau, respectively, indicating that the characteristics of precipitation were almost completely transferred to the E_{s_P-LSH_P} through Eq. (21). In contrast, the correlation coefficients of the Cv of soil moisture and E_{s_P-LSH_θ} were only 0.33 and 0.46 (p<0.001) for the two basins.

Figure 11 Maps of the coefficients of variation of (a, b) five soil moisture datasets and (c, d) five precipitation datasets in (a, c) the Qaidam Basin and (b, d) the Qiangtang Plateau. Each soil moisture or precipitation dataset is the monthly average value from 2003 to 2011. The gray indicates the non-barren areas within the basin.

Table 3 RMSE (mm), Bias (mm), NSE, and R² of the five ET_{P-LSH_θ} and five ET_{P-LSH_P} in comparison with the ET_recon for aggregation of two basins.

From the perspective of the regional average, we used Cv to express uncertainty, considering the magnitude between soil moisture and precipitation. There were some uncertainties in various soil moisture datasets, especially in spring, autumn, and winter (Fig. 12a and b), with Cv of 0.41±0.07 and 0.41±0.08 for the Qaidam Basin and Qiangtang Plateau, respectively. The uncertainty of the peaks in summer tended to be much lower and always at a lower level of soil moisture values. The uncertainty among various precipitation datasets was comparable with soil moisture, with Cv of 0.36±0.20 and 0.55±0.22 for the Qaidam Basin and Qiangtang Plateau. The Cv of the precipitation had a similar temporal pattern, similar to soil moisture, low in summer and high in other seasons (Fig. 12c and d). In terms of E_{s_P-LSH_θ} and E_{s_P-LSH_P}, and considering the same object, we used the interval length of various estimates to express uncertainty. Overall, the uncertainty of E_{s_P-LSH_P} was lower than that of E_{s_P-LSH_θ}, especially in spring and winter in both basins (Fig. 12e and f). The interval length of E_{s_P-LSH_P} were 4.94±3.63 and 14.61±10.45 mm per month, and were 11.41±5.91 and 16.92±7.01 mm per month in E_{s_P-LSH_θ} for the Qaidam Basin and Qiangtang Plateau, respectively. On the Qiangtang Plateau, the higher uncertainty of the precipitation datasets led to a larger interval length in the estimation of E_{s_P-LSH_P} compared with the Qaidam Basin, yet this uncertainty was still smaller than that of E_{s_P-LSH_θ}.

Figure 12 Monthly regional average (a, b) soil moisture datasets, (c, d) precipitation datasets, and (e, f) soil evaporation estimates, and their coefficients of variation over barrens of the Qaidam Basin and the Qiangtang Plateau. The shaded areas indicate the range determined by five datasets/estimates and solid lines represent the mean of them, depending on the left y axis. The light dashed lines represent coefficients of variation of five datasets/estimates, depending on the right y axis. Blue represents results of soil moisture datasets or E_s estimates derived from soil moisture, and red represents results of precipitation datasets or E_s estimates derived from precipitation.

From the perspective of simulation accuracy, the ET_{P-LSH_θ} driven by θ_GLDAS Noah performed best, outperforming estimates derived from any other soil moisture and precipitation (Table 3). The satellite-based θ_ESA CCI showed poor performance, which may be attributed to missing data, and simple temporal linear interpolation weakened the seasonal variation in soil moisture. The ET_{P-LSH_θ} estimates derived from θ_Qu and θ_Yang performed well, where soil moisture came from machine learning and model assimilation, respectively. By contrast, the ET_{P-LSH_P} estimates overall had high and relatively stable precision, with an RMSE of 7.70±0.46 mm per month, while it was 8.39±1.08 mm per month from five ET_{P-LSH_θ} estimates.

5 Discussion and conclusion

This paper compared the applicability and effectiveness of various soil moisture constraint schemes in the existing ET algorithms in typical arid and semi-arid basins of the Tibetan Plateau and then proposed two improved P-LSH algorithms, in which normalized surface soil moisture and the ratio of cumulative antecedent precipitation to cumulative antecedent equilibrium evaporation are used to represent soil water stress, respectively. We further assessed the impacts of uncertainty in the soil moisture and precipitation forcing data on the soil evaporation retrievals. The first part of this study investigated the applicability of six existing coupling algorithms with ET_recon in two basins. The moisture constraints and equilibrium equations for these algorithms were different. The A1 algorithm, which considers cumulative precipitation and equilibrium evaporation in soil evaporation, has the best performance on a monthly scale for both basins, with an RMSE of 4.06 mm per month for the Qaidam Basin and RMSE of 11.13 mm per month for the Qiangtang Plateau. The A5 algorithm, which directly considers soil moisture as a constraint, is second in performance, with an RMSE of 4.36 mm per month for the Qaidam Basin and an RMSE of 11.42 mm per month for the Qiangtang Plateau. The ET estimates from the A2 algorithm hardly match well for both basins because they are significantly affected by high-quality soil properties. The A4 algorithm uses the diurnal temperature range to reflect the apparent thermal inertia and humidity constraints, with moderate performance in both basins. Both algorithms, A3 and A6, use an RH${}^{\mathrm{VPD}/k}$ term to express the sensitivity of the soil water deficit, and take the assumption that the surface moisture status is reflected in the adjacent atmospheric moisture, specifically in the form of evaporative demand of the atmosphere. This method has good applicability for ET estimation (Fisher et al., 2008; K. Zhang et al., 2010; Mu et al., 2011), which may be because it pays more attention to total ET rather than soil evaporation. On the barrens of the TP, vegetation is sparse, and only soil evaporation exists; therefore, defects involving this method are exposed. On the TP, RH has weak seasonality, whereas VPD is high in summer and low in winter, with strong seasonal variability. These phenomena result in RH${}^{\mathrm{VPD}/k}$ being high in winter and low in summer, which is contrary to actual soil moisture. In addition, the relationship between VPD and soil moisture may be decoupled on a daily scale (Purdy et al., 2018; Brust et al., 2021), which will eventually lead to model structural errors involving the A3 and A6 algorithms.

The second part of this study improved the P-LSH algorithm by introducing two schemes for quantifying moisture constraints to ET in terms of surface soil moisture and precipitation. From the perspective of the regional average, the two improved algorithms significantly improved the performance of the P-LSH algorithm, and the simulation accuracy was higher than that of the six existing coupling algorithms. The P-LSH_θ algorithm showed the best performance among all algorithms, indicating that direct soil moisture can adequately express the moisture supply in evaporation estimates for barrens. As a surrogate for moisture supply, precipitation can better express the constraints in barrens evaporation than RH, VPD, ATI, etc. However, the two estimates show some uncertainty in the Qiangtang Plateau, which requires more soil evaporation observations or other means to further estimate their reliability.

The last part of this study tested the applicability of satellite soil moisture and precipitation data for improving ET retrieval and analyzing the influence of soil moisture and precipitation uncertainties on ET estimation on alpine barren areas. In the spatial pattern, the uncertainty of E_{s_P-LSH_θ} was lower because the model structure flattened the magnitude difference in soil moisture. On the Qiangtang Plateau, the uncertainty of E_{s_P-LSH_P} is larger, with 47.4 % of the Cv higher than 0.5, which is mainly due to the underestimation of precipitation by the GPM and MERRA in the northwestern basin. From the perspective of the regional average, the uncertainty of soil moisture is comparable to that of precipitation, yet the uncertainty of E_{s_P-LSH_θ} is higher than that of E_{s_P-LSH_P}. The ET_{P-LSH_θ} derived from θ_GLDAS Noah performs better than those from any other soil moisture and precipitation datasets, and the ET_{P-LSH_P} from all precipitation datasets generally showed high and stable accuracy, suggesting that high-quality soil moisture can optimally express moisture supply to ET, and that more accessible precipitation data can serve as a substitute of soil moisture as an indicator of moisture status for its robust performance in barren evaporation.

There were some uncertainties in this study. Because the revisit rates of various satellites are usually 2–3 d, it is difficult to obtain full daily soil moisture coverage of basins, and the satellite-based θ_ESA CCI faces the risk of spatial or temporal discontinuity. Simple temporal linear interpolation was used in our study, which weakened the seasonality of soil moisture. Although differences in various soil moisture datasets were discussed in this study, more spatially and temporally continuous satellite-based soil moisture datasets would be of significant interest. Considering the coarse spatial resolution, uncertainties in the GRACE data are generally much greater; therefore, the ET_recon estimates derived from it also have a coarse temporal and spatial resolution (monthly and basin scale) and high uncertainty. We matched the pixel-scale and daily remote sensing algorithm outputs with the ET_recon, which may cause errors offset in the algorithms to a certain extent. In addition, various processes for GRACE products are sources of uncertainty in ΔS, which in turn affects the accuracy of the ET_recon. Despite the above uncertainties, the water balance method is still an effective means of providing a benchmark for remote sensing algorithm outputs at a basin scale and is recognized in most studies (Zeng et al., 2012; Long et al., 2014; Hui et al., 2020; Chao et al., 2021; Zhang et al., 2020). In terms of results, almost all algorithms had high uncertainty in the simulation of soil evaporation on the Qiangtang Plateau, especially in the summer of 2006 and subsequent years. Zhang et al. (2017) reported that inland lakes on the Qiangtang Plateau have expanded since the 1990s, whereas static land cover was used in this study. In the future, a dynamic dataset will be necessary to reflect the characteristics of the ground surface for ET estimation.