Attribution of Long-Term Evapotranspiration Trends in the Mekong River Basin with a Remote Sensing-Based Process Model

: Using the Global Land Surface Satellite (GLASS) leaf area index (LAI), the actual evapotranspiration (ET a ) and available water resources in the Mekong River Basin were estimated with the Remote Sensing-Based Vegetation Interface Processes Model (VIP-RS). The relative contributions of climate variables and vegetation greening to ET a were estimated with numerical experiments. The results show that the average ET a in the entire basin increased at a rate of 1.16 mm year − 2 from 1980 to 2012 (36.7% of the area met the 95% signiﬁcance level). Vegetation greening contributed 54.1% of the annual ET a trend, slightly higher than that of climate change. The contributions of air temperature, precipitation and the LAI were positive, whereas contributions of solar radiation and vapor pressure were negative. The effects of water supply and energy availability were equivalent on the variation of ET a throughout most of the basin, except the upper reach and downstream Mekong Delta. In the upper reach, climate warming played a critical role in the ET a variability, while the warming effect was offset by reduced solar radiation in the Mekong Delta (an energy-limited region). For the entire basin, the available water resources showed an increasing trend due to intensiﬁed precipitation; however, in downstream areas, additional pressure on available water resources is ex-erted due to cropland expansion with enhanced agricultural water consumption. The results provide scientiﬁc basis for practices of integrated catchment management and water resources allocation.


Introduction
As a result of the impacts of global warming and human activities, the global water cycle has been intensified, leading to changes in global precipitation [1], evapotranspiration (ET) [2,3] and runoff [4]. More than half of the global absorbed solar radiation is used for ET, which reintroduces approximately 60% of land surface precipitation back into the atmosphere [5]. However, the mechanisms responsible for the variations in ET are spatially heterogeneous; thus, investigating ET variations and their underlying mechanisms can improve our ability to simulate land surface processes and predict future water cycle alterations.
Terrestrial ET variations are impacted by changes in the atmospheric evaporative demand (AED), water supply and physiological properties of vegetation. Observational and experimental evidence has shown that climate factors, including precipitation, solar radiation, temperature and wind speed, control both the AED and the water supply for terrestrial ET. Research has found that in humid areas, long-term ET variations are driven primarily by changes in incident solar radiation, whereas ET variations in arid areas are controlled predominantly by the soil water supply, which is linked to precipitation [6]. In addition to climate factors, elevated atmospheric CO 2 concentrations can also induce the mean annual precipitation in the MRB shows significant spatiotemporal heterogeneity. More than 90% of annual precipitation falls during the rainy season, which spans from mid-May to early October. The highest annual rainfall occurs in the mountainous regions of Laos, which receives mean annual precipitation of 3200 mm, whereas the lowest annual rainfall reaches approximately 1000 mm in northeastern Thailand [27]. The land cover types in the upper catchment are primarily tundra and grassland, whereas farther downstream, the natural vegetation is dominated by evergreen broadleaved forest; in the lower catchment and Mekong Delta, cropland is the major land-use type (Figure 1).

Meteorological Data
Meteorological data, including temperature, surface air pressure, longwave radiation, shortwave radiation, precipitation, relative humidity and wind speed, were collected on a daily scale from 1980 to 2012. The meteorological data were provided by the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) [28,29], with a spatial resolution of 0.5 × 0.5 • . The meteorological data used in this study were resampled to 5 × 5 km with the bilinear resampling method. The air temperature, precipitation and wind speed data were resampled at sea level, and then were corrected to the real elevation based on a digital elevation model (DEM) with a 5 km resolution after bilinear resampling (Appendix A).

