Significant regime shifts in historical water yield in the Upper Brahmaputra River basin

Although evidence of the hydrological response of watersheds to climate change is abundant, reliable assessments of water yield (WY) over mountainous regions, such as the Upper Brahmaputra River (UBR) basin, remain unclear. Here, we examine long-term WY changes during 1982–2013 in the UBR basin, based on multi-station runoff observations. We find that there are significant shifts in hydrological regimes in the late 1990s; WY increases in the range of ∼10 % to ∼80 %, while the directions reverse from increasing to decreasing. Additionally, the double mass curve (DMC) technique is used to assess the effects of climate, vegetation, and cryosphere on WY changes. Results show that cryosphere and climate together contribute to over 80 % of the increase in WY across the entire UBR basin, while the role of vegetation is negligible. The combined effects, however, are either offsetting or additive, thus leading to slight or substantial magnitude increases, respectively. The downward WY trend has primarily been regulated by decreased precipitation in recent years. However, we find that meltwater may alleviate the resulting water shortage in some basins. Therefore, the combined effects of climate and cryosphere on WY should be considered in future water resources management over mountainous basins, particularly involving co-benefits between upstream and downstream regions.

1 Introduction

Water yield (or runoff depth; WY) in mountainous watersheds is crucial for sustaining fragile ecosystems in headwaters, supplying valuable freshwater resources to downstream lowlands, and balancing co-benefits between upstream and downstream areas, especially in large transboundary river systems (Viviroli et al., 2011). In mountainous watersheds, WY changes have been attributed to climate variability (Dierauer et al., 2018; Song et al., 2021), vegetation greening or browning (Goulden and Bales, 2014; Zhou et al., 2021), and glaciers and snow melting (Huss and Hock, 2018; Biemans et al., 2019). These changes are expected to alter the spatial and temporal distribution of water resources (Tang et al., 2019) and further threaten the water supply and food security downstream (Biemans et al., 2019). The Qinghai–Tibet Plateau (QTP; see Fig, 1a), also known as the Asian water tower, supplies water to major rivers in Asia, such as the Brahmaputra, Salween, Mekong, Yangtze, Yellow, and Indus rivers (Kang et al., 2010; Yao et al., 2010, 2019). Therefore, WY changes over this region significantly affect water availability and terrestrial and aquatic ecosystems, which are vital for sustaining the livelihoods of approximately 2 billion people (Immerzeel et al., 2010). Despite some in situ observations and estimates from state-of-the-art remote sensing (Wang et al., 2021), total river runoff has never been reliably assessed in this region, and its response to global warming remains unclear. Therefore, comprehensively assessing long-term changes in WY, particularly magnitude and direction, is of great importance for the sustainable development of water resources in the QTP (Yao et al., 2019).

In recent years, WY has been significantly affected by multiple factors in the QTP. For example, Fan and He (2015) highlighted the important role of precipitation in WY increases in the Salween River basin and Mekong River basin. Li et al. (2020) found that elevated precipitation and meltwater both contributed to substantial WY increases in the Tuotuo River basin. Similarly, Lutz et al. (2014) projected that the increased precipitation near the Salween and Mekong rivers and the accelerated meltwater near the Indus River caused significant WY changes. Vegetation changes have also proven to be vital for mountainous water resources. Li et al. (2017) showed that evaporation, mostly due to grassland restoration, decreased WY in the Yangtze River basin. Furthermore, Li et al. (2021) suggested that vegetation greening may change the seasonality of water resources and increase WY during the dry season in the Upper Brahmaputra River (UBR) basin.

