Long-term monthly 0.05° terrestrial evapotranspiration dataset (1982–2018) for the Tibetan Plateau

Evapotranspiration (ET) plays a crucial role in the water balance of the Tibetan Plateau (TP), often referred to as the “Asian water tower” region. However, accurately monitoring and comprehending the spatial and temporal variations of ET components (including soil evaporation E_s, canopy transpiration E_c, and intercepted water evaporation E_w) in this remote area remains a significant challenge due to the limited availability of observational data. This study generates a 37-year dataset (1982–2018) of monthly ET components for the TP using the MOD16-STM (MOD16 soil texture model). This model utilizes up-to-date soil properties, meteorological data, and remote sensing datasets. The estimated ET results strongly correlate with measurements from nine flux towers, demonstrating a low root mean square error of 13.48 mm per month, a mean bias of 2.85 mm per month, a coefficient of determination of 0.83, and an index of agreement of 0.92. The annual average ET for the entire TP, defined as elevations higher than 2500 m, is approximately 0.93(± 0.037) × 10³ Gt yr⁻¹. The predominant contributor to ET on the TP is E_s, accounting for 84 % of the total ET. Our findings reveal a noteworthy upward trend in ET in most central and eastern parts of the TP, with a rate of approximately 1–4 mm yr⁻¹ (p<0.05) and a significant downward trend with rates between −3 and 1 mm yr⁻¹ in the northwestern part of TP during the period from 1982 to 2018. The average annual increase in ET for the entire TP over the past 37 years is approximately 0.96 mm yr⁻¹. This upward trend can be attributed to the TP's warming and wetting climate conditions. The MOD16-STM ET dataset demonstrates a reliable performance across the TP compared to previous research outcomes. This dataset is valuable for research on water resource management, drought monitoring, and ecological studies. The entire dataset is freely accessible through the Science Data Bank (https://doi.org/10.11922/sciencedb.00020, Ma et al., 2021) and the National Tibetan Plateau Data Center (TPDC) (https://doi.org/10.11888/Terre.tpdc.271913, Yuan, 2021).

1 Introduction

The Tibetan Plateau (TP) (24–40^∘ N, 70–105^∘ E) is often referred to as the “Asian water tower”, owing to its distinctive geographical and ecological characteristics, as acknowledged in studies by Immerzeel et al. (2010, 2020), Yao et al. (2012), and Xu et al. (2019). Within this region, evapotranspiration (ET) plays a vital role in the overall water balance. The TP predominantly features grassland (covering more than 47 % of the area) and sparse vegetation or bare soil (surrounding over 33 %), as indicated by the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover dataset (MCD12C1) (Fig. 1c). Arid or semi-arid conditions mostly characterize this vast expanse. The TP is currently undergoing significant changes in its hydrological cycle, driven by global warming, as documented in studies by Yang et al. (2014), Kuang and Jiao (2016), and Zohaib et al. (2017). Nevertheless, accurately monitoring the spatial and temporal fluctuations in ET on the TP remains a formidable challenge due to the intricate environmental conditions of the TP. Moreover, understanding how ET patterns on the TP will evolve in the context of global warming is essential for assessing the impacts of these changes on the local population's livelihoods.

In recent years, various datasets for estimating ET on the TP have been developed, including the complementary relationship (CR) model (Ma et al., 2019; G. X. Wang et al., 2020), the surface energy balance system (SEBS) model (Chen et al., 2014, 2021; Zhong et al., 2019; Han et al., 2017, 2021), and the Penman–Monteith model with remote sensing (RS-PM) (Wang et al., 2018; Song et al., 2017; Chang et al., 2018; Ma and Zhang, 2022). However, a considerable variance exists among these TP ET products (Peng et al., 2016; Baik et al., 2018; Li et al., 2018; Khan et al., 2018). Studies have utilized eddy-covariance measurements (Shi and Liang, 2014; You et al., 2017; Yang et al., 2019; Ma et al., 2020) and reanalysis datasets (Shi and Liang, 2014; Dan et al., 2017; Yang et al., 2019; de Kok et al., 2020) to investigate ET on the TP. A recent study by Han et al. (2021) produced the region's 18-year ET dataset (2001–2018). Enhancements of the canopy conduction algorithm in the Penman–Monteith model have led to improved ET estimates in previous research (Leuning et al., 2010; Zhang et al., 2010; Li et al., 2015; Y. Zhang et al., 2016, 2019; Gan et al., 2018). However, these ET products tend to perform poorly in TP areas with sparse vegetation and arid to semi-arid climates (Zhang et al., 2010; Li et al., 2014b; Song et al., 2017; Baik et al., 2018; Li et al., 2018; Khan et al., 2018).

The limitations of the MOD16 Penman–Monteith model in arid to semi-arid TP regions are primarily due to its failure to consider the dominant role of topsoil texture and topsoil moisture in governing E_s processes (Yuan et al., 2021). Accurately separating and validating ET components on the TP remains challenging, even though total ET estimates tend to align across different products (Lawrence et al., 2007; Blyth and Harding, 2011; Miralles et al., 2016). The TP is primarily characterized by short and sparse vegetation, and soil moisture is crucial in ET estimation for this region. Several studies have used the Penman–Monteith algorithm to estimate ET on the TP (Wang et al., 2018; Ma and Zhang, 2022). However, these studies have not accounted for the effects of soil moisture (SM) on evaporation resistance and stomatal conductance.

The enhanced Penman–Monteith model, MOD16-STM (MOD16 soil texture model), has been developed to address these limitations. MOD16-STM redefines the modules for E_s to consider the impacts of SM on soil evaporation resistance. This modification is based on eddy-covariance (EC) observations conducted on the TP (Yuan et al., 2021), offering a promising opportunity to estimate ET components in this region accurately. E_s often dominates ET in sparsely vegetated areas, especially in arid and semi-arid regions with large bare soil areas (Wilcox et al., 2003; Kool et al., 2014; Wang et al., 2018; Ma et al., 2015a, b; Ma and Zhang, 2022). Previous studies have highlighted that 20 % to 40 % of global ET is attributed to E_s (Lawrence et al., 2007; Schlesinger and Jasechko, 2014; Wei et al., 2017). Bare soil surface evaporation is a rapid process influenced by shallow surface water (Koster and Suarez, 1996). E_s is primarily controlled by water diffusion in the soil (Good et al., 2015; Yuan et al., 2021). Accurately quantifying and separating E_s is crucial for enhancing our understanding of water and energy cycles on the TP. However, quantifying ET and its components remains challenging due to the influence of atmospheric demand, soil moisture conditions, and complex interactions between heterogeneous vegetation and soil properties (Merlin et al., 2016; Wu et al., 2017; Philips et al., 2017; Lehmann et al., 2018). MOD16-STM holds the potential to generate a remote sensing E_s and ET component dataset covering the satellite era since 1980. In this study, the MOD16-STM model, acknowledging its limitations, was employed to estimate a long-term ET and ET components dataset spanning 37 years (1982–2018) (Yuan et al., 2021).

A preferable approach involves directly estimating ET based on topsoil moisture, significantly impacting the TP's surface water exchange. Thus, leveraging the advantages of the MOD16-STM model for ET estimation on the TP, this study aimed to achieve two main objectives: (1) develop a 37-year (1982–2018) monthly ET dataset for the TP at a 0.05^∘ × 0.05^∘ spatial resolution and (2) quantify the spatial distribution and spatiotemporal variability of ET and its components across the TP.