Land Surface Characterization Data
The land surface characterization data included topographic maps, land-use maps and soil texture maps. DEM data were obtained from the GTOPO30 global DEM (https://lta.cr. usgs.gov/GTOPO30) and resampled to a 5 km resolution with cubic Lagrange interpolation. Land cover data (MOD12Q1) were downloaded from the Moderate Resolution Imaging Spectroradiometer (MODIS) website (http://modis.gsfc.nasa.gov/data/) with a spatial resolution of 500 m and a temporal resolution of one year. Land cover data from 2000 to 2012 were used in this study and were resampled to 5 × 5 km with the majority resampling method [30]. The soil texture map was obtained from the Harmonized World Soil Database v 1.2 provided by the Food and Agriculture Organization of the United Nations.

Remote Sensing Data
The remote sensing data included the Global Land Surface Satellite (GLASS) LAI product, Global Land Evaporation Amsterdam Model (GLEAM) ET data and Gravity Recovery and Climate Experiment (GRACE) data. Based on general regression neural networks, the GLASS LAI product (http://glass-product.bnu.edu.cn/), which has a spatial resolution of 0.05 • and a temporal resolution of 8 days, was retrieved from MODIS and Advanced Very-High-Resolution Radiometer (AVHRR) land surface data. Compared with MODIS LAI and CYCLOPES LAI data, the GLASS LAI product is spatially more complete and temporally more continuous [31]. Monthly actual ET data from GLEAM v3.2 (https://www.gleam.eu/), which have a spatial resolution of 0.25 × 0.25 • , were used to validate the ET a results. Based on the function between the soil moisture stress and potential evaporation calculated with the Priestley-Taylor equation, four evaporation components can be obtained from GLEAM, namely, interception loss, bare soil evaporation, plant transpiration and evaporation from water bodies and regions covered by ice and/or snow. With the bilinear resampling method, the GLASS LAI data and actual ET data from GLEAM from 2000 to 2012 were resampled to 5 × 5 km. In addition to the GLEAM ET a , monthly ET a data calculated from the water balance equation were also used to validate the ET a results at the basin scale: where P and R are the precipitation and the runoff respectively, and ds/dt is the derivative of the terrestrial water storage anomaly (TWSA) with respect to time. The TWSA can be reflected by GRACE RL05 data, which can be retrieved from the website of the Jet Propulsion Laboratory, California Institute of Technology (https:// grace.jpl.nasa.gov/data/get-data/monthly-mass-grids-land/). After removing the effects of the atmosphere, oceans and tides, monthly level-2 GRACE RL05 data from 2003 to 2012 were used in this study. To restore signal losses arising from the sampling and post processing of GRACE data, the data were recorrected with a scaling factor approach [32], and the data gaps were filled using cubic Lagrange interpolation.

Model Introduction
ET a was simulated by the Remote Sensing-Based Vegetation Interface Processes (VIP-RS) model [33,34]. The total ET a is composed of evaporation due to canopy interception (E i ), transpiration from vegetation (E c ) and soil evaporation (E s ). Transpiration (E c ) is estimated based on potential transpiration (E cp ) and is limited by water conditions (f w ) and temperature (f t ), as follows: where f w is the stress function of the atmospheric water vapor pressure deficit calculated using the algorithms proposed by Mu et al. [35] and f t is the stress function of the air temperature calculated using the algorithms proposed by Zhang et al. [36]. E cp is calculated using the Penman-Monteith (PM) equation as follows: where R nc is the net radiation absorbed by the canopy, ∆ represents the slope of the saturated vapor pressure curve versus air temperature,γ, ρ, c p and D represent the psychrometric constant (hPa • C −1 ), air density (kg m −3 ), specific heat capacity of air (J kg −1 • C −1 ) and vapor pressure deficit (hPa), respectively. Additionally, η is the ratio of the minimum stomatal resistance of a natural functional plant type to that of a reference crop, λ is the latent heat of the vaporization of water (J kg −1 ) and r a is the aerodynamic resistance between the canopy and the reference height (s m −1 ). f c is the fractional vegetation cover, and f CO 2 is the stress function of atmospheric CO 2 [37]. These two factors are calculated as follows: where CO 2 is the atmospheric CO 2 concentration (ppm), β is an empirical constant, being 0.6 to 1.2, taken here as 0.7 [38], and VI max and VI min represent the vegetation index values under full vegetation cover and bare soil conditions, respectively (set to 0.95 and 0.01). Additionally, f PAR max and f PAR max are the maximum and minimum f PAR corresponding to VI max and VI min respectively, and are set to 0.95 and 0.001. LAI max is the maximum LAI during the growing period (assumed to be 6.0) [39]; the LAI, which influences the partitioning of energy between the soil and canopy, is retrieved from remote sensing data. The evaporation due to canopy interception (E i ) equals the potential evaporation from the wetted surface. Soil evaporation (E s ) is the minimum value of surface potential evaporation (E sp ) and soil moisture exfiltration (E ex ) [40]: where R ns is the net radiation absorbed by the soil surface, G is the soil heat flux (MJ d −1 ) [41], r as is the aerodynamic resistance between the reference height and the soil surface (s m −1 ), Γ c (set as 0.05) [42] and Γ s (set as 0.315) [43] are the ratios of G to R nc for full vegetation canopy and bare soil, respectively.

Model Implementation and Verification
Driven by the daily meteorological data (average, maximum and minimum air temperatures, atmospheric pressure, humidity, wind speed, precipitation, longwave downwelling radiation and shortwave downwelling radiation), the daily ET a from 1980 to 2012 was simulated by the VIP-RS model at a spatial resolution of 5 km. As the initial condition of soil moisture was not available, the spatial pattern of soil moisture after one-year run was used as the initial state [44], and daily ET a from 1981 to 2012 were used for the subsequent analysis. The tendency of annual ET a and its significances were explored based on the Theil-Sen estimator and Mann-Kendall (M-K) test.
The simulated monthly ET a was consistent with the actual ET from GLEAM ( Figure 2a) and yielded an R 2 value of 0.84. The mean bias and root mean square error (RMSE) between the simulated and GLEAM ET a were 4.5 mm month −1 and 8.17 mm month −1 , respectively. At annual scale, the GRACE-derived ET a and simulated ET a showed acceptable consistency with an R 2 of 0.55. However, the GRACE-derived ET a was 25.72 mm year −1 lower than the simulated ET a , indicating that the VIP-RS model slightly overestimated the yearly ET a in the MRB.

Attribution of ET a to Climate Change
The factor separation methodology [45] is a popular technique for evaluating the contribution of each input factor to a chosen response variable. To quantitatively distinguish the effects of climate change, carbon dioxide enhancement, vegetation greening and their interactions on the ET a trend, thirty simulation experiments were designed based on the factor separation method. In each simulation experiment, one or more variables were varied according to the observation records, while the other variables were varied according to the control conditions [46]. The control simulation (f(con)) is driven by the mean climatology, average atmospheric CO 2 concentration and vegetation dynamics from 1982 to 1990. The ET a result from the control simulation is the expected value under the average climate and vegetation conditions of the 1980s. In the control simulation, the CO 2 concentration is the average CO 2 value from 1982 to 1990. Except for precipitation, the climate drivers and LAI in the control simulation are the multiyear means for individual days (expressed as Julian days) from 1982 to 1990. To eliminate the impact of variation in rainfall patterns on the resulting ET a , the precipitation in each year (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)) is used in the control simulation, and the final ET a result from the control simulation is the average value of each year. The thirty performed simulation experiments included one control simulation, seven simulations (f(var1)), in which only one of the seven factors was varied each time, and twenty-one simulations (f(var1_var2)), in which two of the seven factors were varied for each model run (Appendix B, Table A1). Simulation experiment f(var1) denotes the simulation driven by the variable var1, which was varied according to the ISIMIP records, and the other variables were set to the control conditions. Similarly, simulation experiment f(var1_var2) is the simulation driven by the variables var1 and var2, which were varied according to the ISIMIP records, while the other variables were set to the control conditions. The seven factors used in the simulation experiments are temperature, precipitation, net radiation, vapor pressure, wind speed, LAI and atmospheric CO 2 concentration.
The main effect of each factor on the ET a trend is the difference between the ET a result from the simulation with only one variable factor and the result from the control simulation. For example, the main effects of variable one (var1) and variable two (var2) on the ET a trend are as follows: The two-factor interactive effect of variables one (var1) and two (var2) is calculated by subtracting the main effects of var1 and var2 from the combined effect of var1 plus var2. For example, the two-factor interactive effect of var1 and var2 (E(var1_var2)) is expressed as: where E(var1) and E(var2) are the effects of var1 and var2 on the ET a trend, respectively. The land-use transfer matrix, which is used to analyze spatial-temporal changes in land use, is built with land-use data (with a spatial resolution of 500 m) from 2000 and 2012. To assess the contributions of land-use changes to ET a , the control simulation (f(land_con)) was driven by the land-use data from 2000. For certain land-use types, simulation f(land_type1) was driven by one land-use type that was varied according to the records from 2010, while the other land-use types remained unchanged according to the year 2000. The effect of a change in a particular land-use type (land_type1) on the ET a trend is the difference between f (land_type1) and f (land_con). Considering that ET a was simulated at a 5 km resolution and that the land-use data have a spatial resolution of 500 m, to assess the contributions of land-use changes to ET a , the ET a value in each sub-grid and the contributions of different sub-grid-scale land-use types to the grid box average were calculated (Equation (14)) to reduce the error in the ET a simulation: where ET ai is the ET a value in each sub-grid (5 km resolution) with land-use type i ( Figure 1) and V i is the proportion of land-use type i in each sub-grid with a 5 km resolution. V i is obtained from land-use data with a spatial resolution of 500 m.

Spatial-Temporal Variation of Actual ET
ET a is principally controlled by the AED (represented by ET p ) and precipitation. The annual ET a in the MRB showed remarkable spatial heterogeneity with an average value of 896 ± 34 mm year −1 and increased gradually from north to south with a range of 300 to 2000 mm year −1 . From 1982 to 2012, 70% of the study area showed an increasing ET a trend with an average increase rate of 1.16 mm year −2 .
In the upper basin, where the precipitation is less than 800 mm year −1 and the AED is weak, the annual ET a is less than 400 mm year −1 (Figure 3). Due to rising temperatures, glacial melting in the upper basin has accelerated, which has been accompanied by increased precipitation, and the annual ET a has shown a significant increasing trend over more than 13.7% of the study area (p < 0.01, Figure 4c).
In the central and lower reaches of the basin, the situation is relatively complex. In evergreen broadleaf forest, which have abundant rainfall (annual precipitation > 1700 mm year −1 ), the spatial pattern of the average annual ET a is highly consistent with the average annual ET p (Figure 3), whereas in cropland/natural vegetation mosaic areas, which have an average annual precipitation of less than 1200 mm year −1 , approximately 75% of the rainfall returns to the atmosphere through evaporation; as a result, the spatial pattern of ET a is similar to that of precipitation. A significant positive temporal trend of ET a was found in the Mekong Delta and the southeastern basin (Figure 4a), where precipitation and LAI display increasing trends. Significant positive trends of ET a were found in 36.7% of the MRB (p < 0.05), most of which corresponded to savanna and cropland with high increasing rate of annual ET a (Figure 4c). The highest increasing rate was found in cropland (2.53 mm year −2 ), followed by grassland (1.8 mm year −2 ), savanna (1.68 mm year −2 ) and evergreen broadleaf forest (1.25 mm year −2 ) (Figure 4b). Approximately 16.7% of the area of the MRB had significant negative ET a trends, and most of this area corresponded to mixed forest (with an average decreasing rate of 1.15 mm year −2 ) and cropland/natural vegetation mosaic, because solar radiation and precipitation both showed decreasing trends in these regions (Appendix B, Figure A1).  The spatial pattern of annual available water resources (AWR, the difference between precipitation and ET a ) variation is highly consistent with that of precipitation, as shown in Figure 5 and Appendix B, Figure A1. From 1981 to 2012, the annual AWR in approximately 23.6% of the grids throughout the MRB showed significant positive trends, whereas 6.4% of the grids showed significantly negative trends. The grids with positive trends were distributed mainly in the evergreen broadleaf forest of the central reach of the basin and the savanna of the southeastern basin ( Figure 5). For the entire basin, the average increasing rate of AWR was about 0.32% year −1 . The largest increasing rate was found in evergreen broadleaf forest (approximately 0.88% year −1 ), followed by that in savanna (0.56% year −1 ), cropland and natural vegetation mosaic (0.41% year −1 ) and grassland (0.23% year −1 ). The grids with negative trends were distributed primarily in the mixed forest of the upper basin and the cropland of the Mekong Delta. The precipitation decreased significantly in the mixed forest, which resulted in a decrease in the AWR ( Figure 5). Although the precipitation showed a minor increasing trend in the paddy fields of the Mekong Delta (Appendix B, Figure A1), the ET a increased sharply (Figure 4), resulting in a significant decrease in the AWR. This finding implies that additional water is needed for agricultural development in the Mekong Delta. For the entire basin, the correlation coefficient between the AWR anomaly and precipitation anomaly was 0.96, which is much higher than that between the AWR anomaly and ET a anomaly (0.14) (Figure 6), demonstrating that water consumption has a minor effect on the regional variation of AWR. Thus, the inter-annual variation of AWR is affected mainly by the water supply (precipitation).

Attribution of Actual ET Change to Climate Change and Vegetation Greening
The inter-annual variation of ET a was mainly influenced by changes in precipitation, solar radiation, temperature, vapor pressure, atmospheric CO 2 and the LAI (Figure 7). Although increasing atmospheric CO 2 concentrations reduced stomatal conductance, which could have been responsible for the decreased ET a , the effect of CO 2 fertilization on ET a was very weak (0.0018% year −1 ), indicating that climate change and vegetation greening were the dominant factors driving the terrestrial ET a variation.
From 1981 to 2012, except for net radiation, all the climate factors showed increasing trends in the MRB. The air temperature showed a significant increasing trend in all land-use types with an average increasing rate of 0.23 • C per decade. Accompanied by precipitation increasing (Appendix B, Figure A1), a warming-wetting trend was found throughout the MRB, which is beneficial to vegetation growth in most parts of the MRB. Consequently, the LAI showed an increasing trend in all land-use types except mixed forest (−0.17% year −1 ). The highest LAI increasing rate was found in farmland (1.42% year −1 ), followed by grassland (0.96% year −1 ), savanna (0.68% year −1 ) and evergreen broadleaf forest (0.49% year −1 ). Increases in precipitation, temperature and the LAI had positive effects on the ET a trend that counteracted the negative effects caused by decreases in solar radiation and increases in vapor pressure. For the entire basin, the relative contributions of temperature, precipitation, solar radiation, vapor pressure and LAI to the ET a trend were 43.7%, 21.8%, −16.8%, −2.8% and 54.1%, respectively (Figure 7). Precipitation was the dominant driving factor affecting ET a in approximately 42% of the grids in the MRB, and most of these grids were located in grassland, evergreen broadleaf forest and savanna areas. The LAI had a significant positive effect on the ET a trend in the grassland of the upper basin, the savanna of the central basin and the cropland of the Mekong Delta. Solar radiation and temperature were the dominant factors driving the ET a variation in approximately 35% of the grids throughout the MRB, and most of these grids were located in cropland, grassland and savanna (Figure 8).  Variations in actual ET are principally controlled by variations in the water supply and AED. McVicar [47] divided the environment into three types: energy-limited areas (ET p /p < 0.76), "equitant" areas (0.76 < ET p /p < 1.35), and water-limited areas (ETp/p > 1.35). In the grassland of the upper basin (ET p /p = 1.37), precipitation is the dominant factor affecting the observed ET a variation, but increasing temperatures also play an important role (Figure 9). Glacial meltwater is an important water source in the upper basin. Temperature increasing not only improved the thermal conditions for vegetation growth but also accelerated glacial melting, which may provide more water for vegetation growth. The negative effect of ET a on runoff was offset by the positive effect of glacial meltwater in the Lancang River [48]. Therefore, although the grassland in the upper basin is a water-limited region (precipitation and the LAI are the principal driving factors in approximately 68% of grassland areas), the relative contribution of temperature varied from 5% to about 25% in more than 65% of the grassland grids (Figure 9).
In the central and lower reaches of the basin, the positive effect of increasing temperature on ET a was offset by the negative effect of reduced solar radiation, especially in evergreen broadleaf forest (ET p /p = 0.68), mixed forest (ET p /p = 0.83) and savanna (ET p /p = 0.76) (Figure 9). In more than 75% of grids of the abovementioned vegetation types, the relative contributions of precipitation varied from 25% to about 70%, while the relative contribution of temperature and solar radiation varied from 5% to about 65% and from −25% to about 5%, demonstrating that the water supply (precipitation) and energy (temperature and solar radiation) have equivalent impacts on the ET a variation. Because irrigation is a major source of water in cropland (ET p /p = 0.73) and cropland/natural vegetation mosaic (ET p /p = 0.69), the relative contribution of precipitation to ET a variation is weak (within a range of ±15% in more than 85% of the grids). Thus, temperature and solar radiation are the principal driving factors in more than 80% of the grids, indicating that cropland and cropland/natural vegetation mosaics are energy-limited regions.

Contribution of Land-Use Changes to ET a
Since 2000, the MRB has been characterized by the expansion of cropland and cropland/natural vegetation mosaic. In the central basin, the increase in cropland was mainly attributed to the conversion of cropland/natural vegetation mosaic, and the loss of cropland/natural vegetation mosaic was compensated by the conversion of savanna, which resulted in a decrease in savanna (Table 1). In the Mekong Delta, the increase in cropland benefited from the conversion of wetlands and water bodies. Another noticeable change is that approximately 18.55% and 16.38% of mixed forest were converted into evergreen broadleaf forest and savanna respectively, which resulted in a decrease in mixed forest.
For the entire basin, land-use change enhanced the increase of ET a by 0.014% year −1 , 53.24% and 29.34% of which were attributed to the conversion of mixed forest and savanna respectively, into other land-use types ( Figure 10). Compared with the effects of climate change and vegetation greening on ET a variation (0.15% year −1 ), the relative contribution of land-use change on ET a variation was only 9.3%, indicating that climate change and vegetation greening are the principal factors affecting the regional ET a variation. However, the expansion of cropland in the Mekong Delta enhanced the ET a increase by 15%, demonstrating that land-use change notably affects the water balance in this region.

Impact of Climate Change and Vegetation Greening on ET a
In the upper MRB, ET a is limited by the soil moisture supply. Although the increased rate of precipitation in the upper basin is lower than that in the lower basin ( Figure A1), ET a is sensitive to small changes in precipitation in water-limited regions. In equitant regions, such as evergreen broadleaf forest, the annual precipitation is larger than that in water-limited regions. ET a is not totally limited by the soil moisture supply; thus, the sensitivity and relative contribution of precipitation decline. Similarly, because soil moisture in the cropland of the Mekong Delta is supplied by both precipitation and irrigation (Figure 9), ET a is insensitive to variations in precipitation, while variations in air temperature and solar radiation constitute the predominant driving forces for potential ET, and hence ET a . Therefore, precipitation is responsible for the observed variation in ET a over arid and semiarid regions, whereas air temperature and incident solar radiation are more important over humid regions [6,49,50].
For the entire basin, vegetation greening provides a greater contribution (54.1%) to the variation in ET a than climate change does. Under enhanced CO 2 concentrations, leaf stomatal conductance will decrease, and increasing atmospheric CO 2 concentrations could be responsible for decreases in ET a . Although the direct effect of CO 2 fertilization on ET a is weak (0.02 mm year −1 in MRB and 0.05 mm year −1 at the global scale) [4], it must be noted that CO 2 fertilization explains 70% of the global observed greening trend [51], and its negative effect may be entirely compensated by increases in leaf area. From 1981 to 2011, approximately 55% ± 25% of the observed global terrestrial ET a increase was caused by vegetation greening [52]. The same phenomenon was also found in MRB, increases in the LAI are associated with the intensification of terrestrial ET a (Figure 7). Considering that the greening of the Earth is projected to continue worldwide during the 21st century in most Coupled Model Intercomparison Project Phase 5 (CMIP5) Earth system models involved in the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) [53], the water cycle will continue to be accelerated by greening-induced variations in ET a .