Although a growing body of evidence has shown that WY is significantly changing in the QTP, most previous studies focus on individual basins, which may not fully reveal the spatial variability in WY. Of specific interest is the Upper Brahmaputra River (UBR; see Fig. 1b) basin, which covers an area of over 198 636 km² and has large gradients in elevation, which has influenced climate and vegetation patterns (Li et al., 2019b; Gao et al., 2018). Hence, it is imperative to provide a comprehensive, spatially differentiated study on long-term WY changes in this region. Such studies in the UBR basin, however, are significantly hindered by the sparse network of hydrological observation stations (Li et al., 2019b; Wang et al., 2021; Yao et al., 2019), which leads to large uncertainties in river flow predictions and, thus, water resources assessments. Also, current precipitation estimates are highly uncertain, owing to the complex topography of this region, which limits the ability to accurately model the relationships between precipitation and runoff (Sun and Su, 2020). Last, the current inadequate understanding of hydrological responses to complex interactions among multi-spheres limits the application of hydrological models in these mountainous watersheds (Pellicciotti et al., 2012). In this regard, long-term observed runoff records and recent high-resolution precipitation data in the UBR basin provide a valuable opportunity to estimate runoff responses to warming using statistical methods.

This study uses the double mass curve (DMC) technique to jointly assess historical WY responses to climate warming and the associated environmental changes in the UBR basin. To do this, we collect multi-station runoff observations and detect long-term WY changes during 1982–2013. Then, we use the DMC method to estimate the effects of climate (represented by effective precipitation, P−E, or eP), vegetation (represented by leaf area index – LAI), and cryosphere (e.g., meltwater from glaciers and snow) on the magnitude and direction changes in WY (see Sect. 2). This study is expected to provide essential information for water resources management in the UBR basin and other mountainous watersheds.

2 Data and methods

2.1 Study area

The Brahmaputra River (known as the Yarlung Zangbo River, or YZR, in China), a transboundary river in the southern QTP, originates in the Gyama Langdzom Glacier and flows across China, India, and Bangladesh before emptying into the Indian Ocean. The UBR basin is located above the Nuxia hydrological station (Fig. 1), and its river flow has significant implications for the freshwater resources of South Asia. Here, we divide the UBR basin into the headstream (HYZR), upstream (UYZR), midstream (MYZR), downstream (LYZR), Nianchu River (NCR) basin, and Lhasa River (LSR) basin by the locations of hydrological stations (Table 1 and Fig. 1b).

Table 1 Information of six basins divided according to the locations of hydrological stations. The column Tp indicates the turning point, calculated using the Pettitt method, in which a significant turning point is labeled with an asterisk^*. Glaciers and snow coverage are derived from the land use and cover data in 2000 (see Sect. 2.2). The unit of area is kilometers squared (km²), and the unit of elevation is meters (m).

Figure 1 Location of (a) the Upper Brahmaputra River (UBR) basin in the Qinghai–Tibet Plateau (QTP), which is from Li et al. (2021), and (b) six basins divided by the Lhatse, Nugesha, Shigatse, Yangcun, Lhasa, and Nuxia hydrological stations.

The elevation gradient and the distance to the ocean in the UBR basin together contribute to a large spatial variability in climate (Sang et al., 2016). The annual precipitation in the HYZR basin is less than 400 mm, while that in the LYZR basin is nearly 1000 mm. Similarly, the annual evaporation increases gradually from upstream to downstream areas (Fig. S1a and b). Meanwhile, water and energy availability modulate vegetation conditions (Li et al., 2019a). Vegetation cover increases dramatically from the HYZR to the LYZR basin (Fig. S1c). Additionally, cryospheric meltwater due to atmospheric warming has substantially altered hydrological conditions in this region (Cuo et al., 2019; Yao et al., 2010; Wang et al., 2021).

2.2 Data

Here, we collect multi-station annual runoff observations from 1982 to 2013 and convert the river flow (m³ s⁻¹) into runoff depth (mm). In addition, we acquire high-resolution climate and vegetation data over the same time period and further aggregate these gridded datasets into annual values regionally by considering area-weighted effects (their temporal changes are shown in Fig. S2).

2.2.1 Runoff data

Annual runoff observations between 1982 and 2013 used here come from six hydrological stations along the mainstream and major branches in the UBR basin. WY in the HYZR is equal to the runoff depth at Lhatse station, while WY in other basins is calculated by the difference between runoff depth from the downstream station and that from the upstream and branch stations. For example, WY in the MYZR basin is equal to the difference between runoff depth in Yangcun station and that in Lhasa and Nugesha stations (Fig. 1b).

2.2.2 Climate data

