**1. 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<sup>2</sup> concentrations can also induce

**Citation:** Hu, S.; Mo, X. Attribution of Long-Term Evapotranspiration Trends in the Mekong River Basin with a Remote Sensing-Based Process Model. *Remote Sens.* **2021**, *13*, 303. https://doi.org/10.3390/rs13020303

Received: 1 December 2020 Accepted: 9 January 2021 Published: 16 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

changes in ET by stimulating the photosynthetic rate, thereby increasing plant growth and leaf areas, which accelerates the rate of ET [7]. Due to increasing CO<sup>2</sup> fertilization, nitrogen deposition, climate change and land-use change, the global leaf area index (LAI) increased significantly by 8% during 1982–2011 [8], which had a significant impact on the intensification of terrestrial ET. Numerical modeling has played an important role in simulating terrestrial ET [9,10], and provides a possible solution for quantifying the contributions of different environmental factors on ET [9]. For example, with a land surface model, Shi et al. [10] and Mao et al. [11] pointed out that spatio-temporal variations in ET were determined mainly by climate variables. However, when a remote-sensing ET model with observation-based environmental drivers was applied, it was found that more than half of the accelerated global terrestrial ET was caused by vegetation greening [9], which suggests that the attribution of regional terrestrial ET trends without considering vegetation dynamics is improper. Considering the climate changes and vegetation greening on regional ET and available water supply are still unclear, based on the remote sensed LAI data, the contributions of climate control and vegetation greening to ET<sup>a</sup> can be quantified with an ecological model.

The Mekong River Basin (MRB), which is rich in natural resources, spans large parts of continental Southeast Asia. However, under the impacts of climate change and rapid socioeconomic development over recent decades, the natural resources within the MRB have faced increasing pressure. Since the 1950s, the MRB has been characterized by significant increases in the temperature and spatial-temporal variation of precipitation [12,13]. The hydrological cycle has therefore intensified, resulting in changes in the rates and patterns of ET and river flows [14,15]. Differences in flow variability were found between the upper and lower MRB because floods mostly originate from rainfall and snowmelt in the upper MRB [16], while the hydrology of the lower Mekong is controlled by the western North Pacific monsoon [17]. Moreover, a significant positive correlation was found between discharge, precipitation and snow cover in the upper MRB, especially during the dry season [18]. In contrast, in the lower MRB, changes in discharge are driven primarily by the variation in the Asian summer monsoon [17]. In addition, rapid economic development, population growth and widespread agricultural deforestation have created pressures leading to competition for water [19]. Land-use changes, which manifest mainly as transformation into cropland, have increased the consumption of water through irrigation [20,21] and affected the land surface roughness and albedo, thereby altering both the hydrologic processes [22] and the energy balance at the surface [23,24]. Therefore, although the average variations in ET and runoff at the basin scale are relatively small, the changes in these parameters in the highly irrigated areas of the lower MRB are significant due to the expansion of cropland and additional irrigation [20,25].

According to the results of climate modeling, the monsoon patterns in the MRB will change, resulting in significant increases in the average precipitation and temperature over the basin [26]. However, it remains unclear how these climate changes and the overlapping effects of land-use change will alter the regional ET and available water supply. In this study, the long-term actual ET (ETa) in the MRB was simulated by Remote Sensing-Based Vegetation Interface Processes Model, and we aim to (1) explore the trends and variability of ET<sup>a</sup> and the availability of water resources in recent decades and (2) quantify the contributions of climate control and vegetation greening to ETa. This paper is organized as follows, Section 2 describes the datasets and attribution methods using VIP-RS model, while Section 3 highlights the results for the ET<sup>a</sup> variation and identifies relative contributions of climate variables and vegetation greening to ETa. Finally, Sections 4 and 5 delve into the discussions and conclusions of this work.

#### **2. Materials and Methods**

#### *2.1. Study Region*

As the longest river in Southeast Asia, the Mekong River (4800 km) passes through six countries and boasts a drainage area of 805,604 km<sup>2</sup> . As a result of the Asian monsoon, 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).