2 Materials and methods

2.1 Study area

The Tibetan Plateau, located between 25–40^∘ N and 74–104^∘ E, spans approximately 2.5×10⁶ km² and consists of land above 2500 m in altitude (Fig. 1a). This region, as indicated by the FAO drought index dataset, represents the largest landform unit in Eurasia and encompasses hyper-arid, arid, semi-arid, and sub-humid climate zones (Fig. 1b). The land cover types primarily include mixed forests, grasslands, bare soil, glaciers, and snow-covered areas (see Fig. 1c). The topsoil predominantly consists of sandy loam, loam, and clay (Fig. 1d). The annual average temperature in the region ranges from approximately −3.1 to 4.4 ^∘C. Average annual precipitation gradually increases from less than 50 mm in the northwest to over 1000 mm in the southeast, with the most precipitation occurring during summer (Ding et al., 2017). Over time, the TP has undergone significant environmental changes (Chen et al., 2015), including increased precipitation, decreased wind speed (wind), fewer snow days, reduced radiation, thawing permafrost, glacier melting, and increased vegetation (Kang et al., 2010; Yao et al., 2012; Yang et al., 2014; Kuang and Jiao, 2016; Bibi et al., 2018; Chen et al., 2019).

Figure 1 Maps of the (a) topography (STRM), (b) climate zones (FAO aridity index), (c) land cover types (MCD12C1), and (d) soil textures (HWSD) in the study area. The red dots indicate the flux site locations.

2.2 Generation of a long-term series of monthly ET products

This study introduces a novel dataset comprising a long-term series of monthly ET generated using the MOD16-STM model. The process of calculating monthly ET with the MOD16-STM model and the associated driving datasets is illustrated in Fig. 2.

Figure 2 Workflow of the MOD16-STM evapotranspiration product.

2.2.1 Description of MOD16-STM ET model

The MOD16-STM model computes the ET components using the Penman–Monteith equation as follows (Monteith, 1965; Mu et al., 2007, 2011):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& E_{\mathrm{c}}=\\ & {\displaystyle \frac{\left(\mathrm{\Delta }\times f_{\mathrm{c}}\times (R_{\mathrm{n}}-G_{\mathrm{0}})+{\mathit{\rho }}_{\mathrm{a}}\times C_p\times \frac{\mathrm{VPD}}{r_{\mathrm{a}}}\times f_{\mathrm{c}}\right)\times (\mathrm{1}-F_{\mathrm{wet}})}{\mathit{\lambda }\times \left(\mathrm{\Delta }+\mathit{\gamma }\times (\mathrm{1}+\frac{r_{\mathrm{c}}}{r_{\mathrm{a}}})\right)}}\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& E_{\mathrm{s}}=\\ & {\displaystyle \frac{\left(\mathrm{\Delta }\times (\mathrm{1}-f_{\mathrm{c}})\times (R_{\mathrm{n}}-G_{\mathrm{0}})+{\mathit{\rho }}_{\mathrm{a}}\times C_p\times \frac{\mathrm{VPD}}{r_{\mathrm{a}}}\right)\times (\mathrm{1}-F_{\mathrm{wet}})}{\mathit{\lambda }\times \left(\mathrm{\Delta }+\mathit{\gamma }\times (\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}})\right)}}\\ & \times {\left({\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}}\right)}^{\frac{\mathrm{VPD}}{\mathit{\beta }}}\end{array}\text{(3)}& {\displaystyle }{E}_{\mathrm{w}}=E_{\mathrm{wet}\mathrm{\_}\mathrm{s}}+E_{\mathrm{wet}\mathrm{\_}\mathrm{c}}.\end{array}$$

The total ET combines three distinct components: E_c,E_s, and E_w (wet surface evaporation). For a more detailed explanation of the calculations for E_{wet_s} (evaporation from wet soil) and E_{wet_c} (evaporation from wet canopy), you can refer to Yuan et al. (2021). Here, r_a (s m⁻¹) is the aerodynamic resistance, r_c (s m⁻¹) is the aerodynamic resistance of water vapor of the canopy, and r_s (s m⁻¹) is the surface (or canopy) resistance. Yuan et al. (2021) optimized MOD16 r_a based on the Monin–Obukhov similarity theory (MOST) and calibrated the empirical values of r_c for grassland underlying surfaces. They also pointed out that the topsoil moisture content directly affects the value of r_s, indirectly influencing the E_s process. Therefore, this study extended this optimization algorithm from the site scale to the regional scale. The variables used in the above equations are defined as follows:

• R_n represents the net radiation flux (W m⁻²).

• G₀ denotes the soil heat flux (W m⁻²).

• ρ_a is the density of the air (kg m⁻³).

• C_p stands for the specific heat capacity of the air (J kg⁻¹ K⁻¹).

• VPD represents the vapor pressure deficit (hPa).

• Δ represents the slope of the saturated vapor pressure curve (hPa K⁻¹).

• γ is the psychrometric constant (hPa K⁻¹), calculated as $\mathit{\gamma }=C_p\cdot P_{\mathrm{a}}\cdot M_{\mathrm{a}}/(\mathit{\lambda }\cdot M_{\mathrm{w}})$, where λ is the latent heat of vaporization (J kg⁻¹), and M_a and M_w are the molecular masses of dry air and wet air, respectively.

• r_a signifies the aerodynamic resistance (s m⁻¹).

• r_s represents the surface (or canopy) resistance (s m⁻¹).

• F_wet is the relative surface wetness.

• The vegetation cover fraction (f_c) is estimated using the NDVI (Normalized Difference Vegetation Index):

  $$\begin{array}{}\text{(4)}& f_{\mathrm{c}}={\left({\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{min}}{{\mathrm{NDVI}}_{max}+{\mathrm{NDVI}}_{min}}}\right)}^{\mathrm{2}}.\end{array}$$

R_n and G₀ are calculated as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{R}_{\mathrm{n}}=(\mathrm{1}-\mathit{\alpha })\times \mathrm{SWD}+\mathrm{LWD}-\mathit{\epsilon }\times \mathit{\sigma }\times {\mathrm{LST}}^{\mathrm{4}}\text{(6)}& {\displaystyle }{G}_{\mathrm{0}}=R_{\mathrm{n}}\times (I_{\mathrm{c}}+(\mathrm{1}-f_{\mathrm{c}})\times (I_{\mathrm{s}}-I_{\mathrm{c}}\left)\right).\end{array}$$

Here, SWD is the downward shortwave radiation, α is land surface albedo, LWD is the downward longwave radiation, σ represents the Stefan–Boltzmann constant (5.67×10⁻⁸ W m⁻² K⁻⁴), ε is emissivity, and LST is the land surface temperature. I_c (= 0.05) and I_s (= 0.315) are the ratios of ground heat flux and net radiation for surfaces with full vegetation cover (Su, 2002) and bare soil (determined by NDVI < 0.25 in this study) (Yuan et al., 2021), respectively. When the air temperature (T_a) is below 5 ^∘C, photosynthesis and transpiration processes are not active, and, therefore, E_c is not considered in the calculations. When the land surface temperature is below 0 ^∘C, the sublimation equation is derived by modifying the surface energy balance equation using the Clausius–Clapeyron equation, accounting for the equilibrium of water vapor in both liquid and frozen states. It is important to note that this study did not estimate the evaporation from water surfaces. Previous research has extensively examined water surface evaporation from lakes on the Tibetan Plateau in detail (Wang et al., 2020). Therefore, this study focuses on land ET estimation, excluding water surface evaporation.