The most recent 10 km gridded daily precipitation data, combining topographic and linear correction approaches based on 262 rain gauge observations, is developed for the UBR basin by Sun and Su (2020) and here used to estimate regional annual precipitation (P). The maximum 2 m air temperature is obtained from China Meteorological Forcing Dataset (He et al., 2020). The evaporation (E) with a 0.25^∘ spatial resolution is acquired from Global Land Evaporation Amsterdam Model (GLEAM) version 3.5a (Martens et al., 2017). GLEAM evaporation has been validated in different biome types in China and has shown high correlations with in situ eddy covariance E (Yang et al., 2017).

2.2.3 Vegetation data

Leaf area index (LAI) used in this study is obtained from Global Inventory Modeling and Mapping Studies (GIMMS) LAI3g, with a spatial resolution of 8 km (Zhu et al., 2013). The product is generated using an artificial neural network trained on the collection Terra Moderate Resolution Imaging Spectroradiometer (MODIS) LAI product and the latest version of GIMMS NDVI3g data for the same period, which has been proven to have an improved multi-sensor record harmonization scheme compared to other global LAI products (Forzieri et al., 2020; Gonsamo et al., 2021; Zhu et al., 2016).

2.2.4 Land use and cover

The land use and cover in 2000, with a spatial resolution of 1 km, is used to represent the land cover types in the UBR basin. The data are acquired from the Resource and Environmental Science Data Center and are here divided into seven primary land use types, including cultivated land, forestland, grassland, waterbody, urban land, unused land, and glaciers and snow (Fig. S3).

2.3 Methodology

2.3.1 Trend and turning point analysis

In this study, we use the non-parametric Mann–Kendall test (Kendall, 1938; Mann, 1945), to identify the trend in WY, and Pettitt change detection (Pettitt, 1979), to identify the turning point (Tp) in WY, where the level of significance is set at 0.05. The Pettitt method has the advantages of being non-parametric, rank-based, and distribution-free for finding the turning point in a time series. We then compare the averages of WY before and after Tp to reflect the changes in magnitude and compare the trends in the two periods to represent the direction changes.

2.3.2 Double mass curve technique

The DMC used here is a plot of the cumulative data of one variable versus the cumulative data of another related variable in a concurrent period. It has previously been used to assess the individual effect of climate (Gao et al., 2011), forest disturbance (Wei and Zhang, 2010), wildfire (Hallema et al., 2018), and cryosphere (Brahney et al., 2017) on water resources. For the large and pristine UBR and other mountainous basins, climate, vegetation, and cryosphere (meltwater from glaciers and snow under warming; see Biemans et al., 2019; Huss and Hock, 2018) play important roles in the water cycle, and these three components must be considered together to accurately estimate the hydrological response to warming. It is difficult to directly calculate the supply of meltwater to WY due to the lack of long-term glacier monitoring. On the other hand, runoff observations and high-resolution climate and vegetation data make it possible to use the DMC technique, a data-driven statistical method, to estimate cryospheric contributions to WY.

Appropriate selection of climate and vegetation indices used in the DMC technique is of great importance. Previous studies have shown that effective precipitation (P−E; eP) can adequately represent the impact of climate on WY compared to individual estimates of P or E (Wei and Zhang, 2010; Zhang et al., 2019). LAI quantifies the amount of leaf area in an ecosystem and is an important variable reflecting vegetation structures and biophysical processes (Forzieri et al., 2020), and Li et al. (2021) have used LAI to investigate vegetation effects on seasonal hydrology in the UBR basin. Hence, we consider eP and LAI as being the indices of climate and vegetation, respectively, and use their time series as the inputs in the DMC model.

To obtain cryospheric contributions to WY, we firstly build two types of DMC plots (see Fig. S4) to assess the contributions of climate and vegetation and then subtract the sum of estimated contributions from total WY deviations as cryospheric effects (results are shown in Fig. S5). The mathematical formulas, relevant descriptions, and an example (Fig. 2) are as follows.

Figure 2 The example showing how to estimate the effects of climate, vegetation, and cryosphere on water yield in the MYZR basin (see details in Sect. 2.3), where Tp is 1997.