**Figure 1.** (**a**) Location, (**b**) elevation, (**c**) land-use/cover type and (**d**) soil texture in the Mekong River Basin (MRB) (WB: Water body, EBF: Evergreen broadleaf forest, MF: Mixed forest, SV: Savanna, GL: Grassland, WL: Wetland, CL: Cropland, Urb: Urban and built-up, CL/NV: Cropland and natural vegetation mosaic).

#### *2.2. Data*

#### 2.2.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).

#### 2.2.2. 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.

#### 2.2.3. 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<sup>a</sup> 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 ETa, monthly ET<sup>a</sup> data calculated from the water balance equation were also used to validate the ET<sup>a</sup> results at the basin scale:

$$ET\_a = \left. P \right| - \left. R - \frac{ds}{dt} \right. \tag{1}$$

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.

The runoff data were provided by the Global Runoff Data Centre (https://www.bafg. de/GRDC/EN/02\_srvcs/21\_tmsrs/riverdischarge\_node.html).

#### *2.3. Method*

#### 2.3.1. Model Introduction

ET<sup>a</sup> was simulated by the Remote Sensing-Based Vegetation Interface Processes (VIP-RS) model [33,34]. The total ET<sup>a</sup> is composed of evaporation due to canopy interception (*E<sup>i</sup>* ), transpiration from vegetation (*Ec*) and soil evaporation (*Es*). Transpiration

(*Ec*) is estimated based on potential transpiration (*Ecp*) and is limited by water conditions (*fw*) and temperature (*ft*), as follows:

$$E\_{\mathfrak{c}} = E\_{\mathfrak{cp}} f\_{\mathfrak{w}} f\_{\mathfrak{t}} \tag{2}$$

where *f<sup>w</sup>* is the stress function of the atmospheric water vapor pressure deficit calculated using the algorithms proposed by Mu et al. [35] and *f<sup>t</sup>* is the stress function of the air temperature calculated using the algorithms proposed by Zhang et al. [36]. *Ecp* is calculated using the Penman-Monteith (PM) equation as follows:

$$E\_{cp} = \frac{1}{\lambda} \left(\Delta R\_{nc} + f\_c \rho c\_p D / r\_a\right) / \left(\Delta + \gamma \eta f\_{CO\_2}\right) \tag{3}$$

where *Rnc* is the net radiation absorbed by the canopy, ∆ represents the slope of the saturated vapor pressure curve versus air temperature,*γ*, *ρ*, *c<sup>p</sup>* and *D* represent the psychrometric constant (hPa ◦C −1 ), air density (kg m−<sup>3</sup> ), specific heat capacity of air (J kg−<sup>1</sup> ◦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−<sup>1</sup> ) and *r<sup>a</sup>* is the aerodynamic resistance between the canopy and the reference height (s m−<sup>1</sup> ). *f<sup>c</sup>* is the fractional vegetation cover, and *fCO*<sup>2</sup> is the stress function of atmospheric CO<sup>2</sup> [37]. These two factors are calculated as follows:

$$f\_{\rm CO\_2} = \frac{1}{(-0.001 \text{CO}\_2 + 1.35)}\tag{4}$$

$$f\_c = \left(1 - \left(\frac{VI\_{\text{max}} - VI}{VI\_{\text{max}} - VI\_{\text{min}}}\right)^{\beta}\right. \tag{5}$$

$$VI = VI\_{\rm min} + \frac{f\_{\rm PAR} \left( VI\_{\rm max} - VI\_{\rm min} \right)}{f\_{\rm PAR\_{\rm max}} - f\_{\rm PAR\_{\rm min}}} \tag{6}$$

$$f\_{\text{PAR}} = 1 - e^{\left[\frac{LAl}{LAl\text{max}}\ln(1 - f\_{\text{PAR}\_{\text{max}}})\right]} \tag{7}$$

where CO<sup>2</sup> is the atmospheric CO<sup>2</sup> concentration (ppm), *β* is an empirical constant, being 0.6 to 1.2, taken here as 0.7 [38], and *VImax* and *VImin* represent the vegetation index values under full vegetation cover and bare soil conditions, respectively (set to 0.95 and 0.01). Additionally, *fPARmax* and *fPARmax* are the maximum and minimum *fPAR* corresponding to *VImax* and *VImin* respectively, and are set to 0.95 and 0.001. *LAImax* 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<sup>i</sup>* ) equals the potential evaporation from the wetted surface. Soil evaporation (*Es*) is the minimum value of surface potential evaporation (*Esp*) and soil moisture exfiltration (*Eex*) [40]:

$$E\_s = \min\left(E\_{sp\nu}, E\_{\text{ex}}\right) \tag{8}$$

$$E\_{\rm sp} = \frac{1}{\lambda} \left( \Delta (R\_{\rm ns} - G) + (1 - f\_c) \rho c\_p D / r\_{\rm as} \right) / (\Delta + \gamma) \tag{9}$$

$$\mathcal{G} = \mathcal{R}\_{\mathfrak{m}} \times \left[ \Gamma\_{\mathfrak{c}} + (1 - f\_{\mathfrak{c}}) \times (\Gamma\_{\mathfrak{s}} - \Gamma\_{\mathfrak{c}}) \right] \tag{10}$$

where *Rns* is the net radiation absorbed by the soil surface, *G* is the soil heat flux (MJ d−<sup>1</sup> ) [41], *ras* is the aerodynamic resistance between the reference height and the soil surface (s m−<sup>1</sup> ), Γ*<sup>c</sup>* (set as 0.05) [42] and Γ*<sup>s</sup>* (set as 0.315) [43] are the ratios of *G* to *Rnc* for full vegetation canopy and bare soil, respectively.

#### 2.3.2. 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<sup>a</sup> 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<sup>a</sup> from 1981 to 2012 were used for the subsequent analysis. The tendency of annual ET<sup>a</sup> and its significances were explored based on the Theil-Sen estimator and Mann–Kendall (M-K) test.

The simulated monthly ET<sup>a</sup> was consistent with the actual ET from GLEAM (Figure 2a) and yielded an R<sup>2</sup> value of 0.84. The mean bias and root mean square error (RMSE) between the simulated and GLEAM ET<sup>a</sup> were 4.5 mm month−<sup>1</sup> and 8.17 mm month−<sup>1</sup> , respectively. At annual scale, the GRACE-derived ET<sup>a</sup> and simulated ET<sup>a</sup> showed acceptable consistency with an R<sup>2</sup> of 0.55. However, the GRACE-derived ET<sup>a</sup> was 25.72 mm year−<sup>1</sup> lower than the simulated ETa, indicating that the VIP-RS model slightly overestimated the yearly ET<sup>a</sup> in the MRB.

**Figure 2.** Comparisons of the simulated ETa (actual evapotranspiration) with (**a**) the GLEAM ETa and (**b**) water balancederived ETa.

#### 2.3.3. Attribution of ET<sup>a</sup> 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<sup>a</sup> 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<sup>2</sup> concentration and vegetation dynamics from 1982 to 1990. The ET<sup>a</sup> 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<sup>2</sup> concentration is the average CO<sup>2</sup> 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 ETa, the precipitation in each year (1982–1990) is used in the control simulation, and the final ET<sup>a</sup> 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<sup>2</sup> concentration.

The main effect of each factor on the ET<sup>a</sup> trend is the difference between the ET<sup>a</sup> 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<sup>a</sup> trend are as follows:

$$E(var1) = f(var1) - f(con) \tag{11}$$

$$E(var2) = f(var2) - f(con) \tag{12}$$

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:

$$E(var1\\_var2) = f(var1\\_var2) - f(con) - E(var1) - E(var2) \tag{13}$$

where *E(var1)* and *E(var2)* are the effects of *var1* and *var2* on the ET<sup>a</sup> 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 ETa, 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<sup>a</sup> trend is the difference between *f*(*land\_type1*) and *f*(*land\_con*). Considering that ET<sup>a</sup> 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 ETa, the ET<sup>a</sup> 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<sup>a</sup> simulation:

$$ET\_d = \sum\_{i=1,9} (ET\_{ai} \times V\_i) \tag{14}$$

where *ETai* is the ET<sup>a</sup> value in each sub-grid (5 km resolution) with land-use type *i* (Figure 1) and *V<sup>i</sup>* is the proportion of land-use type *i* in each sub-grid with a 5 km resolution. *V<sup>i</sup>* is obtained from land-use data with a spatial resolution of 500 m.

#### **3. Results**

#### *3.1. Spatial-Temporal Variation of Actual ET*

ET<sup>a</sup> is principally controlled by the AED (represented by ETp) and precipitation. The annual ET<sup>a</sup> in the MRB showed remarkable spatial heterogeneity with an average value of 896 <sup>±</sup> 34 mm year−<sup>1</sup> and increased gradually from north to south with a range of 300 to 2000 mm year−<sup>1</sup> . From 1982 to 2012, 70% of the study area showed an increasing ET<sup>a</sup> trend with an average increase rate of 1.16 mm year−<sup>2</sup> .

In the upper basin, where the precipitation is less than 800 mm year−<sup>1</sup> and the AED is weak, the annual ET<sup>a</sup> is less than 400 mm year−<sup>1</sup> (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<sup>a</sup> 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−<sup>1</sup> ),

the spatial pattern of the average annual ET<sup>a</sup> is highly consistent with the average annual ET<sup>p</sup> (Figure 3), whereas in cropland/natural vegetation mosaic areas, which have an average annual precipitation of less than 1200 mm year−<sup>1</sup> , approximately 75% of the rainfall returns to the atmosphere through evaporation; as a result, the spatial pattern of ET<sup>a</sup> is similar to that of precipitation. A significant positive temporal trend of ET<sup>a</sup> was found in the Mekong Delta and the southeastern basin (Figure 4a), where precipitation and LAI display increasing trends. Significant positive trends of ET<sup>a</sup> 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<sup>a</sup> (Figure 4c). The highest increasing rate was found in cropland (2.53 mm year−<sup>2</sup> ), followed by grassland (1.8 mm year−<sup>2</sup> ), savanna (1.68 mm year−<sup>2</sup> ) and evergreen broadleaf forest (1.25 mm year−<sup>2</sup> ) (Figure 4b). Approximately 16.7% of the area of the MRB had significant negative ET<sup>a</sup> trends, and most of this area corresponded to mixed forest (with an average decreasing rate of 1.15 mm year−<sup>2</sup> ) and cropland/natural vegetation mosaic, because solar radiation and precipitation both showed decreasing trends in these regions (Appendix B, Figure A1).

**Figure 3.** Spatial patterns of annual average (**a**) actual evapotranspiration, (**b**) precipitation and (**c**) potential evapotranspiration over Mekong River Basin for the period of 1981 to 2012.

**Figure 4.** (**a**) Spatial patterns of the temporal variations of average actual evapotranspiration, (**b**) average change rate of actual evapotranspiration in different vegetation types, and (**c**) probability distribution function of significant levels in the actual evapotranspiration trend.

The spatial pattern of annual available water resources (AWR, the difference between precipitation and ETa) 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−<sup>1</sup> . The largest increasing rate was found in evergreen broadleaf forest (approximately 0.88% year−<sup>1</sup> ), followed by that in savanna (0.56% year−<sup>1</sup> ), cropland and natural vegetation mosaic (0.41% year−<sup>1</sup> ) and grassland (0.23% year−<sup>1</sup> ). 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<sup>a</sup> 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<sup>a</sup> 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).