Numerous prior studies have employed optimized conductance to estimate E_c (Jarvis, 1976; Irmak and Mutiibwa, 2010; Zhang et al., 2010; Leuning et al., 2010; Li et al., 2013, 2015), as well as E_s (Sun, 1982; Camillo and Gurney, 1986; Sellers et al., 1996; Sakaguchi and Zeng, 2009; Ortega-Farias et al., 2010; Tang and Riley, 2013). This study computed the r_a using the MOST (Thom, 1972; Liu et al., 2007).

$$\begin{array}{}\text{(7)}& r_{\mathrm{a}}={\displaystyle \frac{\ln \left(\frac{z_h-d_{\mathrm{0}}}{z_{\mathrm{0}h}}-{\mathrm{\Psi }}_h\right)\ln \left(\frac{z_m-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{m}}}-{\mathrm{\Psi }}_m\right)}{k^{\mathrm{2}}u}},\end{array}$$

where k represents von Karman's constant (0.41), z_h and z_m denote the measurement heights for T_a and wind, and d₀ represents the displacement height. The stability correction functions for momentum (ψ_m) and heat transfer (ψ_h) can be computed using universal parts. These correction terms' mathematical expressions (Eqs. 8–12) are as follows (Högström, 1996; Paulson, 1970).

For stable conditions,

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{m}}=-\mathrm{5.3}{\displaystyle \frac{(z_m-z_{\mathrm{0}\mathrm{m}})}{L}}\text{(9)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{h}}=-\mathrm{8.0}{\displaystyle \frac{z_h-z_{\mathrm{0}h}}{L}}.\end{array}$$

For unstable conditions,

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{\mathit{\psi }}_{\mathrm{m}}& =\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}+x}{\mathrm{1}+x_o}}\right)+\ln \left({\displaystyle \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{1}+x_o^{\mathrm{2}}}}\right)-\mathrm{2}{\tan }^{-\mathrm{1}}x\\ & +\mathrm{2}{\tan }^{-\mathrm{1}}{x}_o\end{array}\text{(11)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{h}}=\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}+y}{\mathrm{1}+y_o}}\right).\end{array}$$

For neutral conditions,

$$\begin{array}{}\text{(12)}& {\mathit{\psi }}_{\mathrm{m}}={\mathit{\psi }}_{\mathrm{h}}=\mathrm{0}.\end{array}$$

In Eqs. (8)–(11), the following variables and parameters are defined: $x=(\mathrm{1}-z_m/L{)}^{\mathrm{0.25}}$, $x_o=(\mathrm{1}-z_{om}/L{)}^{\mathrm{0.25}}$, y= (1–$\mathrm{11.6}{z}_h/L$)^0.5, and y_o= (1–11.6 $z_{oh}/L$)^0.5 (Högström, 1996; Paulson, 1970). Here, L represents the Obukhov length (m), calculated as $L=T_{\mathrm{a}}\cdot u_{\ast }^{\mathrm{2}}/$(kg T_∗), where g=9.8 m s⁻², T_∗ is the fractional temperature (K), and u_∗ denotes the friction velocity (m s⁻¹). T_∗ is further defined as $T_{\ast }=-({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{a}})/\left(\right(\ln (z_h/z_{oh})-{\mathit{\psi }}_{\mathrm{h}})$, where θ_s can be approximated using the LST, and ${\mathit{\theta }}_{\mathrm{a}}=T_{\mathrm{a}}+z_h\cdot g/C_p$. The parameterization of u_∗ and L has been successfully applied in previous studies on the TP (Chen et al., 2013). z_0h represents the roughness length for heat transfer (m). A parameterization scheme for z_0h developed by Yang et al. (2008) has been widely utilized in remote sensing land surface fluxes and land surface models (LSMs) across the TP (Biermann et al., 2014; Chen et al., 2013; Ma et al., 2015a). This scheme has also been employed in the current study for consistency.

$$\begin{array}{}\text{(13)}& z_{\mathrm{0}h}={\displaystyle \frac{\mathrm{70}v}{u_{\ast }}}\exp \left(-\mathrm{7.2}{u}_{\ast }^{\mathrm{0.5}}{\left|T_{\ast }\right|}^{\mathrm{0.25}}\right),\end{array}$$

where v is the fluid kinematic viscosity, v=1.328 × 10${}^{-\mathrm{5}}\cdot (P_{\mathrm{0}}/P_{\mathrm{a}})\cdot (T_{\mathrm{a}}/T_{\mathrm{0}}{)}^{\mathrm{1.754}}$, P₀=1013 hPa, and T₀=273.15 K. The roughness height for momentum transfer (z_0m) was determined based on canopy height (h_c), following the method outlined by Chen et al. (2013). The water saturation degree of surface soil (SM$/{\mathit{\theta }}_{\mathrm{sat}}$) is utilized to impose soil classification and soil texture constraints on the r_s and E_s estimates (Yuan et al., 2021), as follows:

$$\begin{array}{}\text{(14)}& r_{\mathrm{s}}=\exp \left(a+b\times {\displaystyle \frac{\mathrm{SM}}{{\mathit{\theta }}_{\mathrm{sat}}}}\right).\end{array}$$

Here, the parameters a and b are empirical coefficients that vary based on different soil textures, as documented in Table 1. The estimation of θ_sat, which considers soil organic content (SOC) and gravel content, can be obtained using the Soc-Vg scheme (Chen et al., 2012; Zhao et al., 2018), and its calculation is as follows:

$$\begin{array}{}\text{(15)}& {\mathit{\theta }}_{\mathrm{sat}}=(\mathrm{1}-V_{\mathrm{SOC}}-V_g)\times {\mathit{\theta }}_{\mathrm{sat},\mathrm{m}}+V_{\mathrm{SOC}}\times {\mathit{\theta }}_{\mathrm{sat},\mathrm{sc}},\end{array}$$

where θ_sat,m represents the porosity of the mineral soil and can be calculated as ${\mathit{\theta }}_{\mathrm{sat},\mathrm{m}}=\mathrm{0}$.489–0.00126⋅ %sand (Cosby et al., 1984). Additionally, θ_sat,sc is the porosity of the SOC. It is assumed to be 0.9 m³ m⁻³ in this study, as per the work of Farouki (1981) and Letts et al. (2000). The variables V_soc and V_g denote the volumetric fractions of the SOC and gravel, respectively, and their calculation is as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}& V_{\mathrm{SOC}}=\\ & {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}\times (\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\times m_{\mathrm{SOC}}}{\begin{array}{c}{\mathit{\rho }}_{\mathrm{SOC}}\times (\mathrm{1}-m_{\mathrm{SOC}})+{\mathit{\rho }}_{\mathrm{p}}\times (\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\\ \times m_{\mathrm{SOC}}+(\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\times \frac{{\mathit{\rho }}_{\mathrm{SOC}}\times m_{\mathrm{g}}}{\mathrm{1}-m_{\mathrm{g}}}\end{array}}}\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& V_{\mathrm{g}}=\\ & {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{SOC}}\times (\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\times m_{\mathrm{g}}}{\begin{array}{c}(\mathrm{1}-m_{\mathrm{g}})\times \left({\mathit{\rho }}_{\mathrm{SOC}}\times (\mathrm{1}-m_{\mathrm{SOC}})+{\mathit{\rho }}_{\mathrm{p}}\times (\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\right.\\ \left.\times m_{\mathrm{SOC}}+(\mathrm{1}-{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}})\times \frac{{\mathit{\rho }}_{\mathrm{SOC}}\times m_{\mathrm{g}}}{\mathrm{1}-m_{\mathrm{g}}}\right)\end{array}}}.\end{array}\end{array}$$