First, the average WY before Tp ($\stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}}$; horizontal black line in Fig. 2b) is defined as follows:

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}}={\displaystyle \frac{\sum _{t=\mathrm{1982}}^{t=\mathrm{Tp}}{\mathrm{WY}}_{\mathrm{o}}\left(t\right)}{\mathrm{Tp}-\mathrm{1982}+\mathrm{1}}}.\end{array}$$

Next, the total WY deviation (ΔWY_s(t); black diamond in Fig. 2c) can be calculated as the difference between WY observations after Tp (WY_o(t); red point in Fig. 2b) and the average WY before Tp ($\stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}}$), as follows:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{\mathrm{WY}}_{\mathrm{s}}\left(t\right)={\mathrm{WY}}_{\mathrm{o}}\left(t\right)-\stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}},t=\mathrm{Tp}+\mathrm{1},\mathrm{Tp}+\mathrm{2},\mathrm{\dots },\mathrm{2013}.\end{array}$$

Second, the regression equation (left panel in Fig. 2a) between the cumulative eP (∑eP) and cumulative WY (∑WY) before Tp can be constructed as follows:

$$\begin{array}{}\text{(3)}& \sum \mathrm{WY}=a_{\mathrm{1}}\sum \mathrm{eP}+b_{\mathrm{1}}.\end{array}$$

Similarly, the regression equation (right panel in Fig. 2a) between cumulative LAI (∑LAI) and cumulative water yield (∑WY) before Tp can be constructed as follows:

$$\begin{array}{}\text{(4)}& \sum \mathrm{WY}=a_{\mathrm{2}}\sum \mathrm{LAI}+b_{\mathrm{2}}.\end{array}$$

Third, WY driven by climate (WY_c(t); blue line in Fig. 2b) can be calculated by using the cumulative eP after Tp as input into Eq. (3), which is built using the cumulative data of WY and eP before Tp. Therefore, WY deviations caused by climate change (ΔWY_c(t); blue bar in Fig. 2c) can be calculated as follows:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}\left(t\right)={\mathrm{WY}}_{\mathrm{c}}\left(t\right)-\stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}},t=\mathrm{Tp}+\mathrm{1},\mathrm{Tp}+\mathrm{2},\mathrm{\dots },\mathrm{2013}.\end{array}$$

Similarly, WY driven by vegetation (WY_v(t); tan line in Fig. 2b) can be calculated via Eq. (4), and WY deviation caused by vegetation (ΔWY_v; tan bar in Fig. 2c) can be calculated as follows:

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}\left(t\right)={\mathrm{WY}}_{\mathrm{v}}\left(t\right)-\stackrel{\mathrm{‾}}{{\mathrm{WY}}_{\mathrm{ob}}},t=\mathrm{Tp}+\mathrm{1},\mathrm{Tp}+\mathrm{2},\mathrm{\dots },\mathrm{2013}.\end{array}$$

Finally, WY deviation caused by cryosphere (ΔWY_g; red bar in Fig. 2c) can be calculated as follows:

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}\left(t\right)=\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{s}}\left(t\right)-\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}\left(t\right)-\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}\left(t\right).\end{array}$$

2.3.3 Attribution analysis on changes in water yield

The average contributions of climate, vegetation, and cryosphere to WY magnitude changes are calculated as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}}={\displaystyle \frac{\sum _{t=\mathrm{Tp}+\mathrm{1}}^{t=\mathrm{2013}}\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}\left(t\right)}{\mathrm{2013}-\mathrm{Tp}}}\\ & \stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}}={\displaystyle \frac{\sum _{t=\mathrm{Tp}+\mathrm{1}}^{t=\mathrm{2013}}\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}\left(t\right)}{\mathrm{2013}-\mathrm{Tp}}}\\ & \stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}}={\displaystyle \frac{\sum _{t=\mathrm{Tp}+\mathrm{1}}^{t=\mathrm{2013}}\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}\left(t\right)}{\mathrm{2013}-\mathrm{Tp}}}.\end{array}\end{array}$$