Variation in Available Water Resources
As the Mekong River crosses several international boundaries, the water resources of the downstream countries are strongly connected with those of the upstream countries. The upper basin flows are affected mainly by snowmelt and precipitation [17], the negative effect of ET a on runoff was offset by the positive effect of glacial meltwater in the upper basin (Lancang River in China) [54]. Because the AWR are calculated herein as the difference between precipitation and ET a , which ignores the contribution from snowmelt, the AWR in the upper basin were found to exhibit a minor increasing trend ( Figure 5). Although the AWR increased in the upper basin, the maximum flows at the Chiang Saen gauging station, which is located in the upper MRB, showed a decreasing trend, and the decreasing rate of flow was accelerated due to the completion of dams in the upstream region [55]. Other studies have further argued that climate change is not the main mechanism responsible for the variation in AWR in the Mekong River; instead, the development of dams on the Lancang River may be the major influencing factor, especially in the dry season and upstream region [55,56]. Therefore, although the AWR increased in the upper basin, the conditions of the actual water resources both upstream and downstream remain somewhat unclear.
In the central basin and Mekong Delta, cropland expansion has placed great pressure on the AWR; for example, a dramatic increase in ET a and a significant decrease in AWR were found in the Mekong Delta and southeastern Thailand (Figures 4 and 5). In the Mekong Delta, the regional population grew by nearly 45% from 1980 to 2000 [57]. Additionally, rice cropping systems have undergone remarkable changes from single-cropping to doubleor triple-cropping systems [58]. Rapid population growth and agricultural development have placed additional pressures on cropland and water for food production. In Lao, Thailand, Cambodia and the central highlands of Vietnam, rain-fed rice is the primary crop [25], and irrigation is intended mainly for rice during the rainy season due to its low water requirements [59]. Due to changes in the climate and dam construction, changes in seasonal flow patterns have strongly influenced land-use patterns [1]. Fortunately, a significant increase in dry-season flows and a decrease in wet-season flows occurred when the Xiaowan dam (constructed in 2010) and Nuozhadu dam (constructed in 2014) were erected in the upper reaches of the Mekong River (Lancang River) [54]. Increased dryseason flows may benefit agricultural irrigation downstream, and decreased wet-season flows may benefit downstream flood and drought management.