In these equations, ρ_p represents the mineral particle density and is set at 2700 kg m⁻³, while ρ_soc is the bulk density of organic matter, maintained at 130 kg m⁻³. Also, m_soc and m_g denote the percentages of organic matter and gravel within the topsoil layer.

There are many parameters in Eqs. (7)–(17). Some parameters have already been assessed for their importance in ET estimation by Yuan et al. (2021). There are too many studies on investigating the empirical parameters in Eqs. (7)–(13). We will not repeat these analyses again. The parameterization method of θ_sat in the estimation of r_s in this study is composed of various empirical parameters (ρ_p, ρ_soc, and θ_sat,sc) for different soil types. We have conducted a uncertainty analysis of the estimated θ_sat and sensitivity of its uncertainty to the changes of empirical parameters. The impact of the empirical parameters on the estimation of ET is illustrated in Fig. A1a–c. The results indicate that with a 20 % uncertainty range in the estimated parameters ρ_p, ρ_soc, and θ_sat,sc for θ_sat, the loss in estimating ET is only below 3 %. Thus, the conclusion is drawn that the estimation of ET is not sensitive to uncertainty in the three parameters. Figure A2 also shows the accuracy of the estimated θ_sat by the method used in this paper. Additionally, a sensitivity analysis is conducted on the empirical parameters a and b for calculating r_s. Keeping θ_sat and SM constant, r_s exhibits exponential changes with variations in a and b, leading to significant fluctuations in the estimation of ET. Within a 20 % range of variation in a and b, the maximum loss in ET exceeded 50 % in Fig. A1d–e. Therefore, it is essential to perform significance tests on the fitting results of the empirical parameters a and b, as well as independent validation of the final ET estimates. The performance of soil surface resistance r_s estimated by the MOD16-STM model at the site scale is demonstrated in Fig. 3. The observations at the 10 stations show that the soil surface resistance exponentially decreases with the increasing SM$/{\mathit{\theta }}_{\mathrm{sat}}$. The MOD16-STM has caught this exponential law. It has a coefficient of determination (R²) higher than 0.34, which may enable the model to estimate the TP ET reasonably.

This demonstrates that the parameterized method of θ_sat maintains a high level of consistency with the observed values. The sensitivity test show that the factors that have a significant impact on r_s and ET are the topsoil moisture and soil organic matter content. Figure A3 present the impact of soil organic matter content on θ_sat and ET estimation at different soil types. Hereby, we have collected the most updated soil texture and soil moisture data to estimate the soil evaporation resistance.

Table 1 Parameters a and b were based on different soil textures used to calculate surface soil resistances.

Figure 3 The scatter point relationship between soil surface resistance (r_s) and SM$/{\mathit{\theta }}_{\mathrm{sat}}$ observed at QOMS (sandy loam), NAMOR (sandy), NASDE (loamy sand), Arou (silt loam), Maqu (loam), IT-CAS (clay loam), US-IHO (sandy clay loam), US-Arm (clay), CH-Oe2 (silty clay), and US-Ib2 (silty clay loam). The red curves show the equations used by the MOD16-STM model. The site information is given in Tables 3 and A1.

2.2.2 Input data for calculating the TP ET

The MOD16-STM model relies on various remote sensing datasets, reanalysis datasets, and meteorological forcing datasets to estimate monthly ET across the TP. Specific datasets are carefully selected to minimize spatial and temporal gaps in the final product (Table 2). Here is a breakdown of the critical datasets and their sources:

• Monthly meteorological forcing data, including wind, T_a, air specific humidity (q), P_a, SWD, and LWD, were obtained from the China Meteorological Forcing Dataset (CMFD) with a 0.1^∘ spatial resolution from 1982–2018. This data source was accessed from the National Tibetan Plateau Data Center (Yang et al., 2010; He et al., 2020). CMFD can be downloaded from TPDC (https://data.tpdc.ac.cn/, last access: 1 February 2024).

• LST and precipitation data were sourced from ERA5-Land, which provides data with a spatial resolution of 0.1^∘ and a monthly temporal resolution. These datasets were obtained from the European Centre for Medium-Range Weather Forecasts (ECWMF).

• Albedo (α) data, with a spatial resolution of 0.05^∘ and an 8 d temporal resolution, were derived from the Global Land Surface Satellite (GLASS) dataset (Liang et al., 2021).

• A long-term NDVI dataset, with a spatial resolution of 0.05^∘ and daily temporal resolution, was acquired from the National Oceanic and Atmospheric Administration's National Centers for Environmental Information (NOAA-NCEI) (https://www.ncei.noaa.gov/products/climate-data-records/normalized-difference-vegetation-index, last access: 1 February 2024). This dataset calculates the canopy height and leaf area index (LAI) (Chen et al., 2013).

• A topsoil moisture dataset for the 0–10 cm depth, with a spatial resolution of 0.25^∘ and a monthly temporal resolution, was obtained from the Global Land Evaporation Amsterdam Model (GLEAM) (Miralles et al., 2011).

• Upward longwave radiation (LWU) was derived from LST using the Stefan–Boltzmann law. The emissivity of mixed pixels was calculated based on the specific emissivity values for vegetated and bare land surfaces, following Sobrino et al. (2004).

• Soil texture and soil property data were obtained from the Harmonized World Soil Database v1.2 (HWSD) (Wieder et al., 2014). These data were used to calculate soil evaporation resistance.

Daily and 8 d input data were averaged over the temporal scale to create monthly datasets to ensure consistency. The average value was considered invalid if the ratio of valid data in any given month was below 90 %. Additionally, the spatial resolutions of all input datasets were interpolated to a standard 0.05^∘ spatial resolution using a widely used bilinear interpolation method.

Table 2 Input datasets are used to calculate the ET on the Tibetan Plateau.

2.3 Validation methods

2.3.1 Model validation at site scale

Limited stations on the TP make it impossible to collect ET observations at all kinds of soil textures. We have collected datasets from 17 flux sites outside the TP (Table A1 in Appendix A). Five sites are used to verify the relationship between soil surface resistance and SM$/{\mathit{\theta }}_{\mathrm{sat}}$. This result is presented in Fig. 3. The other 12 verification sites include 10 different soil textures (sandy loam, sand, loamy sand, silt loam, loam, clay loam, sandy clay loam, clay, silty clay, and silty clay loam) and three surface cover types (grassland, evergreen forest, and cropland) (Table A1). These 12 sites are used to do model validation at the site scale. When evaluating the MOD16-STM at the site scale, the meteorological forcing data come from the station measurement. This helps us to minimize the simulation uncertainty due to the errors in the model forcing datasets. This methodology can allow us to diagnose the model's limitation in representing the evapotranspiration process. Figure A4 shows that MOD16-STM can capture the ET variations at the 12 sites. Table A2 also lists the statistical values of the daily ET estimation. Since these sites include all kinds of soil textures and different canopy covers, we believe the MOD16-STM model could be applied to the TP regional scale.

2.3.2 ET product evaluation

The remote sensing ET product is validated through comparison with flux tower observations on the TP. Table 3 lists details of nine flux stations on the TP used for the ET product evaluation. These stations belong to ChinaFLUX (Dang-Xiong site (DX), Hai-Bei site (HB), Yu et al., 2006; L. M. Zhang et al., 2019), the Tibetan Observation and Research Platform (TORP) (BJ, NADORS, SETORS, QOMS, NAMORS, and Shuang-Hu (SH), Ma et al., 2020), and the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) (Arou, Liu et al., 2011, 2018; Che et al., 2019) networks. The nine stations are in areas with three different land cover types: alpine meadow, alpine steppe, and Gobi. Half-hourly flux data measured by eddy covariance from the nine stations are collected. It is important to note that the energy balance closure ratio (ECR) indicates whether the sum of sensible heat (H), latent heat (LE), and soil heat flux (G₀) matches the R_n. Half-hourly data are screened and corrected accordingly to ensure the reliability of eddy-covariance measurements. Half-hourly LE data are corrected using the Bowen ratio energy balance correction method (Chen et al., 2014).