The relative contribution (RC), ranging from 0 to 1, of climate, vegetation, and cryosphere can be calculated as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\mathrm{RC}}_{\mathrm{c}}={\displaystyle \frac{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}}\right|}{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}}\right|}}\\ & {\mathrm{RC}}_{\mathrm{v}}={\displaystyle \frac{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}}\right|}{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}}\right|}}\\ & {\mathrm{RC}}_{\mathrm{g}}={\displaystyle \frac{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}}\right|}{\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{c}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{v}}}\right|+\left|\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{WY}}_{\mathrm{g}}}\right|}}.\end{array}\end{array}$$

Additionally, we use the Pearson correlation coefficient (r) to quantify the relationships between total WY deviation (ΔWY_s(t)) and its components, where the WY deviation is caused by climate (ΔWY_c(t)), vegetation (ΔWY_v(t)), and cryosphere (ΔWY_g(t)). The Student's t test with the Bonferroni correction (Bonferroni, 1935) is used to detect the statistical significance of the correlation coefficient at the level of 0.05.

3 Results

3.1 Long-term changes in historical water yield

The detection of long-term WY changes in the UBR basin from 1982 to 2013 is illustrated in Fig. 3. There is a high degree of spatial variability in WY (Fig. 3a) across the entire basin. The average WY is the highest in the LYZR basin (over 600 mm), followed by that in the LSR basin (nearly 400 mm). However, the average WY in the HYZR and NCR basins is less than 100 mm. This spatial variability is consistent with that of precipitation (Fig. S6), which is mainly determined by elevation and distance to the ocean (Sang et al., 2016). In addition, WY generally increases during the period, as shown by the positive slope in Fig. 3b, even though the significant trend is only detected in the UYZR and MYZR basins (hatched areas in Fig. 3b). Tp mainly occurs during the late 1990s (Fig. 3c), but the abrupt change detected in some basins is not statistically significant (Table 1). Similarly, the cumulative anomaly curve (Fig. 3d) shows that WY decreases prior to the late 1990s and then increases in the entire UBR basin, which further supports the results obtained from the Pettitt method.

Figure 3 Long-term water yield changes in the entire UBR basin. (a) The annual values calculated by averaging water yield from 1982 to 2013. (b) The temporal trends detected by the Mann–Kendall Sen's slope method. The black hatching represents statistically significant (p<0.05) trends. (c) The turning points detected by the Pettitt method. (d) The cumulative anomaly curve of water yield. The solid green line represents the ensemble expectation of cumulative anomaly curves of water yield in the entire UBR basin (green shading).

3.2 Regime shifts in historical water yield

Based on the results of the abrupt detection, we divide the period from 1982 to 2013 into the before and after Tp period and analyze the magnitude and direction of WY changes in the entire UBR basin. Figure 4a shows that WY increases from 9.5 to 130.9 mm, with a high spatial variability. The slight increase observed in the HYZR and LYZR basins accounts for less than 10 % of the mean WY before Tp, while a substantial WY increase of 61.6 % and 80.5 % is found in the UYZR and MYZR basins, respectively. In addition, higher variability is detected after Tp, suggesting a more dramatic variation in recent years. With respect to the direction of change in WY, we find that the trend is positive before Tp but becomes negative afterward in most basins (Figs. 4b and S2). The significantly decreasing trend is detected in the UYZR, NCR, and LSR basins. In contrast, although the WY in the MYZR basin increases during two periods, the rate has slowed; the positive trend after Tp (3.64 mm yr⁻¹; p>0.05) is less than that before Tp (8.95 mm yr⁻¹; p<0.05). Overall, WY regime shifts occur in the late 1990s in the entire UBR basin because the magnitude generally increases, but the direction of WY has reversed or slowed.

Figure 4 Water yield regime shifts in the entire UBR basin. (a) Magnitude of water yield changes. The black cross (x) shows the mean of the water yield, and the relative change is labeled in number (%) in each box plot. (b) Direction of water yield changes. The black hatching represents the statistically significant trend (p<0.05). The color of the boxes represents the period before (light color) and after (dark color) the Tp period.

3.3 Attribution analysis of magnitude increases in water yield