**Figure 5.** Spatial patterns of the temporal variation of annual average available water resources (AWR) for the period of 1981 to 2012.

**Figure 6.** Inter-annual variations of anomalies of (**a**) the actual ET (ETa) anomaly, potential ET (ETp) anomaly and precipitation (P) anomaly, and (**b**) available water resources (AWR) anomaly from 1981 to 2012.

#### *3.2. Attribution of Actual ET*

3.2.1. Attribution of Actual ET Change to Climate Change and Vegetation Greening

The inter-annual variation of ET<sup>a</sup> was mainly influenced by changes in precipitation, solar radiation, temperature, vapor pressure, atmospheric CO<sup>2</sup> and the LAI (Figure 7). Although increasing atmospheric CO<sup>2</sup> concentrations reduced stomatal conductance, which could have been responsible for the decreased ETa, the effect of CO<sup>2</sup> fertilization on ET<sup>a</sup> was very weak (0.0018% year−<sup>1</sup> ), indicating that climate change and vegetation greening were the dominant factors driving the terrestrial ET<sup>a</sup> 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−<sup>1</sup> ). The highest LAI increasing rate was found in farmland (1.42% year−<sup>1</sup> ), followed by grassland (0.96% year−<sup>1</sup> ), savanna (0.68% year−<sup>1</sup> ) and evergreen broadleaf forest (0.49% year−<sup>1</sup> ). Increases in precipitation, temperature and the LAI had positive effects on the ET<sup>a</sup> 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<sup>a</sup> trend were 43.7%, 21.8%, −16.8%, −2.8% and 54.1%, respectively (Figure 7). Precipitation was the dominant driving factor affecting ET<sup>a</sup> 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<sup>a</sup> 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<sup>a</sup> variation in approximately 35% of the grids throughout the MRB, and most of these grids were located in cropland, grassland and savanna (Figure 8).