$$\begin{array}{}\text{(18)}& {\displaystyle }\mathrm{ECR}={\displaystyle \frac{H+LE}{R_{\mathrm{n}}-G_{\mathrm{0}}}}\text{(19)}& {\displaystyle }L{E}_{\mathrm{cor}}={\displaystyle \frac{\mathrm{1}}{\mathrm{ECR}}}\times LE\end{array}$$

Table 3 Details of the nine flux observation stations on the TP used for the ET product evaluation.

In the validation process, the half-hourly LE_cor data obtained from all the flux sites are subjected to further processing, including conversion to daily and monthly averages, while employing a stringent quality control procedure. Daily values are null if derived from valid data points amounting to less than 80 % in a single day. Similarly, monthly average values are disregarded in the validation if they are derived from valid data points accounting for less than 80 % of observations for that month. This approach ensured the robustness of the validation process.

2.4 Accuracy metrics

The accuracy of the modeled ET was assessed by comparing the pixel values (M_i), corresponding to the latitude and longitude of the flux site, with the flux tower measurements (G_i). Several statistical metrics are employed for validation, including the coefficient of determination (R²), a measure of the proportion of the variance in the observed data (G_i) that is explained by the modeled data (M_i). A higher R² value indicates a stronger linear relationship between the two datasets. The mean bias (MB) represents the average difference between the modeled ET (M_i) and the observed flux tower measurements (G_i). Positive MB values suggest overestimation by the model, while negative values indicate underestimation. The root mean square error (RMSE) measures the standard deviation of the differences between modeled and observed values (M_i−G_i). A smaller RMSE implies greater accuracy in the model's predictions. The index of agreement (IOA) indicates the degree of agreement between modeled and observed data, with a value of 1 indicating perfect agreement. Higher IOA values indicate better agreement between the two datasets. The equations for these parameters are as follows:

$$\begin{array}{}\text{(20)}& {\displaystyle }{R}^{\mathrm{2}}={\displaystyle \frac{{\left({\sum }_{i=\mathrm{1}}^n(M_i-\stackrel{\mathrm{‾}}{M})(G_i-\stackrel{\mathrm{‾}}{G})\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(M_i-\stackrel{\mathrm{‾}}{M}{)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n(G_i-\stackrel{\mathrm{‾}}{G}{)}^{\mathrm{2}}}},\mathrm{0}\le R^{\mathrm{2}}\le \mathrm{1}\text{(21)}& {\displaystyle }\mathrm{MB}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^n(M_i-G_i)\text{(22)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{n}}}\sum _{i=\mathrm{1}}^n(M_i-G_i{)}^{\mathrm{2}}}\text{(23)}& {\displaystyle }\mathrm{IOA}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(M_i-G_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(\mathrm{|}{M}_i-\stackrel{\mathrm{‾}}{G}\mathrm{|}+\mathrm{|}{G}_i-\stackrel{\mathrm{‾}}{G}\mathrm{|}{)}^{\mathrm{2}}}}.\end{array}$$

The subscript i denotes individual samples, and n is the total number of samples used in the assessment. The significance of each parameter helps evaluate the model's performance in estimating ET.

3 Results

3.1 Evaluation of ET products against flux tower measurements

The reliability of remote-sensing-based ET estimates is often questioned without ground measurement verification. This study compares the simulated monthly ET rates from the 0.05^∘ grid where each EC site is located with the flux tower observational data to validate the MOD16-STM ET results. The validation outcomes for monthly MOD16-STM ET, using flux tower data, are illustrated in Fig. 4. The modeled ET exhibits excellent performance and high consistency across the TP compared to ET observations.

Specifically, the grassland sites (SETORS, Arou, DX, and HB) display strong agreement, with R² and IOA values exceeding 0.82 and 0.95, respectively. The NAMORS site has the lowest performance, with the highest RMSE (17.8 mm per month) and the lowest R² and IOA (0.63 and 0.87, respectively). On average, the mean R² and IOA values exceed 0.83 and 0.93, respectively. All R² values pass the significance test at the p<0.05 level. The mean values of the absolute bias (| MB |) and RMSE values are less than 3 and 14 mm per month. It is important to note that positive MB values indicate an overestimation of ET, particularly during the dry season over barren land (QOMS, DX, SH, and NADORS) (Fig. 4). Conversely, underestimation occurs at higher ET rates in the summer, likely because the soil is close to saturation, leading to an overestimation of r_sand underestimation of E_s and ET. Generally, the time series of ET variations at the nine flux tower stations exhibit seasonal and annual periodic variations (Fig. 5). The site-scale validation results demonstrate that MOD16-STM ET provides accurate estimates in the TP region.

Figure 4 The validation of the MOD16-STM monthly ET at (a) SETORS, (b) Arou, (c) HB, (d) QOMS, (e) DX, (f) NAMORS, (g) BJ, (h) SH, (i) NADORS, and (j) all sites.

Figure 5 Time-series variations in the MOD16-STM simulated ET (solid blue line with “*” marks) and flux-tower-observed ET (red circles) at (a) SETORS, (b) Arou, (c) HB, (d) QOMS, (e) DX, (f) NAMORS, (g) BJ, (h) SH, and (i) NADORS.

3.2 Spatial pattern of the multiyear averaged ET across TP