In Fig. 5, we quantify the contributions from climate (eP), vegetation (LAI), and cryosphere on WY magnitude increases in the entire UBR basin. We find that cryosphere changes, on average, contribute to over half of the WY increases in the HYZR, UYZR, NCR, and MYZR basins. However, climate plays a more important role in WY magnitude increases in the LSR and LYZR basins, with average relative contributions of 55.4 % and 46.0 %, respectively. In contrast to the dominant roles of climate and cryosphere, vegetation has a consistently positive contribution to the WY increases, although its relative contribution is much less than those from climate and cryosphere.

Figure 5 Attribution analysis of magnitude increases in water yield due to climate (ΔWY_c; blue box), vegetation (ΔWY_v; tan box), and cryosphere (ΔWY_g; red box) and their relative contributions (the bar with colors on the top) in each basin. The black cross (x) shows the average of water yield deviations (see Fig. S5) in each box plot.

Climate and cryosphere – two important factors affecting WY – together contribute to over 80 % of the average increases in WY magnitudes over the entire UBR basin. However, they play either additive or offsetting roles in different basins (Fig. 5) and thus result in slight or substantial WY increases (Fig. 4a). For example, although cryospheric changes contribute to an increase in WY in the HYZR (28.3 mm) and NCR basins (30.3 mm; black cross (x) in Fig. 5a and c), the negative contributions from climate offset a considerable part of these increases, leading to only slight increase in WY after Tp in these regions. Similarly, the positive contribution from climate offsets the negative contribution from cryosphere in the LSR and LYZR basins, which also results in a slight WY increase. In addition, the additive effects from climate and cryosphere lead to substantial increases in WY, from 162.6 to 293.5 mm in the MYZR basin and from 164.9 to 266.5 mm in the UYZR basin (Fig. 4a).

3.4 Attribution analysis of direction shifts in water yield

The Pearson correlation coefficient is applied to determine the role of climate, vegetation, and cryosphere in the reversal or slowing of WY trends after Tp, as shown in Fig. 4b. Results in Fig. 6 show that, although the correlation varies greatly across basins ranging from 0.11 to 0.93 after Tp, the climate typically is positively associated with total WY, for which the correlation is significant in half of the basins (p<0.05), again revealing the major role of climate in influencing the hydrological trends across the entire UBR basin. Further analysis shows that precipitation is much more important because it exhibits a stronger reversal in trend compared with the trend in evaporation (Fig. S7), which is also similar to the directional changes in WY (Fig. 4b). Additionally, despite the weak contribution of vegetation (Fig. 5), its positive role in WY changes is more apparent in the drier basins (such as the HYZR, UYZR, and NCR basins), while the correlation is negative in the relatively humid LYZR basin.

Figure 6 The correlation between the time series of total water yield deviation (ΔWY_s(t); x axis) and its components (y axis) caused by climate (ΔWY_c(t); blue point), vegetation (ΔWY_v(t); tan point), and cryosphere (ΔWY_g(t); red point), respectively. The fitting line and its 95 % confidence interval are shown only when the p value < 0.05. n indicates the number of years after Tp, which is determined by the Pettitt method (see Table 1 and Fig. 3c).

In contrast to the positive contributions of climate, we find that WY caused by cryosphere exhibits a negative association with reduced total WY deviations in recent years in the UYZR ($r=-\mathrm{0.39}$; p>0.05) and LSR ($r=-\mathrm{0.36}$; p>0.05) basins. The negative but weak relationship indicates that meltwater from cryospheric loss may compensate for low flow and even mitigate water shortage risks. Also, the compensating effect from cryosphere is stronger in the MYZR (r=0.47; p>0.05) and, together with climate contributions, contributes to the increasing WY trend (Fig. 4). In contrast to other regions, however, the HYZR basin shows a significantly positive relationship between cryospheric contributions and total WY deviations (r=0.76; p<0.05), indicating that the cryosphere instead of climate leads to the downward trend in headwaters. This signifies that in this region, cryospheric contributions may have already passed its potential for supplying to river flow, due to decreased glaciers and snow under continuous warming. This is further verified by the relationship of cryospheric contributions to total WY (RC_g) with temperature (Fig. S8a). The cryospheric contribution to WY increases with temperature in the beginning, until a maximum is reached, and then begins to decrease. In addition, the compensating effect of meltwater can be seen in the UYZR, MYZR, and LSR basins, i.e., WY caused by cryospheric changes keeps a positive relationship with the increase in temperature, further supporting the observed higher correlation in these basins (Fig. 6).

