Hydrological responses to the synergistic climate and land-use changes in the upper Lancang river basin

As the largest transboundary river in Asia, the Lancang-Mekong River originates from the Qinghai-Tibetan Plateau and traverses six riparian countries including China, Laos, Thailand, Myanmar, Cambodia, and Vietnam. The focus on this paper is the upper LRB, which constitutes an area of 52 508 km<sup="">2</sup> and serves as the source area of the LMR. The upper LRB is characterized by high altitudes ranging from 3000 to 6000 meters and topographic ruggedness. The basin information are obtained from multiple, openly accessible datasets (table S1).

Detailed climate information encompassing multiple variables is imperative to evaluate the influence of climate change on watershed hydrology. With the release of the Coupled Model Intercomparison Project Phase 6 (CMIP6), three widely adopted General Circulation Models (GCMs): GFDL-CM4.0 [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib25" id="fnref-erlad1347bib25"="">25</a>], MIP-ESM1.2-LR [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib26" id="fnref-erlad1347bib26"="">26</a>], and MIROC6 [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib27" id="fnref-erlad1347bib27"="">27</a>], are employed for climate change projections in the Tibetan Plateau.

The GCMs generate daily climate data (precipitation, maximum temperature, and minimum temperature) for the upper LRB across various horizontal resolutions under two emission scenarios: Shared Socioeconomic Pathways (SSP) 2.45 and SSP 585. These two scenarios respectively represent moderate and high levels of radiative forcing, reflecting future climate fluctuations under different intensities of human activities. To reduce the contingencies in basin-scale impact assessment, we propose a generalized joint bias correction (JBC) method (Note S1) that couples the Quantile Mapping method [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib28" id="fnref-erlad1347bib28"="">28</a>] and the JBC method [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib29" id="fnref-erlad1347bib29"="">29</a>] to correct the simulation datasets derived from GCMs projections. The corrected grid-scale corrected datasets are further subjected to statistical downscaling using the Empirical Bayesian Kriging geographic interpolation approach [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib30" id="fnref-erlad1347bib30"="">30</a>], reducing the resolution to 0.5° <b="">&#x00d7;</b> 0.5° and facilitating the interpolation to meteorological stations in the upper LRB.

The pixel-wise MLR model allows for a quantitative analysis to investigate the relationship between climate changes and vegetation cover dynamics at the pixel level in the upper LRB. In this study, the normalized difference vegetation index (NDVI) derived from satellite data is chosen as an indicator of live green vegetation growth in the high-mountain area, which is valuable in understanding the ongoing vegetation dynamics and the potential implications for land use. The terrestrial vegetation dynamics is influenced by both long-term changes in climate averages and extreme climate events. Here, five representative climatic variables (precipitation, maximum temperature, minimum temperature, RX1 day, and R10 mm) are employed as explanatory variables for NDVI estimation. R10 mm represents the number of days per year characterized by heavy precipitation (⩾10 mm d<sup="">−1</sup>); RX1 day describes the annual maximum 1 day precipitation. The MLR model incorporates downscaled climatic variables as inputs and employs NDVI images as target outputs during the period 2000–2015:

$$\begin{equation}\begin{array}{*{20}{c}} {{\text{NDV}}{{\text{I}}_{ix}} = {b_{0x}} + \sum\limits_{j = 1}^p {b_{jx}}{x_{ijx}} + {e_{ix}}} \end{array}.\end{equation} \tag{ 1 }$$

In equation (1), considering the observation $i{\text{ }} \in \left\{ {1, \ldots ,n} \right\}$ for each year, the variable ${x_{ijx}}$ denotes the observed value of predictor $j$ in pixel $x$, and ${\text{NDV}}{{\text{I}}_{ix}}$ represents the real-valued NDVI in pixel $x$. The parameters ${b_{0x}}$ and ${b_{jx}}$ refer to the intercept and regression slope associated with predictor $j{\text{ }}$in pixel $x$, respectively. Additionally, ${e_{ix}}$ represents a Gaussian error term.

To estimate the NDVI for the year 2050, we incorporate the predicted average climate variables within five-year intervals preceding and following 2050 (i.e. 2048–2052) as inputs for the MLR model. These average climate variables are utilized as estimations for the anticipated climate conditions around the year 2050 (figure S2).

The RF approach is a robust supervised machine learning algorithm used for non-linear classification and regression tasks. RF builds and collects multiple decision trees <em="">via</em> a bootstrap aggregating method on the training data to generate reliable predictions [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib31" id="fnref-erlad1347bib31"="">31</a>]. In this study, we developed a RF classifier for land use classification and projection in the upper LRB. Given the limited coverage of barren land, urban areas, and water bodies in the study area (which collectively accounted for approximately 5% of the total area), these land use types are grouped as 'others'. Therefore, the land use types considered in the analysis include forest, grass, cropland, and others. The RF classifier utilizes annual NDVI data from 2000 to 2015, along with static geographic information (i.e. elevation and slope) as input variables to simulate the land use maps from 2000 to 2015. Figure S3 presents boxplots illustrating various explanatory variables in the RF classifier segregated by each land use class, which emphasize the impacts of environmental factors on the classification in the study area.

In this study area, the observed land use data is class-imbalanced, with grassland accounting for over 70% of the total samples. To mitigate the potential bias in classification, an oversampling procedure is implemented to balance the priors of land use classes. To evaluate the experimental performance, we conduct a comparison between the observed land use map and the projected land use map for the year 2015 (figure S4). Additionally, a confusion matrix and compile statistics on classification performance metrics are created for the validation of the RF classification (tables S2 and S3). The RF classifier achieves with an overall accuracy of 0.88, with a 95% confidence level, suggesting a satisfactory performance of the algorithm in the study region. With the RF classifier, we project the land use map in 2050, using NDVI projected in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#erlad1347s2-3"="">2.3</a>, as well as elevation and slope, as input variables.

The SWAT model is a semi-distributed and physically based tool, which has been extensively used to tackle hydrological and environmental across various basin scales and geo-climatic conditions [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib32" id="fnref-erlad1347bib32"="">32</a>]. In this study, the calibrated SWAT model [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib33" id="fnref-erlad1347bib33"="">33</a>] was used for simulating the various scenarios of climate and land use in the region. In the upper LRB, the SWAT model delineated the study area into 150 sub-basins, which were further divided into 958 hydrological response units (HRUs) based on land use, soil type, and slope.

To estimate soil erosion and sediment yield within each spatial unit (i.e. Hydrological Response Unit)-HRU, the SWAT model employs a Modified Universal Soil Loss Equation (MUSLE) [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib34" id="fnref-erlad1347bib34"="">34</a>]. The MUSLE predicts annual soil loss as a result of six factors, including rainfall erosivity, soil erodibility, topography, land cover and management, and support practice. The rationale of the model for the study purpose and details of the model establishment can be seen in the Note S2.

The complex topography and atmospheric circulation dominated by Asian monsoons and westerly perturbations jointly influence the climate characteristics in the region of upper LRB. As a result, climate change trends in the region are especially prominent on a global scale [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib35" id="fnref-erlad1347bib35"="">35</a>]. Figures S5 and S6 illustrate the time series for climate variables throughout both the historical and future period, encompassing a range of climate projections. To examine the trends in climatic variables, the nonparametric Mann-Kendall (M-K) trend test (Note S3) was performed annually to the region, and the statistics results are shown in table S4. An overall increasing trend is found in the annual minimum and maximum temperature in the region. The existing change in the maximum and minimum temperature is 0.52 °C and 0.57 °C per decade during 1978–2020, respectively; the projected changes of the two variables are 0.69 °C and 0.51 °C per decade during 2021–2080, with a larger change in the maximum temperature than that of the minimum. Changes in precipitation, in terms of concentration, volume, and intensity, are also considerable. There was a non-significant downward trend in the annual precipitation during 1978–2020; while the projected annual precipitation during 2021–2080 shows a statistically upward trend. Furthermore, the frequency and intensity of extreme climate events had increased in the past and the trend will continue across much of the globe [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib36" id="fnref-erlad1347bib36"="">36</a>]. In this study, the moderate and extreme precipitation events are indicated using two indicators: R10 mm and RX1day on an annual scale throughout the study period (1978–2080). Results of the M-K test indicate a rise in RX1 day by 6.4 mm and R10 mm by 9.6 d over the future period of 2021–2080. The <em="">Z</em> values in the trend analysis of the two indicators are greater than 1.96, showing a prominent increase in the frequency of extreme climate variables at a high confidence level.

Land use in the upper LRB has undergone considerable transformations, in the forms of vegetation type alterations, land use transitions, and the spatial distribution of both vegetation and land use. Climate change was proved as the key factor driving the long-term change in vegetation cover [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib37" id="fnref-erlad1347bib37"="">37</a>]. In particular, the vegetation cover evolution under climate change has greatly impacted the spatial and temporal distribution of land use in this region. The results of a pixel-wise MLR model, as depicted in figures S7/8, show how vegetation cover dynamics are affected by the selected climate variables. The increase in NDVI can be attributed to favorable climatic conditions, specifically the rise in minimum temperature and precipitation. Using the MLR model and climate change projections, we estimate the NDVI in 2050 as well as its statistical distribution (figure S9).

These changes facilitate the necessary moisture and warmth required for vegetation growth, ultimately resulting in increased vegetation cover and land use change, along with socioeconomic drivers. Figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f1"="">1</a> shows the reclassified land use maps derived from the European Space Agency [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib38" id="fnref-erlad1347bib38"="">38</a>] for 1992(A) and 2015(B). The earliest (1992) and latest (2015) years of the historical period with data availability are chosen to represent the past and current land use conditions. The increase of cropland and decrease of grassland (as shown in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f1"="">1</a>(D)) reflect the joint impact of climate change, natural landscape, and the associated land use decisions over the past decades. The land use transition matrix depicted in table S5 delves into the land use transition processes between 1992 and 2015. In the study area, the primary transitions involved changes from grasslands to both forests and croplands in places with higher elevation and substantial transition from forest to cropland in places with lower elevation, with forests acting as an intermediary role, under both natural and socioeconomic driving forces [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlad1347bib39" id="fnref-erlad1347bib39"="">39</a>]. Thus, in an overarching view, statistical observation shows a primary loss with grassland and a primary gain with cropland, and a complicated change with forest (see more in the <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#erlad1347s4"="">discussion</a> section). The forest area increased in size at the end of last century due to more favorite climate condition but diminished after the year of 2000 due to dramatic human land development. The forest area increased in size at the end of last century due to more favorite climate condition, but diminished after the year of 2000 due to dramatic human land development.

Zoom In Zoom Out Reset image size

To assess the spatiotemporal dynamics in land use change and predict the future land use under both climate change and human driving forces, we predict the future land use map (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f1"="">1</a>(C)) using the RF algorithm. The midpoint year (2050) is chosen as a representative year for the future period during 2021–2080. The land use record (1992, 2015) and projection (2050) are used for subsequent comparative land use change analysis. According to the prediction, under more favorable climate to vegetation and agriculture as the primary driving force, along with socioeconomic development as the secondary driving force, the cultivated land area is to be expanded by 13 052 km<sup="">2</sup>, and the forest area is to increase by 237 km<sup="">2</sup> in 2050; meanwhile, grassland area is to decrease by 15 538 km<sup="">2</sup>, compared to the land use in 1992 (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f1"="">1</a>(A)). Overall, the human influences on the land will be rising through the transformation from herding to cultivation in the upper LRB. The interactions between climate, vegetation, and land use can cause unintended environmental consequences such as increased sediment discharge, as to be described in the following.

The spatiotemporal variation of sediment yield is assessed <em="">via</em> three SWAT modeling scenarios: one baseline and two benchmark scenarios (BMS) of alternative climate and land use conditions. The baseline scenario (BAS) assumes future climate (2021–2080) projected as a recurring cyclic pattern and uses the observed land use map of 2015 to represent the future land use, i.e. no change with neither climate nor land use from now to the future. BAS is used as a basis for comparison even though the assumptions are not likely realistic. BMS1 assumes future climate (following the GCMs projections), as shown in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#erlad1347s3-1"="">3.1</a>, and still uses the historical land use map in 2015. This scenario represents a situation that despite future climate change, there will be no change in land use due to mandatory land use control policies. BMS2 assumes future climate change and and land use change, i.e. without land use controls as enforced for BMS1. The land use modeling was based on the historical dataset, with land use as the result of both climate change and human intervention. Therefore, when we simulate future land use, the ongoing human intervention will persist into the future. In this way, the land use change will be the result of both future climate change and human intervention. In the SWAT model. The land use projections are periodically revised at 20 yr intervals to incorporate the changing landscape dynamics.

Figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f2"="">2</a> displays the 5-yr moving average of the simulated annual sediment yield under the three scenarios. The detailed 10-yr moving average of sediment yield under three scenarios across various climate change projections is presented in figure S10. The sediment yield increase under BMS1 compared to the BAS is estimated as 29.8% reflecting the individual impact of future climate change; and 40.6% under BMS2, indicating the synergistic impacts of climate change and associated land use continuing the historical changing trend. An estimated 6.9% increase in sediment yield is observed when comparing BMS2 to BMS1 under the same future climate conditions but differing in land use conditions (i.e. without or with land use control). Thus, the land use change is anticipated to result in a notable increase in sediment yield, corresponding to the projected climate change and land use change in the area. It is interesting to notice that BMS2 results in larger sediment yield increase than the sum of the increase between BMS1 and BAS and that between BMS2 and BMS1, i.e. 40.6% &gt; 36.7% (=29.8% + 6.9%). This shows the synergetic effect of climate change and land use change, as observed in the past and possibly continuing to the future. That is to say, climate change and land use change associated with climate change (if not controlled) will have the most significant impact on sediment yield increase in the area (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f3"="">3</a>(C)), with major contribution from cropland and barren land. The conversion from grassland and forest to cropland usually destroys dense root networks, leading to sparser vegetation cover and land surface vulnerable to erosion. In addition, more intense precipitation (i.e. floods) caused by climate change will result in a great magnitude of soil erosion. The M-K test confirms that the sediment discharge is likely to increase at the highest rate under BMS2 involving the synergistic climate and land use changes (table S6). Such a synergistic effect indicates that it may not be appropriate to separately assess the impact of climate change and land use, the two independent environmental variables, on the hydrological processes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides the increasing sediment yield trend of the entire basin as shown in figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f2"="">2</a> and <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlad1347f3"="">3</a> portrays the spatial distribution of sediment yield change under the BAS and BMS1 and BMS2. An overall increasing trend in sediment yield is found for most sub-areas with increased cropland and decreased grassland and forestland. Exceptions are observed for a small number of sub-areas where there is increased vegetation cover intensity due to warmer and wetter climate. For example, forest expansion is found in the central sub-areas, with relatively low elevation and slope, which benefits soil erosion control under changing environmental conditions. However, since the increase of cropland is far more than the increase of forest, the net sediment yield from the area will increase significantly under the project climate change and associated land use change without control.

The heterogeneity of sediment yield across the sub-areas highlights the necessity for distributed land management and adaptation planning. Areas exposed to the reductions of forest or grassland may require prioritization for land acquisition and protection efforts. The mid-south sub-basins and eastern sub-basins differ significantly in their hydrological responses to climatic and land use scenarios; the mid-south sub-basins require adaptation strategies to prevent high rates of soil losses.