Figure 6 displays the spatial distribution of the average annual ET and its three components across the TP. ET exhibits a decreasing trend from southeast to northwest, with the highest values exceeding 1000 mm yr⁻¹ in the southeastern TP (the Hengduan Mountains) and the lowest values of less than 100 mm yr⁻¹ in the Qaidam Basin and northwestern TP. This spatial pattern of annual ET closely mirrored that of the aridity index (Fig. 1b), which is influenced by atmospheric demand and water supply. The sub-humid zone, covering approximately 32.9 % of the TP, contributes the highest proportion (43 % of the TP's total ET) compared to other climate zones. E_s dominates the central and western TP, with its spatial distribution closely resembling the overall ET. The spatial distributions of E_c and E_w are in line with the distribution of vegetation. High values of E_c (> 200 mm yr⁻¹) and E_w (> 50 mm yr⁻¹) are primarily concentrated in densely vegetated areas such as the Hengduan Mountains in the southeastern TP.

The multiyear average ET for each season on the TP is depicted in Fig. 6, covering spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February). The estimated ET reflects the general seasonal patterns quite accurately. During spring, the average ET is higher than in autumn, ranging from 20 to 250 mm per month in spring and from 20 to 150 mm per month in autumn. This difference can be attributed to the increase in surface water generated by the thawing of permafrost and snow and ice melting as temperatures rise in spring, intensifying surface evaporation processes. Additionally, vegetation transpiration increases during the growing season. In summer, ET exceeds 200 mm per month over most TP, except for large areas in the northwestern TP where ET remains below 100 mm per month. Conversely, in winter, lower ET values are observed primarily in the densely vegetated southeastern region of the TP due to reduced water availability (precipitation) and lower T_a across the entire TP during this season.

Over the TP, the multiyear seasonal ET averages across the entire TP are as follows: 90.79 ± 3.16 mm yr⁻¹ (0.23(± 0.0081) × 10³ Gt yr⁻¹) in spring, 152.05 ± 8.44 mm yr⁻¹ (0.38(± 0.021) × 10³ Gt yr⁻¹) in summer, 71.96 ± 2.86 mm yr⁻¹ (0.18(± 0.0074) × 10³ Gt yr⁻¹) in autumn, and 30.54 ± 1.85 mm yr⁻¹ (0.077(± 0.0047) × 10³ Gt yr⁻¹) in winter. The multiyear average ET is 346.5 ± 13.2 mm yr⁻¹, representing both the mean and standard deviation, which characterizes interannual variability. This corresponds to approximately 0.88(± 0.034) × 10³ Gt yr⁻¹. Among its components, E_s accounted for 292.36 ± 10.39 mm yr⁻¹ (0.74(± 0.027) × 10³ Gt yr⁻¹), E_c amounted to 47.85 ± 3.34 mm yr⁻¹ (0.12(± 0.006) × 10³ Gt yr⁻¹), and E_w is 7.07 ± 2.89 mm yr⁻¹ (0.02(± 0.001) × 10³ Gt yr⁻¹). Notably, E_s constitutes the majority of ET on the TP, exceeding 84 %. Wang et al. (2020) accurately estimated that the water evaporated from all plateau lakes is 0.0517 Gt yr⁻¹. Therefore, utilizing the area-weighted average method, the annual average water evaporation across the entire TP is approximately 0.93(± 0.037) × 10³ Gt yr⁻¹. Furthermore, the TP has an average annual rainfall of about 1.8 × 10³ Gt yr⁻¹, estimated by Jiang et al. (2023). Approximately 53 % of the TP's precipitation returns to the atmosphere through ET (Zhang et al., 2020).

Figure 6 Spatial distributions of the multiyear (1982–2016) mean seasonal ET in (a) spring, (b) summer, (c) autumn, and (d) winter across the TP.

Figure 7 Spatial pattern of the multiyear (1982–2018) mean annual (a) ET (evapotranspiration), (b) E_s (soil evaporation), (c) E_c (canopy transpiration), and (d) E_w (intercepted water evaporation) across the TP.

3.3 Temporal variations in ET across TP

Quantifying variations in ET, both inter- and intra-annual, holds significant importance in understanding monsoon phenomena and studying climate change patterns on the TP. Figure 7 presents the spatial distribution of annual ET and its component trends from 1982 to 2018. These trends exhibit spatial heterogeneity across the TP. The annual ET has seen a significant increase, with rates ranging from 1 to 4 mm yr⁻¹ (p<0.05), primarily in the central and eastern TP, encompassing more than 86 % of the TP. Conversely, there has been a notable decrease in annual ET, with rates ranging from −3 to −1 mm yr⁻¹ in the northwestern TP. Similarly, the trends for E_s mirror those of ET, albeit with slightly lower magnitudes (1–3 mm yr⁻¹, p<0.05). E_c and E_w have shown slightly increasing trends of 0–2 mm yr⁻¹ (p<0.05). When averaged across the entire TP, ET, E_s, and E_c exhibited significant increases during the period from 1982 to 2018, with rates of 0.96 mm yr⁻¹, 0.64 mm yr⁻¹, and 0.44 mm yr⁻¹, respectively (p<0.05; see Fig. 8). Seasonally, positive, and significant trends are observed in all seasons for ET (Fig. 9), with the strongest trends occurring in summer (0.46 mm yr⁻¹). Furthermore, multisource ET products indicate that most regions of the TP have exhibited consistent ET changes over the past 30 years (Yin et al., 2013; Peng et al., 2016; Wang et al., 2018; Ma et al., 2019; G. X. Wang et al., 2020; Li et al., 2021; Ma and Zhang, 2022).

Figure 8 Spatial patterns of the trends (1982–2018) of the annual (a) ET (evapotranspiration), (b) E_s (soil evaporation), (c) E_c (canopy transpiration), and (d) E_w (intercepted water evaporation) across the TP. The stippling on the maps indicates the statistically significant trends (p<0.05).

Figure 9 Time series of the (a) anomalies in the annual ET and its three components and (b) anomalies in seasonal mean ET. The least-squares fitted linear trends are demonstrated by the dashed colored lines.

The rise in ET across the entire TP from 1982 to 2018 can be attributed to the concurrent warming and increased precipitation experienced in this region during the same period. Since the 1980s, the TP has undergone a general trend of greening, warming, and heightened precipitation, as illustrated in Fig. 10. ET has consistently increased over the past 4 decades, but there was a notable shift in climate factors around 2000. From 1982 to 2000, ET showed a continuous increase, accompanied by a rapid decline in wind speed, while the R_n remained relatively stable. However, between 2000 and 2018, there was a sharp decrease in R_n alongside an unchanged wind speed, but ET continued to rise during this period.

Consequently, R_n and wind speed are not the dominant factors driving annual variations in ET. The significant increases in T_a, SM, and precipitation have coincided with the greening of the land surface over the last 2 decades. These factors collectively contribute to the observed increase in ET. In the most recent decade, the substantial growth in SM has emerged as the primary control factor for ET growth.

Figure 10 Time series of the annual anomalies in the (a) NDVI, (b) T_a, (c) R_n, (d) u, (e) SM, and (f) precipitation and their least-squares fitted linear trends during two periods of 1982–2000 and 2000–2018.

In summary, the increase in ET over the TP can be attributed to multiple factors. The rise in available surface water plays a significant role throughout the study. Additionally, there is evidence of a general increase in precipitation across the TP (Fig. 10). The combined impact of warming (shown by T_a in Fig. 10) and vegetation greening (shown by NDVI in Fig. 10) further facilitates the opening of vegetation stomata, promoting increased vegetation transpiration. Warming the land surface and increased wind speeds enhance the efficiency of turbulent water exchange between the land and atmosphere. Furthermore, land surface warming accelerates the melting permafrost and glaciers on the TP. The surface wetting and the thickening of the active soil layer facilitate water transport from the lower soil layers to the upper layers. These environmental changes, such as water availability, precipitation patterns, vegetation dynamics, and temperature trends, all contribute to the increase in ET over the TP.