**Figure 7.** Contributions of changes in climate variables, the atmospheric CO<sup>2</sup> concentration and the leaf area index (LAI) to the ETa variation and their two-factor interactive effects (Net: total effect of all the variables, C: atmospheric CO<sup>2</sup> concentration, T: temperature, P: precipitation, R: radiation, H: vapor pressure, W: wind speed, Sum\_EF: sum of all the two-factor interactive effects of variables).

**Figure 8.** Spatial patterns of the relative contributions of the variation of (**a**) LAI, (**b**) temperature, (**c**) precipitation, (**d**) radiation, (**e**) vapor pressure, (**f**) wind speed, (**g**) atmospheric CO<sup>2</sup> concentration variation to ETa and (**h**) probability distribution function of relative contributions of the main climatic variables and LAI changes to ETa trends in Mekong River Basin.

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 (ETp/*p* < 0.76), "equitant" areas (0.76 < ETp/*p* < 1.35), and water-limited areas (ETp/*p* > 1.35). In the grassland of the upper basin (ETp/*p* = 1.37), precipitation is the dominant factor affecting the observed ET<sup>a</sup> 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<sup>a</sup> 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<sup>a</sup> was offset by the negative effect of reduced solar radiation, especially in evergreen broadleaf forest (ETp/*p* = 0.68), mixed forest (ETp/*p* = 0.83) and savanna (ETp/*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<sup>a</sup> variation. Because irrigation is a major source of water in cropland (ETp/*p* = 0.73) and cropland/natural vegetation mosaic (ETp/*p* = 0.69), the relative contribution of precipitation to ET<sup>a</sup> 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.

**Figure 9.** Probability distribution function of relative contributions of climatic variables and LAI changes to ETa trends in different vegetation types.

#### 3.2.2. Contribution of Land-Use Changes to ET<sup>a</sup>

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<sup>a</sup> by 0.014% year−<sup>1</sup> , 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<sup>a</sup> variation (0.15% year−<sup>1</sup> ), the relative contribution of land-use change on ET<sup>a</sup> variation was only 9.3%, indicating that climate change and vegetation greening are the principal factors affecting the regional ET<sup>a</sup> variation. However, the expansion of cropland in the Mekong Delta enhanced the ET<sup>a</sup> increase by 15%, demonstrating that land-use change notably affects the water balance in this region.


**Table 1.** Land-use transformation matrix (%).

WB: Water body, EBF: Evergreen broadleaf forest, MF: Mixed forest, SV: Savanna, WL: Wetland, GL: Grassland, CL: Cropland, Urb: Urban and built-up, CL/NV: Cropland and natural vegetation mosaic.

**Figure 10.** Effects of land-use changes on the ETa trend (Net: total effect of all the land-use changes).

#### **4. Discussion**

#### *4.1. Impact of Climate Change and Vegetation Greening on ET<sup>a</sup>*

In the upper MRB, ET<sup>a</sup> 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<sup>a</sup> 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<sup>a</sup> 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<sup>a</sup> is insensitive to variations in precipitation, while variations in air temperature and solar radiation constitute the predominant driving forces for potential ET, and hence ETa. Therefore, precipitation is responsible for the observed variation in ET<sup>a</sup> 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<sup>a</sup> than climate change does. Under enhanced CO<sup>2</sup> concentrations, leaf stomatal conductance will decrease, and increasing atmospheric CO<sup>2</sup> concentrations could be responsible for decreases in ETa. Although the direct effect of CO<sup>2</sup> fertilization on ET<sup>a</sup> is weak (0.02 mm year−<sup>1</sup> in MRB and 0.05 mm year−<sup>1</sup> at the global scale) [4], it must be noted that CO<sup>2</sup> 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<sup>a</sup> 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<sup>a</sup> (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 ETa.

#### *4.2. 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<sup>a</sup> 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 ETa, 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<sup>a</sup> 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.

#### **5. Conclusions**

The spatial-temporal patterns of ET<sup>a</sup> 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<sup>a</sup> trends were quantitatively determined with 30 model experiments. Due to vegetation greening and increasing temperature and precipitation, ET<sup>a</sup> showed an increasing trend at a rate of 1.16 mm year−<sup>2</sup> from 1981 to 2012. In most of the basin, the water and energy had equivalent impacts on the ET<sup>a</sup> 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<sup>a</sup> 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.

**Author Contributions:** Conceptualization, X.M. and S.H.; methodology, X.M.; writing—original draft preparation, S.H.; writing—review and editing, X.M. and S.H.; project administration, X.M.; funding acquisition, X.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Key R&D Program of China, grant number 2017YFA0603702, and National Natural Science Foundation of China, grant number 41971232.

**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:

$$T\_{sl} = T + 22 + \frac{0.68 \times (\text{DEM} - 5000)}{100} \left(\text{DEM} > 5000\right) \tag{A1}$$

$$T\_{sl} = T + \frac{0.44 \times \text{DEM}}{100} \text{ (DEM } \le 5000) \tag{A2}$$

where *Tsl* ( ◦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−<sup>1</sup> , it was corrected to sea level as follows. If the precipitation is lower than 3 mm day−<sup>1</sup> , it was directly resampled to 5 × 5 km with the bilinear resampling method.

$$P\_{s\,\,l} = P - \frac{0.055 \times \text{DEN}}{100} \,\, (P \,\, \ge \,\, 3) \tag{A3}$$

where *Psl* (mm day−<sup>1</sup> ) 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:

$$\mathcal{W}\_{\rm sl} = \,^\text{W} \times \frac{4.87}{\log\_{10}(67.8 \times \text{DEM} - 5.42)}\tag{A4}$$

where *Wsl* (m s−<sup>1</sup> ) is the wind speed at 2 m, and *W* (m s−<sup>1</sup> ) is the wind speed at real elevation.

#### **Appendix B**

**Table A1.** Summary of simulation experiments.


#### **Table A1.** *Cont.*


**Figure A1.** Spatial patterns of the temporal variations of the annual (**a**) leaf area index (LAI), (**b**) precipitation, (**c**) temperature and (**d**) radiation from 1981 to 1012.

#### **References**

