Spatio-temporal dynamics of evapotranspiration on the Tibetan Plateau from 2000 to 2010

Evapotranspiration (ET) is a key process of the climate system because it links water, energy and carbon cycles. In this study we modified a Penman–Monteith based algorithm to estimate ET on the Tibetan Plateau at a 1 km spatial resolution for the period 2000–2010 using meteorological and satellite remote sensing data. The results showed that the average annual ET on the Tibetan Plateau was 350.3 mm year−1 and decreased from the southeast toward the northwest. The highest ET value was found in open water bodies (680.9 mm year−1) and the lowest ET value was found in open shrubland (254.0 mm year−1). Overall, the inter-annual ET decreased from 2000 to 2010 and there was significant negative ET trend over 42% of the region, primarily in the northwest of the Tibetan Plateau. Relative humidity was the dominant factor in controlling long-term variations of ET in the arid northwest plateau. But under moist conditions, leaf area index or temperature drove ET. In addition, P-ET on the Tibetan Plateau significantly increased and about 37% of the region showed strong positive P-ET trend primarily in the central of the Tibetan Plateau. The positive P-ET trend in four seasons suggested that the Tibetan Plateau might have become wetter during the past decade.


Introduction
Evapotranspiration (ET), which is the sum of water lost from various land surfaces through evaporation, from the stomata of plants through transpiration, and from snow cover through sublimation, plays a key role in land surface energy balance and hydrological cycle (Oki and Kanae 2006, Trenberth et al 2009, Vinukollu et al 2011. The hydrological cycle is expected to accelerate with global warming and there has been an increasing trend in ET as revealed by a review of historical trends in hydrologic variables (Huntington 2006). However, our understanding of the relative magnitudes and spatial patterns of ET changing trends is limited due to sparse ET observation networks. An accurate assessment of temporal and spatial dynamics of ET in response to climate change is essential for projecting potential changes in the global hydrological cycle.
The water transfer processes involved in ET depend on atmospheric demand and moisture supply (Jung et al 2010, Seneviratne et al 2010. And there are numerous factors driving ET, including radiation, vapour-pressure deficit, temperature, soil moisture and stomatal conductance (influenced by leaf area and vegetation coverage) (Hu et al 2013, Wang et al 2010, Jung et al 2010, Chen et al 2014, Seneviratne et al 2010, making accurate assessment of ET a challenge. ET is routinely measured at flux towers, but the high spatial heterogeneity of ET limits the application of flux tower measurements to regional or global assessments (Anderson et al 2003, Wang et al 2010.
Several satellite-based approaches are available to quantify ET at large scales, including (1) empirical methods (Jung et al 2010, Jin et al 2011, Zhang et al 2007, Yin et al 2013; (2) surface energy balance methods, such as Surface Energy Balance System (SEBS; Su 2002); (3) physical models, such as Penman-Monteith algorithm (Monteith 1965, Penman 1948. The empirical methods, based on statistical relationships between measured ET and eco-environmental parameters, are convenient for ETestimation on global or regional scales. However, uncertainties in the empirical coefficients for different ecosystems or regions pose a challenge for the accurate assessment of ET (Jin et al 2011). The surface energy balance models are theoretically sound, but their performances are constrained by the estimation of sensible heat flux and soil heat flux (Cleugh et al 2007). Among these modeling methods, the Penman-Monteith algorithm and its varieties have been widely used to estimate ETat regional or global scales (Cleugh et al 2007, 2011, Zhang et al 2009, Yuan et al 2010, due to its more mechanistic representation of ET processes by introducing aerodynamic and surface resistances. Mu et al (2011) took a further step to improve the Penman-Monteith algorithm by developing more realistic formulations of stomatal conductance, aerodynamic resistance and boundary layer resistance. This modification of the Penman-Monteith algorithm had been proved to be biophysically sound and rigorous for estimating ET at regional scales (Jung et al 2010, Vinukollu et al 2011, Li et al 2014, Chen et al 2014. The Tibetan Plateau (TP) is very sensitive to climate change (Zheng 1996, Immerzeel et al 2010. As indicated by previous studies (Liu and Chen 2000, Wu et al 2007, the warming trend over the plateau might have led to glacier retreat (Yao et al 2004) and permafrost degradation Zhang 2010, Yang et al 2010). These warming-induced changes will greatly affect the ET processes over the plateau. Until now, most of the ET studies on the TP focused on ET at river basin (Zhang et al 2007, Xue et al 2013 or site scales (Hu et al 2009, 2013, Zhu et al 2013, Yao et al, 2008, Gu et al 2005. Only a few studies conducted large-scale ET assessments on the TP. For example, Yin et al (2013) investigated ET dynamics with Lund-Potsdam-Jena dynamic vegetation model and concluded that the increasing ET trend over the last 30 years was linked with increased precipitation. To get a more accurate ET assessment on the TP, a continuous ET record with a high spatio-temporal resolution is needed. Furthermore, the assessment of ET dynamics on the TP is also limited by the lack of our knowledge on the ET controls and net precipitation (precipitation minus ET) at regional scales.
The objectives of this study are: (1) to evaluate the performance of the modified Penman-Monteith based algorithm; (2) to quantify the spatial and temporal patterns of ETamong different land-cover types on the TP from 2000 to 2010; (3) to investigate the dominant environment variables; and (4) to explore the response patterns of both ET and net precipitation to climate change during the period.
2. Methods and data sources 2.1. Study site The Tibetan Plateau (TP) (75°E-105°E, 25°N-40°N), the highest (average elevation >4000 m) and largest (>2.5 million km 2 ) highland in the world, is located in western China and covers most of the Tibet Autonomous Region and Qinghai province, as well as small portions of western Sichuan province, southwestern Gansu province, northern Yunnan province and Xinjiang Autonomous Region (figure 1). Most areas of the TP are characterized with an arid/ semiarid climate. The mean annual precipitation is 473 mm, ranging from 1000 to 50 mm. The mean annual temperature is 3.8°C, ranging from À15°C to 10°C. The patterns of precipitation and temperature result in an increase in aridity from southeast to northwest of the TP. More than 60% of the plateau is alpine grasslands, which includes alpine meadows, alpine steppes and alpine desert. The TP is rich in lakes, glaciers, and wetlands, and is the main source of several major rivers in Asia.

The modified Penman-Monteith algorithm.
Our ET algorithm is a revision of the algorithm proposed by Mu et al (2011) which is based on the well-known Penman-Monteith equation (Penman 1948, Monteith 1965) (hereafter called PM-Mu (2011)) where lET (W m À2 ) is the latent heat flux (LE), l (J kg À1 ) is the latent heat of vaporization, R net (W m À2 ) is the net incoming radiation to the land surface, G flux (W m À2 ) is the soil heat flux, r (kg m À3 ) is the air density, C p (J kg À1 K À1 ) is the specific heat capacity of air, VPD (Pa) is the vapour-pressure deficit, D (Pa K À1 ) is the slope of the saturate vapor pressure curve, g is the psychrometric constant, r a (s m À1 ) is the aerodynamic resistance and r s (s m À1 ) is the surface (or canopy) resistance. Similar with PM-Mu (2011), we estimated evaporation from soil (lE SOIL ), wet plant (lE wet À c ) and transpiration from plant (lE trans ) separately. And details of estimation are as follows.
The R net was calculated as equation (2). where a is MODIS albedo, R sw (W m À2 ) is the incoming shortwave radiation and s (5.67 Â 10 À8 W m À2 k À4 ) is Stefan-Boltzmann constant. Different from PM-Mu (2011), we modified the G flux with following equations recommended by Frempong (1983): For water surfaces, For other land-cover types, where G is the ratio of soil heat flux to the net radiation. The values of G for soil and canopy are 0.315 (Kustas and Daughtry 1990) and 0.05 (Monteith 1973), respectively. The fractional vegetation cover information (f c ), estimated from MODIS FPAR (the Fraction of Absorbed Photosynthetically Active Radiation), was used to partition the available energy (R net À G flux ) between soil and vegetation. As As in PM-Mu (2011), we separated dry canopy surface from the wet by F wet to calculate evaporation from the wet canopy surface and moist bare soil surface.
The evaporation from wet canopy surface was calculated as equation (7).
where rhc (s m À1 ) is the wet canopy resistance to sensible heat, rrc (s m À1 ) is the resistance to radiative heat transfer through air, gl_sh (s m À1 ) is leaf conductance to sensible hear per unit leaf area index (LAI), gl_e_wv (m s À1 ) is leaf conductance to evaporated water vapor per unit LAI, and T air is the air temperature. The plant transpiration was calculated as equation (9).  Environ. Res. Lett. 12 (2017) 014011 The value of rr was equal to the rrc in equation (8).
where C c is the canopy conductance, which is derived from stomatal conductance (Gs1), leaf cuticular conductance (Gcu) and leaf boundary-layer conductance (Gs2), P a is the atmospheric pressure, m(T min ) and m(VPD) are multipliers that limits potential stomatal conductance, which were calculated using the same equations in Mu et al (2011). Based on the characteristics of TP, we calibrated VPD close , the biome-specific critical value of VPD at which canopy stomata are completely closed, and C L , the mean potential stomatal conductance per unit leaf area, to be 3600 Pa and 0.0065 m s À1 for grasslands based on the in-situ measurements of ET. The other biome property parameters were taken from Mu et al (2011). Soil evaporation (lE SOIL ) was calculated in equation (12).
where b was set as 200 in PM- Mu (2011), and was revised as 50 in the improved algorithm. Details in calculating r as , r tot and other parameters are the same with Mu et al (2011).

Evaporation from snow.
Since snow pack is important for water resources and hydrology, we considered the evaporation over snow based on the characteristics of TP. And Penman equation was used: where r a is the aerodynamic resistance and expressed as a power function of wind speed at height 2 m (U 2 ): The transpiration over snow covered plants was assumed to be negligible since stomata closes at freezing temperature (Vinukollu et al 2011).

Evaporation from water bodies.
Besides the snow evaporation, we also considered the evaporation over water bodies. And the Priestley-Taylor equation (Priestley and Taylor 1972) was used: where a is calculated by the equation recommended by Er-Raki et al (2010): Finally, the total daily ET was the sum of different parts: For water surfaces, For other land-cover types, The model time step was set to 8 days because the satellite products were calculated as 8-day composites.

Eddy covariance data.
To validate the revised model, we used in-situ measurements of LE at two eddy covariance sites. One is located at Haibei alpine grassland station in Qinghai (37°36 0 N101°20 0 E; 3160 m a.s.l.). The other is located at Damxung grassland station, north of Lhasa city, Tibet (30°51 0 N, 91°05 0 E; 4333 m a.s.l.). The study periods extended from January 2003 to December 2004 for Haibei station and from July 2003 to December 2005 for Damxung station. We aggregated the half-hourly flux observations into daily values, and those days with more than 20% missing half-hourly data were excluded.
2.3.4. Pre-processing input data. The satellite products contain cloud-contaminated or missing data. All missing or unreliable data were processed using method proposed by Zhao et al (2005) with a software package MATLAB10.0. The daily values of meteorological measurements and R sw were averaged to 8-day interval and spatially interpolated to 1-km with inverse distance weighting method (Lu and Wong 2008).
The tower measured data are given every 30 min, and for each 30-min time period, ET (mm 30 min À1 ) was calculated as: l ¼ ð2:501 À 0:002361 Â T air Þ Â 10 6 ð20Þ 2.3.5. Model performance and trend analysis. Root mean square error (RMSE), mean bias error (MBE) as well as observation (R 2 ) were calculated to analyze the difference between the modeled simulations and observations.
where ET obs and ET mod is observed and modeled value, respectively; and n is the total number of observations. Linear trend analysis was used to quantify the inter-annual, seasonal and each pixel trend of the ET (y t ).
Statistically significant differences were set as p < 0.1 unless otherwise stated. Correlation (r) was calculated to determine which environmental factor explained most of the variation of ET.

Model performance
The ET averaged over the 5 Â 5 1-km 2 pixels surrounding each site is compared with the tower ET observations. Two tower sites with available LE are used to compare ET seasonal variability between our modeled LE flux and MOD16A2. Compared with MOD16A2, the ET estimates from our revised model are much lower in winter, but much closer to the observations in summer. For these two tower sites, the revised model shows generally favorable agreement with the tower observations and captures the seasonality of observed LE (figure 2). Overall, the modified model reduces MBE from 0.03 W m À2 using the Mu et al (2011) algorithm to 0.007 W m À2 , improves LE accuracy (RMSE) from 25.7 W m À2 to 15.8 W m À2 , and enhances LE correlation (R 2 ) from 0.36 to 0.76 (figure 3).

Spatial patterns of ET
The mean annual ET on the TP estimated by the modified PM-Mu (2011) is 350.3 ± 139.9 mm year À1 and the spatial ET distribution corresponds to variations in precipitation, available energy and temperature on the plateau. The largest ET exhibits in the southeast whereas small ET exhibits in the dry and lower temperature regions especially in the northwest. Spatially, the distribution ET estimated by modified model matches well with MOD16A2 (figure 4), although spatially-averaged ET of the MOD16A2 is about 111.1 mm higher.
Furthermore, the spatial ET distribution is also associated with distributions of the major land-cover types (figure 5). As expected, the highest average ET value is found in open water bodies which is 680.9 mm year À1 ; and open shrubland have the lowest ET, which is 254.0 mm year À1 . The ET values for other landcover types lie between the two extremes. ET of croplands is generally higher than that of grasslands, which due to irrigation in croplands. Temporally, in spring (March-May), most parts of southeast TP show significantly high ET values because the forest ecosystems, distributing in these regions with mild temperatures and moist conditions during the spring, start to grow. ET goes to the highest level in summer (June-August) because of high precipitation, dense vegetation and intensive radiation. In autumn (September-November), the ET drops as vegetation senesced and radiation declined. The majority of the TP has very low ET in winter (December-February) because of low available energy, temperature, and stomatal conductance.      7(a)). As a whole, ET significantly decreases on TP for the 2000-2010 period (4.7 mm year À1 , p < 0.001) ( figure 7(b)). We calculated the correlation between environmental factors (R sw , LAI, RH, T air , FPAR, and Albedo) and ET in each pixel. And the results show that mean annual RH (r = 0.72), LAI (r = 0.62) and T air (r = 0.58) are important factors that affect ET variations overall the TP from 2000 to 2010. Generally, for dry conditions (northwestern Tibetan Plateau), ET contributes to RH, but under moist conditions (southeastern Tibetan Plateau), LAI or T air drives ET (figure 8).

Changes in net precipitation
Theoretically, ETshould be less than precipitation over a relative long time period. To validate the ET results and assess the water balance changes, we calculated the net precipitation (P-ET). The results show that ET is less than P spatially and a mean annual P-ET value of 112.2 mm year À1 from 2000 to 2010 is estimated on the TP. Spatially, P-ET presents a similar distribution pattern with respect to ET. The P-ET is generally low in arid regions of the northwest TP, whereas relatively high P-ET is located downstream of Brahmaputra (figure 9).
Contrary to ET, the inter-annual P-ET increases significantly (5.0 mm year À1 , p < 0.1) and varies between 68.5 and 150.1 mm year À1 from 2000 to 2010 ( figure 10(b)). Regionally, about 37% areas showing strong positive P-ET trends primarily occur in the central of the TP, and only 2% areas present strong negative P-ET trends ( figure 10(a)). Positive P-ET trends generally occur over grasslands and open shrublands, while negative trends occur within forests regions of the domain (table 1).
In addition, both regional average ET and P-ET show seasonal changes during the 11-year period (figure 11). The ET shows significant (p < 0.1) negative variability for all four seasons, with the largest ET decreases in summer followed by spring and autumn. Opposite to ET seasonal pattern, P-ET presents significant positive variability (p < 0.01) in winter, with slightly increasing P-ET in spring, summer and autumn. The positive P-ET trends in four seasons indicate that the TP became wetter during the past decade.

Model performance and uncertainties
The accurate ET estimation at the diurnal time scale in this study demonstrated satisfactory predictive capabilities of the modified PM-Mu (2011) model. Compared with the MODIS ET product, our modified model illustrated better performance at the tower sites (figures 2 and figure 3). There were two reasons that resulted in the better performance of the modified model. On the one hand, the parameters in this study were locally validated while those of the MODIS ET were only tested by AmeriFlux eddy covariance tower sites. On the other hand, the different input driving variables were used, i.e. meteorological data from CMA was used in this study while re-analyses meteorological data from GMAO was used in MODIS ET. The GMAO climate data tended to have high uncertainties in winter (Zhao et al 2005).
Potential sources of uncertainties in ET estimates could come from the algorithm itself. Although many studies had confirmed the good performance of the Penman-Monteith model and contributed to its parameterization (Cleugh et al 2007, 2011, Zhang et al 2009, Yuan et al 2010, Li et al 2014, Chen et al 2014, not claiming completeness, some improvements can be made to the modified PM- Mu (2011). For example, as the source of water cycle, P was not an input parameter to constrain ET. In addition, since soil moisture was the main limiting factor for the plant physiological processes, the effect of soil water condition on canopy stomatal conductance should add to the model (Hu et al 2013). Another limitation of the ET algorithm was that the effect of high CO 2 and nitrogen induced stomatal closure on stomatal resistance had not been considered in estimating ET. Furthermore, some key parameters were given with empirical values, which may cause large deviation.
In addition, uncertainties in ET estimates could also come from the input driving data. Due to inaccessibility, the complex terrain, sparse meteorological stations and higher elevation (>4500 m) in the western TP, the ET estimates may be negatively impacted by the inherited interpolation meteorological inputs. In addition, uncertainties from MODIS LAI/FPAR, albedo and land-cover data can introduce biases to ET estimates. For example, MODIS LAI tended to be higher than measurements  and the bias of MODIS land-cover data were in range of 20%-30% (Strahler et al 2002). Moreover, evaporation from open water bodies were likely underestimated due to the inability of 1-km MODIS land-cover information to resolve rivers on the TP.
Generally, validation of the estimations with eddy flux LE measurements posed challenges. Firstly, uncertainties may be introduced by comparing measured ET with the estimated from the 5 Â 5 1-km 2 pixels due to differences in tower footprints (Anderson et al 2003). Second, the eddy covariance flux towers had a problem with energy balance closure because of complexity in wind variation, footprint representation, and sampling variability (Twine et al 2000). The energy balance was rarely closed resulting in non-closure on the order of 20%-30%. Therefore, the best model would subsequently explain only an R 2 of around 0.70-0.80, which was generally what our model produces.

Variations of ET on the Tibetan Plateau
The distribution pattern and ranges of variations in ET in this study were consistent with previous studies (Chen et  Although the assessments of ET had been made by numerous models, ET inter-comparison studies showed large uncertainties in global or regional ET estimates (Vinukollu et al 2011, Chen et al 2014. For example, Chen et al (2014) compared eight ETmodels, including five empirical and three process-based models, and found that the process-based models performed better than the empirical models. Although eight models indicated similar spatial patterns, there were substantial differences in the magnitude and inter-annual variability of ET, which may be because of differences in model-dominated variables. The modified PM- Mu (2011)     Generally, the ET variability is either induced by atmospheric energy demand or by terrestrial moisture-supply (Jung et al 2010). ET will respond to atmospheric demand (e.g. temperature) in wellwatered regions, where warming is expected to promote increase in ET and precipitation leading to a general acceleration of the water cycle (Huntington 2006). However, moisture-supply will turn to be the dominant driving factor if the soils are too dry. Study of Jung et al (2010) demonstrated that the warminginduced increasing ET disappeared after the last big El Niño event in 1998. Similar with previous studies (Wang et al 2010, Yao et al 2013, Li et al 2014, Chen et al 2014, the declining ET in our study was not consistent with the expected acceleration of the hydrological cycle. Our TP-ET estimates suggested that the inter-annual ET decreased from 2000 to 2010 with a linear trend of 4.7 mm year À1 , coinciding with regional permafrost thawing Zhang 2010, Yang et al 2010) and vegetation degradation (Harris 2010). It was found that the largest regional contributions to the declining trend since 2000 originated from the northwest of the TP in which ET was limited by moisture-supply, indicated here by RH ( figure 8). The decreasing RH also contributed to inter-annual and seasonal decreased variations in ET. In these dry regions, lower ET was expected to feed back to the atmosphere and increased atmospheric dryness. However, under relatively moist conditions (southeast of TP), dependencies of ET on LAI and T air appeared to be largely independent of RH (figure 8). Lower LAI and T air indicated less leaves and lower canopy conductance, thus decreased in ET (Wang et al 2010, Chen et al 2014. The net precipitation (P-ET) describes moisture flux between the land surface and the atmosphere. Changes in P-ET could affect local atmospheric moisture content and regional climate. Our results showed that P-ET close to zero in northwest, which confirmed the droughts in this region (figure 9). And relatively small difference between P and ET implies a strong sensitivity of water supplies to relatively small changes in Por ET. In addition, our findings of increasing water supply since 2000 coincide with increase in runoff and lake level on the TP (Zhang et al 2011).

Conclusions
ET remains one of the greatest unknowns within the global energy, water and carbon cycles. To better understand the spatio-temporal dynamics of ET on the Tibetan Plateau, we have modified the PM-Mu (2011) model to estimate the Tibetan Plateau's ETon an 8-day time scale. Overall, the mean annual ET was 350.3 mm year À1 and decreased from southeast to northwest of Tibetan Plateau. The inter-annual ET showed a decreasing variability from 2000 to 2010. High negative ET rates were found in the northwest of the Tibetan Plateau. RH was the dominant factor in controlling long-term variations of ET in arid northwest plateau. However, in the moist southeast plateau, variations of ET were primarily controlled by T air and LAI. Contrary to the ET variability, P-ET presented significant positive trends all season and an increasing pattern from 2000 to 2010. This indicates that the Tibetan Plateau became wetter during the past decade.