3.4 Comparison of the MOD16-STM product to other ET products over the TP

We have compared the accuracy of the MOD16-STM product and other available TP region datasets. It is shown in Fig. 11. The MOD16-STM ET model demonstrates high performance on the TP, with an average R² value of 0.87 and an average RMSE of 13.48 mm per month. Wang et al. (2018) evaluated a modified PML model for ET estimation on the TP, reporting R² values exceeding 0.85 and RMSE values lower than 14 mm per month. The spatially averaged ET for 1982–2012 is 378.1 mm yr⁻¹. G. X. Wang et al. (2020) assessed the performance of the generalized nonlinear complementary principle for ET estimation based on flux tower observations from the TP. Their results indicated an R² of 0.93 and a RMSE of 0.40 mm d⁻¹. The spatially averaged ET during 1982–2014 is 398.3 mm yr⁻¹. Han et al. (2021) used a combination of the effective aerodynamic roughness length and the surface energy balance model to estimate ET for the entire TP from 2001 to 2018 (Han-ET). They found good agreement between modeled and in situ-measured values (R²>0.81, RMSE < 14.5 mm per month), and the average annual ET is approximately 496 ± 23 mm, which is higher than the 346.5 ± 13.2 mm obtained in this study (Fig. 12). The discrepancy can be attributed to differences in models and periods used in the two studies. Ma and Zhang (2022) also employed the PML_V2 model to estimate ET on the TP (PML), yielding R² and RMSE values ranging from 0.4 to 0.9 and 0.3 to 0.8 mm d⁻¹, respectively. The 35-year mean annual ET rates from PML-Ma resulted in an average value of 353 ± 24 mm yr⁻¹ for the entire TP. Notably, the proportion of soil evaporation estimated by PML-Ma was approximately 64 % of the total ET, which is lower than the estimated 84 % in this study. The primary reason for this difference may be attributed to variations in land cover classification. The MODIS land cover classification largely categorizes the land surface in the northwestern TP as bare soil, increasing the proportion of soil evaporation.

Figure 11 Taylor diagram of the monthly scale ET dataset validated with flux ET observations.

4 Discussions

4.1 Cross-comparison of the spatial distribution of ET on the TP

A cross-comparison of the multi-year average values of various ET products is conducted to assess the differences and consistency in their spatial patterns. From the spatial distribution of annual average ET (Fig. 12), all the ET products for the TP exhibit a decreasing trend from southeast to northwest, consistent with the transition in surface types from forests to grasslands and bare soil. In the region of the Hengduan Mountains, all products show high values (> 600 mm yr⁻¹) due to the dense vegetation and ample precipitation. However, significant absolute differences are observed among these 15 ET products. There are high differences among the products in the sparsely vegetated western and central TP regions. In the central TP region, the Han-ET product exhibits the highest annual ET (> 600 mm yr⁻¹), while GLDAS-VIC has the lowest (approximately 35 to 50 mm yr⁻¹). In the northwestern TP, EB-ET, GLDAS, MERR-2, and GLEAM-v3.5a products display low values (< 50 mm yr⁻¹), while others range between 100 and 300 mm yr⁻¹. In the extremely arid Qiangtang Plateau, all products show low values due to limited available surface water. ERA-Interim, ERA5-Land, PML-Zhang, CR-Ma, MOD16-STM, and GLEAM-v3.5b have relatively balanced distributions in the central and western TP regions (200–350 mm yr⁻¹). There are high differences in the distribution of ET among the products in the downstream area of the Yarlung Tsangpo River. The spatial resolution of our product is 0.05^∘. This might be the reason that MOD16-STM has low ET for this topographically complex region. It is worth noting that MOD16 ET product has many missing values in the northwestern TP region, making it inadequate for a comprehensive assessment of ET across the entire TP.

Figure 12 Spatial distribution of annual averaged ET on the TP during 2000 to 2014 derived from 15 products.

4.2 ET components' partitioning

It is also necessary to conduct a comparative validation of ET components to enhance the practicality of the data generated in this study. Unfortunately, there are no measured ET component data publicly available at the moment. Comparative validation can be conducted based on existing research findings. Cui et al. (2020) estimated the E_c $/$ ET at the Nagqu station (31.37^∘ N, 91.90^∘ E; 4509 m above sea level) in the central region of the TP using laser spectroscopy and chamber methods. During the observation period, the isotopic-based E_c $/$ ET ranged from 15 % to 73 %, with an average value of 43 %. We calculated E_c $/$ ET from our dataset at the same location and time period (June and July). The values of E_c $/$ ET from MOD16-STM are in the range of 13.1 % to 62.6 %. The average of E_c $/$ ET is 38.4 % ± 4.7 %, which has a difference of 4.6 % relative to isotopic estimation. Our E_c $/$ ET estimation is close to the observation at Nagqu. Guo et al. (2017) also pointed out that E_c constituted less than half of total ET (41 % annually, 29 % during monsoon) in the Magazangbu catchment over the TP.

Moreover, we assess the similarities and differences between MOD16-STM and GLEAM-v3.5a ET components on the TP. Figure 13 shows that GLEAM's E_s values are generally smaller than our estimation throughout the TP region. The most recent results from Zheng et al. (2022) also suggest that the GLEAM product underestimates global E_s outputs. Conversely, GLEAM's E_c values are overestimated in the central and eastern TP. The differences in E_w are minimal because the values in that region are inherently small. Previous research has indicated that in the central TP region, E_s $/$ ET accounts for over 60 % (Cui et al., 2020), and the average E_s $/$ ET ratio across the entire region exceeds 65 % (Wang et al., 2018). The reason for the relatively higher E_s in the central TP is that this region primarily consists of high-altitude grassland as the underlying surface. In the summer, the dominant processes are E_c and E_s, but in the winter E_s becomes the predominant process. Consequently, the proportion of E_s is higher over the entire year. GLEAM's results show that the E_c process predominates in the central TP, which differs somewhat from the findings of this study. Zheng et al. (2022) also indicate that the E_s process predominates in the central TP, exceeding 300 mm yr⁻¹. Therefore, the ET components in this study, when compared with previous research, are more in line with the actual conditions in the TP.

Figure 13 Spatial comparison of ET components and their differences (MOD16-STM minus GLEAM-v3.5a).

4.3 Discrepancy in the estimation of annual ET over the TP