4 Discussion

4.1 Climate and cryosphere cause water yield regime shifts

Previous studies have reported the increasing WY trend in the LSR (Lin et al., 2020), LYZR (Zhang et al., 2011), and UBR basins (Li et al., 2021). Here, we not only provide further evidence on the long-term trend of WY changes in the above regions but also conduct trend analyses in other regions that have received less attention. As such, the study comprehensively indicates a general increase in WY (Fig. 3a) in the entire UBR basin. Furthermore, we find that regime shifts in WY occur during the late 1990s in the entire UBR basin, where the magnitude in WY increases (Fig. 4a), but the direction has reversed or slowed after Tp (Fig. 4b). The result agrees with climate shifts (Fig. S9), drought changes (Li et al., 2019b), and lake area changes (Zhang et al., 2017) in this region.

Attribution analysis of WY magnitude shifts highlights the dependency of spatial gradients of climate and cryosphere accounting for WY changes. Climate explains a greater increase in WY in downstream regions, while meltwater is more important in upstream regions (Fig. 5). This may be attributed to divergent water supply sources. Biemans et al. (2019) indicate that meltwater supplies over 40 % of river flow in the upper regions, but the contribution is less than 30 % in the downstream regions of the UBR basin (Fig. S9). Vegetation does not significantly affect the total river flow in the UBR (Fig. 5), despite its remarkable role in seasonal runoff detected by Li et al. (2021).

Climate, especially precipitation, still controls the declining WY trend after Tp in most regions (Figs. 6 and S7) and may become an important factor in occurrence of turning points (Fig. 3c and d). This suggests the importance of precipitation and its projections on future hydrological process in mountainous watersheds (Lutz et al., 2014). The cryospheric contribution is also important for water yield regime shifts, as meltwater from glaciers and snow can alleviate water resources shortages, which have mainly been caused by decreased precipitation in recent years (Figs. 6 and S7). This finding is also supported by observed glacier runoff data (Yao et al., 2010) and several modeling studies (Lutz et al., 2014; Zhang et al., 2020; Wang et al., 2021). However, after glacier runoff reaches a maximum (defined as peak water; Gleick and Palaniappan, 2010), cryospheric mass loss cannot sustain the rising meltwater with atmospheric warming (e.g., the HYZR basin in Fig. S8), which is in agreement with Huss and Hock (2018).

4.2 Uncertainties and limitations

This study has some limitations pertaining to the DMC model used for partitioning the effects of climatic and cryospheric changes on the hydrological regime shifts in the UBR basin. The DMC method is a useful statistical alternative to physical modeling approaches, especially in alpine river basins (e.g., UBR basin) where there is inadequate knowledge of the complex hydrological mechanisms. However, the method is still dependent on our prior understanding of hydrological responses to warming and related environmental changes, such as glaciers melting and vegetation changes. For the UBR basin, besides climate change, cryosphere (Biemans et al., 2019; Yao et al., 2019) and vegetation (Li et al., 2021, 2019a) are two major factors which drive hydrological changes. Here, the cryospheric contribution is estimated as being the deviation between total water yield and the climate and vegetation contribution estimated by the DMC method. However, in some mountainous basins, human activities such as urbanization, dam regulation, and irrigation may consume a significant portion of water resources or induce changes in seasonal runoff patterns. Therefore, it is necessary to account for anthropogenic impacts when assessing river flow changes via the statistical models in these regions, in addition to the factors considered in this study. In addition, the DMC method assumes a linear relationship between two variables, and thus it may fail to capture nonlinear interactions among water, vegetation, and cryospheric melting. For example, the role of vegetation in influencing hydrological responses is expected to be complex due to biophysical (e.g., via transpiration and albedo changes) and biochemical (e.g., via CO₂ uptake and release) feedbacks (Bonan, 2008; Krich et al., 2022; Li et al., 2021). Thus, with the availability of long-term in situ observations and high-resolution remote sensing datasets in the UBR basin (Wang et al., 2022), other powerful statistical models which account for nonlinear and casual relationships should be applied to identify the causes of water yield changes (Runge et al., 2019).