Conclusions
The spatial-temporal patterns of ET a and AWR in the MRB were estimated from 1981 to 2012 using the VIP-RS model. The contributions of climate change and vegetation greening to the ET a trends were quantitatively determined with 30 model experiments. Due to vegetation greening and increasing temperature and precipitation, ET a showed an increasing trend at a rate of 1.16 mm year −2 from 1981 to 2012. In most of the basin, the water and energy had equivalent impacts on the ET a trends. Although the AED increased, the AWR showed an increasing trend due to increasing precipitation. In the paddy fields of the Mekong Delta region, the ET a increased significantly, resulting in additional pressure on agricultural development. In the central basin where rain-fed rice is planted, the AWR is decreasing following the reduced precipitation. Therefore, the MRB may face water shortages in the dry season.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

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

Appendix A
The air temperature was corrected to sea level as follows: where T sl ( • C) is the air temperature at sea level, and T ( • C) is the air temperature at real elevation. If the precipitation is higher than 3 mm day −1 , it was corrected to sea level as follows. If the precipitation is lower than 3 mm day −1 , it was directly resampled to 5 × 5 km with the bilinear resampling method.
where P sl (mm day −1 ) is the precipitation at sea level, and P (mm/day) is the precipitation at real elevation. The wind speed was corrected to the value at the 2 m level as follows: where W sl (m s −1 ) is the wind speed at 2 m, and W (m s −1 ) is the wind speed at real elevation.