Figure 14 provides a comprehensive overview of the periods covered by various ET datasets and their annual ET estimations for the TP. Yao et al. (2013) estimated TP's ET (PT-Yao) in China using a satellite-driven modified Priestley–Taylor algorithm, constrained by NDVI and apparent thermal inertia derived from temperature changes. Their reported mean annual ET for the TP is approximately 320 mm yr⁻¹. Song et al. (2017) estimated TP's ET (PM-Song) from 2000 to 2010 using the improved Penman–Monteith method and meteorological and satellite remote sensing data at a 1 km spatial resolution. They concluded that the average annual ET on the TP is 350.3 mm yr⁻¹. Additionally, 18 mean annual ET values are estimated using existing ET products (PML-Zhang (Y. Zhang et al., 2019), EB-ET (Chen et al., 2019), CR-Ma (Ma et al., 2019), CMIP6-ssp126 (Eyring et al., 2016), GLDAS-Noah (Rodell et al., 2004), GLASS (Liang et al., 2021), GLEAM-v3.5b (Miralles et al., 2011, 2016), ERAR-Land (Muñoz-Sabater et al., 2021), MTE (Jung et al., 2010), PM-Li (Li et al., 2014a, b), LPJ-Yin (Yin et al., 2013)) and are included for comparison. Han-ET, ERA5-Land, and CMIP6 produce the highest values (> 400 mm yr⁻¹), while LPJ-Yin, GLASS, EB-ET, GLDAS, and GLEAM values are less than 300 mm yr⁻¹. The results demonstrate substantial variability in the TP's estimated mean annual ET values. These differences are influenced by objective factors such as data accuracy, limitations of validation method, and algorithm flaws. The ensemble mean of all datasets yields an annual ET of 348.6 mm yr⁻¹, with the MOD16-STM model's estimation (346.5 mm yr⁻¹) being the closest to this ensemble mean. Overall, the MOD16-STM ET model again demonstrates acceptable TP performance.

Figure 14 (a) The annual mean ET values of 18 datasets. The x axis is the time coverage of the ET datasets, and the y axis is the multiyear mean value. (b) The bars denote the mean values and variations of the annual ET.

4.4 Errors caused by objective factors

The MOD16-STM and other models rely on remote sensing and reanalysis data as primary input sources. However, it is essential to acknowledge the inherent uncertainty in these datasets (Ramoelo et al., 2014). For example, the topsoil water content from satellite data includes some errors (Liu et al., 2021). This indicates that SM from GLEAM may introduce uncertainties into our ET estimation. Some studies have documented the greening of the TP. Figure 10a demonstrated a significant decrease in NDVI after 2000, contrasting with the NDVI changes reported by Wang et al. (2022). This inconsistency highlights the considerable uncertainty in the NDVI data.

Additionally, LST plays a fundamental role in calculating the surface energy balance. Consequently, errors in ERA5 LST can also bring uncertainty to the ET estimation. This study used a threshold value of NDVI (0.25) to categorize land surfaces as bare soil or canopy-covered pixels. This threshold value may misclassify bare soil and grassland on the TP. The land cover mismatches between the reality and the land surface types in the MOD16-STM ET model can also introduce errors into the model simulation.

It is worth noting that flux towers used for validation typically cover areas ranging from a few hundred square meters to several square kilometers. These validation sites' representativeness depends on observation instrument height, turbulence intensity, topography, environment, and vegetation conditions. While site-scale evaluations of the MOD16-STM ET are conducted in this study, it is essential to recognize the uncertainties stemming from the limited number of validation sites. Future research should consider validation across various land cover types, climate zones, elevations, and seasons to provide a more comprehensive assessment of model performance.

While the MOD16 model directly estimates ET, bypassing the need for calculating sensible heat, it still relies on empirical coefficients, particularly those redefined for different soil textures. However, the remaining empirical parameters, such as C_L (the mean potential and stomatal conductance per unit leaf area), can introduce uncertainties into simulation results. Thus, future studies should prioritize parameterizing these empirical factors based on physical processes to reduce simulation uncertainties. It is crucial to consider the influence of physical processes related to deeper soil water and heat transfer on resistance. The MOD16-STM algorithm's accuracy is highly dependent on higher-precision soil moisture products. Since a substantial portion of the TP is covered by permafrost and seasonally frozen soil, assessing soil moisture conditions during freezing and thawing periods becomes challenging. Consequently, it is essential to employ observations during freeze–thaw periods to validate the model's applicability.

In summary, enhancing the model by incorporating physical parameterizations, especially for empirical coefficients, and accounting for the complexities of soil moisture variations in frozen regions will reduce uncertainties in ET simulation results in future research.

5 Data availability

The monthly ET dataset presented and analyzed in this article has been released. It is freely available at the Science Data Bank (https://doi.org/10.11922/sciencedb.00020, Ma et al., 2021) and the National Tibetan Plateau Data Center (TPDC) (https://doi.org/10.11888/Terre.tpdc.271913, Yuan, 2021). The dataset is published under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.

6 Conclusions

In this study, we have developed a 37-year (1982–2018) monthly ET dataset with a high spatial resolution (0.05^∘) for the TP using the newly developed MOD16-STM model. This dataset covers the entire study area with high spatial resolution and a long time series, making it a valuable resource for climate studies. Then, we investigated ET's spatial distribution and temporal trends across the TP. Key findings are summarized in the following.

The ET product generated by the MOD16-STM model exhibits strong performance on the TP. Compared to flux tower observation data, the model achieves high R² and IOA values of 0.83 and 0.93, respectively, with an RMSE of 13.48 mm per month and a modest bias (MB) of 2.58 mm per month. This ET dataset holds potential applications in water resource management, drought monitoring, and ecological studies.

The ET on the TP displays spatial heterogeneity and temporal variations driven by a combination of atmospheric demand and water supply. Generally, annual ET decreases from the southeastern to the northwestern regions of the TP. E_s accounts for over 84 % of the annual ET, and the estimated multiyear mean annual ET on the TP for 1982–2018 is 346.5 ± 13.2 mm. This corresponds to an annual water evaporation of about 0.93 ± 0.037 Gt from the entire TP.

Significant temporal trends are observed in the ET. Most parts of the central and eastern TP exhibit increasing trends of about 1 to 4 mm yr⁻¹ (p<0.05), whereas the northwestern TP shows a decreasing trend of −3 to −1 mm yr⁻¹ (p<0.05). Averaged across the entire TP, the ET increased significantly at a rate of 0.96 mm yr⁻¹ (p<0.05) from 1982 to 2018. This increase in ET over the entire TP from 1982 to 2018 can be attributed to the warming and wetting of the climate during this period.

These findings contribute to a better understanding of the ET dynamics on the Tibetan Plateau and provide a valuable dataset for climate research and related applications.

Appendix A: MOD16-STM model validation at flux site out of the Tibetan Plateau

Table A1 Basic information about the five test sites (which are used to test the relationship between soil surface resistance (r_s) and SM$/{\mathit{\theta }}_{\mathrm{sat}}$ in the MOD16-STM model) and 12 verification sites (used for the MOD16-STM model evaluation at site scale). All the stations are located outside of the Tibetan Plateau.

Table A2 Assessment results of the daily ET (mm d⁻¹) simulated by the MOD16-STM model at the 12 verification sites. The in situ meteorological observation data drive this simulation.

Figure A1 A sensitivity analysis of the impact of the uncertainty in the empirical parameters (ρ_p, ρ_soc, θ_sat,sc, a and b) on the estimation of ET (test in August 2018).

Figure A2 Validation of the consistency between the estimated and the observed values for θ_sat over the TP.

Figure A3 The sensitivity of LE and θ_sat to the changes of m_soc content at different sites.

Figure A4 Time-series comparisons of the daily ET estimated by the MOD16-STM model and observations at the 12 verification sites, which include three grassland sites (US-DK1, US-Fwf, and US-Wkg), three cropland sites (US-IB1, ES-ES2, and IT-Bci), and six evergreen forest sites (CA-Obs, CA-Ojp, CA-Ca2, CA-Ca3, US-DK3, and US-Fuf), respectively.