The data used in the study are another source of uncertainty in our results. The precipitation dataset used here is generated by topographical and linear corrections based on observations. As Sun and Su (2020) have pointed out, the results of the linear correction approach vary significantly with station density. For example, increasing the number of stations from 4 to 10 in the basin results in a decrease in the mean annual precipitation estimate by about 20 mm. Hence, the reconstruction of the precipitation dataset will rely on the density of the observed stations. Besides the topographic correction, the effects of basin size and climate seasonality should be considered in the work of precipitation reconstruction in the UBR basin due to its complex climate and environment (Sun et al., 2019). Compared with precipitation, the estimation of evaporation may be much more challenging in high mountains. Although GLEAM evaporation shows good agreement with in situ eddy covariance records (Yang et al., 2017), its model structure does not include wind speed and solar radiation, which may affect the estimation of sublimation and, thus, total evaporation (X. Li et al., 2019). In addition, the coarse spatial resolution (0.25^∘) of GLEAM may be insufficient to estimate regional evaporation in the UBR basin. GIMMS LAI has the advantage of capturing the ecosystem structure and thus is widely used to assess vegetation conditions and their effects on hydrological changes (Zhu et al., 2016; Forzieri et al., 2020; Gonsamo et al., 2021). However, LAI ignores the vegetation's physiological process (Fang et al., 2019). Hu et al. (2022) have indicated that LAI can cope with hydroclimatic fluctuations in arid environments, but the tradeoff between ecosystem structure (LAI) and physiology (photosynthesis per unit leaf area) is stronger in humid climates. Thus, using LAI products in energy-limited regions may result in some biased assessments of the vegetation effects on water yield. In spite of the uncertainties in the methods and data used here, this study provides essential information for water resource management in the UBR basin based on multi-station runoff observations. With the help of the Second Tibetan Plateau Scientific Expedition and Research, the observation networks in meteorology, cryosphere, and hydrology will be expanded, which can potentially lead to a better understanding of the hydrological processes in this region and make developing physically based cryosphere–hydrological modeling possible (Wang et al., 2022).

4.3 Broad implications for mountainous water management

Understanding the hydrological regime shifts and their causes in high mountain regions is especially important in managing water resources. Of specific importance is balancing the co-benefits between mountains and downstream lowlands (Viviroli et al., 2011). In this study, the combined (both offsetting and additive) effects of climate and cryosphere are detected (Fig. 5). They further lead to either slight or substantial increases in WY across the entire UBR basin (Fig. 4a). The combined effects often hinder the role of each driver in influencing hydrological changes (Wei et al., 2018; Zhang and Wei, 2021) and thus should be considered when designing water management strategies in large transboundary river systems. For example, the additive effect may be beneficial for mitigating droughts and water shortage during droughts, but it may exacerbate flood risks due to increased precipitation and accelerated melting of the cryosphere in the future (Immerzeel et al., 2013). In addition, in the headwaters, WY resulting from cryosphere melting continues to increase with temperature until a maximum is reached, beyond which cryospheric contribution to total WY begins to decrease (Fig. S8a). This shows that the meltwater from glaciers may have already surpassed the peak water time. Such hydrological changes will substantially affect future water resources management, and thus accurate projections of the occurrence time of peak water will be important in managing mountainous water resources.

5 Conclusions

In this study, regime shifts in historical WY are detected during the late 1990s in the UBR basin. The WY magnitude generally increases, but its direction has reversed or slowed in recent years. Attribution analysis based on the DMC method shows that the combined effects of climate and cryosphere are either offsetting or additive, leading to slight or substantial WY increases, while the role of vegetation is much weaker. The declining or slowing WY trend after Tp is mainly driven by precipitation in most regions, while meltwater may alleviate drought and water shortage. In headwaters, however, cryospheric contributions to WY have declined due to reduced glaciers under warming. These findings suggest that the combined effects of climate and cryosphere should be considered in the sustainable development of water resources, especially involving co-benefits in upstream and downstream regions.