Changes in nonlinearity and stability of streamflow recession characteristics under climate warming in a large glaciated basin of the Tibetan Plateau

The accelerated climate warming in the Tibetan Plateau after 1997 has profound consequences in hydrology, geography, and social wellbeing. In hydrology, the change in streamflow as a result of changes in dynamic water storage that originated from glacier melt and permafrost thawing in the warming climate directly affects the available water resources for societies of the most populated nations in the world. In this study, annual streamflow recession characteristics are analyzed using daily climate and hydrological data during 1980–2015 in the Yarlung Zangbo River basin (YRB) of the southern Tibetan Plateau. The recession characteristics are examined in terms of $\mathrm{d}Q/\mathrm{d}t$ = − aQ^b and the response/sensitivity of streamflow to changes in groundwater storage. Major results show that climate warming has significantly increased the nonlinearity of the response (b) and streamflow stability [log(a)] in most subbasins of the YRB. These changes in the recession characteristics are attributed to the opposite effects of increases in the available water storage and recession timescale on the recession. Climate warming has increased subbasin water storage considerably due to more recharge from accelerated glacier melting and permafrost thawing after 1997. Meanwhile, the enlarged storage lengthens recession timescales and thereby decreases the sensitivity of discharge to storage. In the recession period when recharge diminished, increased evaporation and the decreased buffering effect of frost soils under warmer temperatures accelerate the initial recession of streamflow. By contrast, enlarged storage and lengthened recession timescales slow down the recession. While reservoir regulations in some basins have helped reduce and even reverse some of these climate warming effects, this short-term remedy can only function before the solid water storage is exhausted should the climate warming continue.

1 Introduction

The warming rates of air temperature in high latitudes and high altitudes are greater than the rate of change in global average near-surface air temperature (e.g., Pepin et al., 2015). The greater climate warming has accelerated glacier melting and permafrost thawing in cold alpine regions, causing significant glacier and permafrost retreats (e.g., Yao et al., 2004, 2007) and landscape alternations (e.g., Niu et al., 2019). These changes undoubtedly alter the hydrodynamics of streamflow and groundwater storage in the alpine regions and their downstream tributaries (Bense et al., 2012; Walvoord and Striegl, 2007; Walvoord and Kurylyk, 2016; Li et al., 2018; Wang et al., 2018; Yi et al., 2021). A recent study of Wang et al. (2021) has shown that such changes have also caused changes in the precipitation–streamflow relationship. The societal impacts of these changes are profound because they redefine freshwater availability and its seasonality for populations of billions in the downstream tributaries (e.g., Cuo et al., 2014; Zhang et al., 2013; Wang et al., 2020).

Many studies have found that the compound effects of glacier and permafrost retreats in the past few decades have reshaped the groundwater flow and hydrological cycle (e.g., Bring et al., 2016; Forster et al., 2014; Ji et al., 2020; Walvoord and Kurylyk, 2016). The accelerated glacier melting and permafrost thawing have increased the soil active layer thickness (ALT) and therefore enlarged groundwater storage and allowed the exchange of surface water and groundwater (Xu et al., 2017; Forster et al., 2014; Ji et al., 2020). Such an exchange further alters streamflow composition in arctic catchments (Chang et al 2008; Walvoord and Kurylyk, 2016; Bring et al., 2016) and in the northeastern and southern Tibetan Plateau (TP; Li et al., 2018; Wang et al., 2018; Yi et al., 2021).

In those catchments, changes in groundwater storage and the subsurface moisture profile due to permafrost retreat could reroute subsurface flow paths during low flows (Koch et al., 2014; Payn et al., 2012). The thickened ALT allows infiltration through the previously permafrost layer into aquifers and thus increases the recharge in permafrost basins. It has been reported that the permafrost loss in the past decades has enhanced regional groundwater circulation (shortening its timescale) with more discharge to stream flows (Ji et al., 2020; Walvoord et al., 2012). As an example, the baseflow in the source region of the Yangtze River increased at a rate of 1.35 mm a⁻¹ during 1962–2012, following the annual temperature rise of 1.32 ^∘C (Yi et al., 2021) and 1.09 mm a⁻¹ during 1979–2013 in glacierized basins in the TP, with its annual temperature rise of 0.98 ^∘C (Lin et al., 2020). Walvoord and Striegl (2007) found that the groundwater contribution to streamflow in an arctic basin increased by 0.7 %–0.9 % per year from the 1950s to 2005, following its annual temperature rise of 1.24 ^∘C during that period.

The effects of these changes in water budget and subsurface moisture profile on hydrographs are complicated in high altitude/frozen areas. The increase in ALT could reduce the buffering effect of soils on streamflow variability and thereby increase the baseflow recession rate (Lyon et al., 2009; Lyon and Destouni, 2010; Brutsaert and Hiyama, 2012). On the other hand, the increase in ALT enlarges groundwater storage and subsequently strengthens aquifer regulations on groundwater flow and slows the recession rate (Lin et al., 2020; Mao and Wang, 2016). These effects on the recession rate can result in the strongly nonlinear behavior of streamflow in time and space.

During periods of little or no precipitation, the baseflow recession ($\mathrm{d}Q/\mathrm{d}t$ vs. Q, where Q is discharge) is typically described by a power law $\mathrm{d}Q/\mathrm{d}t$ = − aQ^b (Brutsaert and Nieber, 1977; Tallaksen, 1995). The depletion of baseflow in relation to the parameters a and b contains valuable information concerning storage properties and aquifer characteristics of basins (Tallaksen, 1995). The recession-scale parameter a is a function of the hydraulic and geometric properties of the aquifer of a basin and can be used as a proxy for determining the effective depth of permafrost in frozen areas (Lyon and Destouni, 2010). The parameter b as reflected in the concavity of the hydrograph or the nonlinearity of recession (Dralle et al., 2017) is a function of boundary conditions to describe the equivalent water depth profile of an aquifer (Brutsaert and Nieber, 1977; Tashie et al., 2020). So, b can be interpreted as a measure of the diversity of water transport timescales throughout various parts of a catchment (Harman et al., 2009). Therefore, the variations in a and b in time and space can describe recession characteristics of a basin (e.g., Brutsaert and Nieber, 1977; Kirchner, 2009).

For the effects of a and b on the recession characteristics, Tashie et al. (2019) defined log(a) as the stability of streamflow and b as recession nonlinearity. For individual events, an increase in the relative value of a between events indicates decreased streamflow stability, while an increase in b indicates increased nonlinearity. Changes in the recession characteristics in time reflect their vulnerability to climatic and anthropogenic forcings (Berghuijs et al., 2016; Brooks et al., 2015; Buttle, 2018). Streamflow stability, log(a), has a significant seasonal cycle in over 99 % of the basins in the continental United States (Tashie et al., 2020). Moistening climate in catchments could increase the diversity of flow paths and the nonlinear relationships between storage/recharge and discharge (Brutsaert and Nieber, 1977; Hinzman et al., 2020). In cold climate regions, reduced glacier size can lead to considerable amplification of the seasonality of streamflow (Juen et al., 2007; Vuille et al., 2008). Hinzman et al. (2020) reported a widespread increase in the nonlinearity of recessions in northern Sweden due to climate warming. In addition, they found that this nonlinearity is significantly higher in warm winters than in cold winters. These analyses were performed under the assumption that a and b are effectively decorrelated. However, because changes in a in time are dependent on discharge Q, log(a) may be strongly correlated with b in individual events (Dralle et al. 2015; Biswal, 2021). Dralle et al. (2015) proposed a rescaling technique that eliminates the scale dependence of fitted power law parameters. Biswal (2021) selected the median of the values of b in individual events as the fixed b value and used it to find the event values of a.

Figure 1 (a) Geographical location of the Yarlung Zangbo River basin (YRB; the entire basin is above the hydrological station of Nuxia) and its five subbasins (NGS for Nugesha, YC for Yangcun, NX for Nuxia, YBJ for Yangbajain, and LS for Lhasa) from upstream to downstream. (b) Distributions of permafrost and seasonally frozen ground in 2012 and active tensile faults (the red dotted lines).

In the southern TP, the Yarlung Zangbo River basin (YRB; Fig. 1a) has decades of observations and offers an opportunity to estimate variations in the recession characteristics of streamflow in glaciated areas of TP. Using those data, some recent studies have shown that the climate of the area has become warmer and wetter from 1980–2015 (e.g., Wang et al., 2021). Climate warming has reduced the buffering effect of glacial and permafrost on streamflow, leading to catchment property change with a shorter streamflow response time to precipitation (Wang et al., 2021). These changes must have affected streamflow recession characteristics.

The objective of this study is to investigate temporal and spatial variations in streamflow recession characteristics driven by climate and landscape changes in the glaciated basin of the YRB. The changes in these characteristics are examined using comparisons and contrasts of streamflow recessions indifferent subbasins and time periods according to the power law. To describe the temporal variability in the recession characteristics under climate warming, the decorrelation method by Dralle et al. (2015) is used to obtain a recession parameter a^′ independent of b from individual recessions over the period of 1980–2015. Then, the a^′ and b are regressed with mean temperature in recession periods for each subbasin of YRB. A sensitivity analysis reveals the effect of climate warming on the recession parameters, recession rates, and storage/recharge–discharge relations of different subbasins of YRB. They show the extent of the nonlinearity in the variation of the streamflow recessions in this glaciated basin of TP.

2 Study region and data

The YRB (28.2–31.2^∘ N, 82.0–94.9^∘ E) is the largest river basin in TP (Fig. 1). The main stem of YRB is formed by major suture zones in southern TP from the collision of the Indian plate and the Eurasian plate. The modern YRB flows along the suture from the west to the east before bending to the south at the eastern Himalayan syntaxes with an average west–east gradient of about 2.63 ‰ (Fig. 1a; Tan et al., 2021). In this study, we selected the upstream area of the great gorge of YRB (main stem about 1100 km, with an area of 2.0 × 10⁵ km²). The elevation of the study area drops drastically from 6234 m in the west to 2030 m in the east (Yao et al., 2007; Wang et al., 2021).

Climate in YRB is heavily influenced by the Indian monsoon in summer and the westerlies in winter (Ren et al., 2018; Tian et al., 2020). From the west to the east of the basin, mean annual temperature varies from −9.3 to 22.0 ^∘C, and the mean annual precipitation from 300 to 1050 mm. Nearly 90 % of the annual precipitation falls during June to September. As a result, the mean annual total streamflow of the entire basin, 289.7 mm is highly unevenly distributed in seasons. The summer streamflow is derived from monsoon rainfall and glacier meltwater. Groundwater accounts for about 55 % of the annual streamflow upstream and 27 % downstream of the YRB (Yao et al., 2021).

Table 1 Information of the data used in this study.

There are four hydrological stations along the main stem of YRB, namely LZ, NGS, YC, and NX, as shown in Fig. 1a. There are two additional hydrological stations, YBJ and LS, located in the major tributaries of the Lhasa River, which originates from the Nyainqêntanglha Mountains north of the YRB (Fig. 1a). Daily streamflow data from 1980 to 2015 are available at these hydrological stations, except for LZ. Accordingly, we divide YRB into five subbasins, with three nested subbasins of NGS, YC, and NX in the main stem of YRB and two subbasins of YBJ and LS in the tributaries of the Lhasa River.

There are four main dams/reservoirs (marked by the purple squares in Fig. 1a) in the YRB above the hydrological station of NX. The reservoirs, ML, ZK, PD, and ZM (Manla, Zhikong, Pangduo, and Zangmu, respectively), were built in 1999, 2003, 2007, and 2014, respectively. The reservoirs ML, PD, and ZM are operated daily, while ZK is operated seasonally. The impact of reservoir regulations on streamflow is minor for the subbasins NGS, YC, and NX in the main stem of the YRB because the reservoirs are operated daily and affect less than 10 % of the areas of the tributaries. In the tributary of the Lhasa River, the subbasin YBJ has no reservoirs while, the subbasin LS has two reservoirs, PD and ZK, which have impacts on streamflow.

Daily gridded data (0.1^∘ × 0.1^∘ spatial resolution) of precipitation (P) and mean surface air temperature (T) during 1980–2015 were provided by the National Tibetan Plateau Data Center (Yang and He, 2019; He et al., 2020; http://data.tpdc.ac.cn, last access: 26 November 2020). The subbasin-averaged P and T are calculated by the geometric mean of the gridded data.

Data of the glacier area and permafrost area, and the normalized difference vegetation index (NDVI) were collected from National Tibetan Plateau Data Center (Table 1; http://data.tpdc.ac.cn, last access: 26 November 2020). Glacier areas are located at altitudes from 3370 to 6460 m above sea level (Fig. 1a). Before 2000, the glacier and permafrost area accounted for 1.88 % and 41.8 % of the YRB area, respectively. These coverage percentages have reduced significantly since 2000 (Table 2). The annual mean NDVI was calculated using the maximum value synthesis method from the Global GIMMS NDVI3g v1 dataset, with a 15 d temporal resolution and $\mathrm{1}/\mathrm{12}$^∘ spatial resolution. The vegetation types are mainly alpine meadow, alpine steppe in the upstream (LZ), alpine shrubs and grasslands in the middle region (LZ-NGS), and alpine grassland and forest in the lower YRB (NGS-NX; Liu et al., 2014).

The annual depth of glacial melt data (G) during 1980–2015 is estimated using the degree day model (Su et al., 2015; Liu and Zhang, 2018). The calculation procedures are detailed in Wang et al. (2021). The annual ALT is estimated using a linear statistical function of the air freezing index (FI_air) described in Xu et al. (2017), where FI_air is calculated according to the cumulative value of daily mean temperature below 0 ^∘C in a year.

After the warm season (June–September), there is little precipitation (Hayashi, 2020) in YRB, and the melting of the snow and glacier is minor due to cold temperatures (< −5 ^∘C; Fig. 3). Subsequently, the flow discharge recedes from September to February of the following year. In this study, we use the daily discharge (Q_t) in this recession period, which is defined as being from 1 October to 15 February of the following year, in our analysis of the recession process.

3 Methodologies

3.1 Detection of changes in annual climate and hydrological series

The Mann–Kendall (MK) method (Mann, 1945; Kendall, 1975) combines a trend-free prewhitening treatment (TFPW-MK; Yue and Wang, 2002) with a Sen slope (Sen, 1968). The TFPW-MK is used in this study to test the temporal trend of annual variations in meteorological and hydrological elements at the specified significance level of α = 0.05.

The Pettitt method is applied to detect the change point (year) of annual hydrological and meteorological variables. the Pettitt method is nonparametric and has been widely used in mutation point detection (Pettitt, 1979; Mallakpour and Villarini, 2016; Wang et al., 2021).

3.2 Streamflow recession analysis

Based on analytical solutions to the Boussinesq equation, the relationship of streamflow (Q in units of mm d⁻¹) and streamflow change ($-\mathrm{d}Q/\mathrm{d}t$; mm d⁻²) in a recession period can be expressed by the power law, as follows (Brutsaert and Nieber, 1977):

$$\begin{array}{}\text{(1)}& -\mathrm{d}Q/\mathrm{d}t=a{Q}^b,\end{array}$$

where a (${\mathrm{mm}}^{\mathrm{1}-\mathrm{b}}{\mathrm{d}}^{\mathrm{b}-\mathrm{2}}$) and b (dimensionless) are the recession coefficients (Brutsaert and Nieber, 1977; Tashie et al., 2020).

Based on Eq. (1), the relationship between groundwater storage (S) and streamflow can be derived, as follows:

$$\begin{array}{}\text{(2)}& S=K{Q}^m,\end{array}$$

where K = $[a\cdot (\mathrm{2}-b){]}^{-\mathrm{1}}$ (${\mathrm{mm}}^{\mathrm{b}-\mathrm{1}}{\mathrm{d}}^{\mathrm{2}-\mathrm{b}}$), and $m=\mathrm{2}-b$.

The recession timescale (τ) measures the recession rates of individual recessions (Kirchner, 2009) and is defined as follows:

$$\begin{array}{}\text{(3)}& \mathit{\tau }={\displaystyle \frac{\mathrm{d}S}{\mathrm{d}Q}}={\displaystyle \frac{\mathrm{1}}{a{Q}^{b-\mathrm{1}}}}.\end{array}$$

From Eqs. (1)–(3), the storage sensitivity of discharge (λ_S) for the recession curve (Berghuijs et al., 2016) is as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\lambda }}_{\mathrm{S}}={\displaystyle \frac{\mathrm{d}Q/Q}{\mathrm{d}S}}={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }Q}}=a{Q}^{b-\mathrm{2}},\end{array}$$

where λ_S (mm⁻¹) is a measure of the sensitivity of instantaneous discharge values to water storage changes and indicates the fractional increase in discharge for each unit of increase in storage. The larger (or smaller) the values of a (or b) are, the more sensitive the discharge is to water storage.

Both the relationships of $(-\mathrm{d}Q/\mathrm{d}t)$–Q and S–Q are linear if b=1 and nonlinear if b≠1. When b≠1, the discharge recession is as follows:

$$\begin{array}{}\text{(5)}& Q_{\mathrm{t}}=Q_{\mathrm{0}}{\left(\mathrm{1}+Q_{\mathrm{0}}^{b-\mathrm{1}}a(b-\mathrm{1})t\right)}^{\mathrm{1}/(\mathrm{1}-b)},\end{array}$$

where Q₀ and Q_t are the initial (t=0) discharge and discharge at time t. For any specific initial discharge Q₀, the larger the b is, the faster the hydrograph recession is for high discharge and the more stable the recession is for low discharge (Tashie et al., 2020).

The parameters a and b can be determined by fitting the daily observation data of ($\mathrm{\Delta }Q/\mathrm{\Delta }t)\sim Q$ in a log–log diagram using linear least squares regression. The fitted values of a and b are used to estimate $-\mathrm{d}Q/\mathrm{d}t$ and Q_t using Eqs. (1) and (5). The accuracy of the estimated $-\mathrm{d}Q/\mathrm{d}t$ values is evaluated by the root mean square logarithmic error (RMSLE) as follows (Bekele and Nicklow, 2007):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \text{RMSLE}=\\ & {\left[{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left(\text{log}\right(-\mathrm{d}{Q}_{\text{est}}\left(i\right)/\mathrm{d}t)-\text{log}(-\mathrm{d}{Q}_{\text{obs}}\left(i\right)/\mathrm{d}t){)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where Q_obs.(i) and Q_est.(i) are the observed and estimated discharges, respectively. The terms $\mathrm{d}{Q}_{\text{est}}\left(i\right)/\mathrm{d}t$ and $\mathrm{d}{Q}_{\text{obs}}\left(i\right)/\mathrm{d}t$ in Eq. (6) are derived from taking the derivatives of Eq. (1). In practice, their finite difference $\mathrm{\Delta }Q/\mathrm{\Delta }t$ is determined from the observed recession segments ΔQ in time interval Δt. N in Eq. (6) is the number of data points of $-\mathrm{d}Q/\mathrm{d}t$ in individual recessions.

According to Dralle et al. (2017), by minimizing the following:

$$\begin{array}{}\text{(7)}& E_{\text{MAP}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|{\displaystyle \frac{Q_{\text{obs,}i}-Q_{\text{est},i}}{Q_{\text{obs,}i}}}\right|,\end{array}$$

where E_MAP is the absolute relative error between Q_obs.i and Q_est.i over the recession period, a fitting in each recession hydrograph would ensure that the estimated volume of recession discharge approaches to that of the observed recession.

To avoid recession-scale parameter dependence, we use the method of Dralle et al. (2015) and rescale the discharge Q by $Q=k\widehat{Q}$. The revised power law for the rescaled discharge $\widehat{Q}$ is as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}\widehat{Q}}{\mathrm{d}t}}=-a{k}^{b-\mathrm{1}}{\widehat{Q}}^b=-a^{\prime }{\widehat{Q}}^b,\end{array}$$

where k is a constant, and $a^{\prime }=a{k}^{b-\mathrm{1}}$ is a new recession parameter independent of b. The unit of a^′ is days (Dralle et al., 2015).

Furthermore, to minimize the correlation of the fitted recession exponent and log-transformed fitted recession-scale parameters for a unique value of k, we use the following equation to compute the scaling factor k, according to Bergner and Zouhar (2000), as follows:

$$\begin{array}{}\text{(9)}& k=\text{exp}\left(-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(b_i-\stackrel{\mathrm{‾}}{b}\right)\left(\text{log}\left(a_i\right)-\stackrel{\mathrm{‾}}{\text{log}\left(a\right)}\right)}{{\sum }_{i=\mathrm{1}}^n{\left(b_i-\stackrel{\mathrm{‾}}{b}\right)}^{\mathrm{2}}}}\right).\end{array}$$

In Eq. (9), $\stackrel{\mathrm{‾}}{b}$ and $\stackrel{\mathrm{‾}}{\text{log}\left(a\right)}$ are the arithmetic mean of annually fitted recession exponents $\mathit{\{}{b}_{\mathrm{1}},b_{\mathrm{2}},\mathrm{\dots },b_n\mathit{\}}$ and the log-transformed fitted recession intercepts $\mathit{\left\{}\text{log}\right(a_{\mathrm{1}}),\text{log}(a_{\mathrm{2}}),\mathrm{\dots },\text{log}(a_n\left)\mathit{\right\}}$, respectively, and i is the number of annual values from 1980 to 2015.

3.3 Changes in recession characteristics under warming climate

In cold climate regions, recession coefficients are closely related to the thickness of the active layer in the soil profile above the permafrost layer (Bense et al., 2012; Brutsaert and Hiyama, 2012). Changes in these catchment properties depend on daily, seasonal, and annual temperature variability. For example, the transition from unfrozen to frozen ground for temperatures varying between 0 and −0.5 ^∘C coincides with a reduction in hydraulic conductivity of several orders of magnitude for saturated porous media (Burt and Williams, 1976). On the other hand, when the temperature rises, the thawing front in the active soil layer moves progressively downward as summer proceeds, leading to increased water storage in the active layer. If the frozen soil beneath the thawing front is ice saturated (thus relatively impermeable), the active soil layer can function as a very shallow perched aquifer that controls streamflow response to snowmelt and summer precipitation (Carey and Woo, 2005; Yamazaki et al., 2006; Wright et al., 2009; Koch et al., 2014). To include these variations in a study, recession coefficients are considered to vary, depending on the watershed state and temperatures (Tashie et al., 2019).

In this study, the variability in the parameters a^′ and b are expressed as a function of temperature (T), i.e., a^′(T) and b(T) as follows:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}{a}^{\prime }\left(T\right)=\mathit{\alpha }\cdot \exp \left({\mathit{\alpha }}_{\mathrm{1}}T\right)\\ b\left(T\right)=\mathit{\beta }\cdot \exp \left({\mathit{\beta }}_{\mathrm{1}}T\right),\end{array}\right.\end{array}$$

where α and α₁ are coefficients for a^′(T), β and β₁ are coefficients for b(T), and T is the mean surface air temperature in the recession period (T_re). These coefficients can be obtained by fitting the recession parameters of individual recession events with known T_re in each subbasin.

From Eq. (2), the temporal change in S (ΔS) can be described as follows:

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }S={\displaystyle \frac{\partial S}{\partial K}}\mathrm{\Delta }K+{\displaystyle \frac{\partial S}{\partial m}}\mathrm{\Delta }m+{\displaystyle \frac{\partial S}{\partial Q}}\mathrm{\Delta }Q.\end{array}$$

Incorporating $K=k^{b-\mathrm{1}}/[a^{\prime }\cdot (\mathrm{2}-b\left)\right]$ and $m=\mathrm{2}-b$, we can write ΔS as follows:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathrm{\Delta }S& ={\displaystyle \frac{\partial S}{\partial K}}{\displaystyle \frac{\partial K}{\partial a^{\prime }}}\mathrm{\Delta }{a}^{\prime }+\left({\displaystyle \frac{\partial S}{\partial K}}{\displaystyle \frac{\partial K}{\partial b}}+{\displaystyle \frac{\partial S}{\partial m}}{\displaystyle \frac{\partial m}{\partial b}}\right)\mathrm{\Delta }b\\ & +{\displaystyle \frac{\partial S}{\partial Q}}\mathrm{\Delta }Q\end{array}\end{array}$$

or

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{\Delta }S=& {\mathit{\lambda }}_{a^{\prime }}\mathrm{\Delta }{a}^{\prime }+{\mathit{\lambda }}_b\mathrm{\Delta }b\\ & +{\mathit{\lambda }}_Q\mathrm{\Delta }Q,\end{array}\end{array}$$

where ${\mathit{\lambda }}_{a^{\prime }}$, λ_b, and λ_Q are the sensitivity coefficients of ΔS to a^′, b, and Q, respectively, and can be derived as

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{\mathit{\lambda }}_{a^{\prime }}=\frac{\partial S}{\partial K}\frac{\partial K}{\partial a^{\prime }}=\frac{-k^{b-\mathrm{1}}}{{a^{\prime }}^{\mathrm{2}}(\mathrm{2}-b)}{Q}^{\mathrm{2}-b}\\ {\mathit{\lambda }}_b=\frac{\partial S}{\partial K}\frac{\partial K}{\partial b}+\frac{\partial S}{\partial m}\frac{\partial m}{\partial b}=\frac{k^{b-\mathrm{1}}}{a^{\prime }(\mathrm{2}-b)}{Q}^{\mathrm{2}-b}\left(\frac{\mathrm{1}}{\mathrm{2}-b}-\ln Q\right)\\ {\mathit{\lambda }}_Q=\frac{\partial S}{\partial Q}=\frac{k^{b-\mathrm{1}}}{a^{\prime }{Q}^{b-\mathrm{1}}}=\mathit{\tau }\end{array}\right..\end{array}$$

Because the recession parameters a^′ and b are functions of T in Eq. (8), S is a function of T and Q. The change in S (ΔS) can therefore be expressed as follows:

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }S\left(T\right)={\displaystyle \frac{\partial S}{\partial T}}\mathrm{\Delta }T+{\displaystyle \frac{\partial S}{\partial Q}}\mathrm{\Delta }Q\end{array}$$

or

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }S\left(T\right)={\mathit{\lambda }}_T\mathrm{\Delta }T+{\mathit{\lambda }}_Q\mathrm{\Delta }Q,\end{array}$$

where λ_T is the sensitivity coefficient of ΔS to T. λ_T can be derived as follows:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathit{\lambda }}_T& ={\displaystyle \frac{\partial S}{\partial T}}={\displaystyle \frac{k^{b-\mathrm{1}}}{a^{\prime }(\mathrm{2}-b)}}{Q}^{\mathrm{2}-b}\\ & \times \left({\displaystyle \frac{b({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\alpha }}_{\mathrm{1}})-\mathrm{2}{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}-b}}-b{\mathit{\beta }}_{\mathrm{1}}\ln Q\right).\end{array}\end{array}$$

Figure 2 Variations in (a) annul mean temperature (T), (b) mean temperature in a recession period (T_re), (c) precipitation (P), (d) discharge (Q), (e) discharge in a recession period (Q_re), (f) glacier meltwater (G), (g) NDVI, (h) active layer thickness (ALT), (i) percentage of permafrost area (PPA), (j) the total number of days with the mean temperature above 0 ^∘C in a year (MTD_a), and (k) the total number of days with the mean temperature above 0 ^∘C for recession period (MTD_re) from 1980 to 2015 in the five subbasins. The subscripts 1 and 2 refer to the early period from 1980 to 1996 and the recent period from 1997 to 2015, respectively. NGS–YC refers to the area between the nested subbasins NGS and YC, and YC–NX is for the area of the nested subbasins YC and NX.

4 Results

4.1 Spatial and temporal variations in climate and hydrological variables

4.1.1 Spatial variations

The mean values of the observed climate and hydrological variables during 1980–2015 in the subbasins have shown that the region's climate has become warmer and wetter (Fig. 2a–2c and Table 2). The wet trend is largest in the YBJ and LS subbasins. We note that the mean temperature of YBJ is lowest among the subbasins because most of its area is at higher altitudes. It has the largest fractional glacier coverage area of about 10 %. The percentage of glacier area in the other subbasins is 1.63 %, 1.52 %, 1.92 %, and 0.75 % in NGS, YC, NX, and LS subbasins, respectively. The percentage of permafrost area (PPA) ranges from 41.8 % to 47.7 % in the five subbasins.

The mean annual streamflow (Q) and streamflow in the recession period (Q_re) (both in units of mm) increased (Fig. 2d and e), resulting in higher runoff coefficient Rc $(R/P)$ towards the wetter downstream of the main stem of the YRB. However, the wettest subbasins of YBJ and LS do not have the greatest discharge and high Rc, possibly because of the icy environment. The daily coefficient of variation (CV) in streamflow is higher in YBJ and LS and upstream of NGS. CV decreases towards the wetter downstream of YRB, partially because of strengthened watershed regulation as the subbasin areas increase, and the dams are included in the area of the analysis.

Table 2 Summary of subbasin characteristics in the YR basin.

^∗ Glacier area is from China's second glacier catalogue data in 2009. The mean of A_ice in the years 1976 and 2001 is used as a reference value in the subperiod before 1997, and the mean of A_ice in 2001, 2009, and 2013 is used as a reference value in the subperiod after 1997. The subscripts 1, and 2 represent 1980–1996 and 1997–2015, respectively.

4.1.2 Annual variations in climate and hydrological variables during 1980–2015

Figure 2a–f show variations in the observed climate and hydrological variables from 1980 to 2015. Tested by TFPW-MK, annual mean temperature (T) and temperature in the recession period (T_re) rose significantly (p < 0.05) in all subbasins of the YRB (Fig. 2a and b). Annual T rose at a rate of 0.045–0.075 ${}^{\circ }\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ for the five subbasins, which is smaller than the rate of 0.070–0.097 ${}^{\circ }\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ for annual T_re. Annual precipitation (P) also increased, and was significant (p < 0.1), in the subbasins of YC, NX, and LS, but insignificant in NGS and YBJ (Fig. 2c). Annual glacier meltwater G and the total number of days with the mean temperature above 0 ^∘C in a year (MTD_a) increased significantly in all subbasins (Fig. 2f and j). Meanwhile, the total number of days with the mean temperature above 0 ^∘C in the recession period (MTD_re) also increased, and it was significant in the main stream of the YRB (e.g., NGS, YC, and NX) and insignificant in the two subbasins of YBJ and LS (Fig. 2k). The rate of G increase is from 0.46 mm a⁻¹ in LS to 2.86 mm a⁻¹ in YBJ. The rate of increase in annual MTD_a and MTD_re is 0.48–0.82 and 0.12–0.32 d a⁻¹, respectively, in the five subbasins. Under the warmer and wetter climate, the vegetation coverage tended to increase. The NDVI increased at a rate of 5.1⁻⁴–8.03⁻⁴ a⁻¹ during 1980–2015 in the four subbasins YC, NX, NGS, and YBJ but decreased at a rate of −1.03⁻⁴ a⁻¹ in LS (Fig. 2g).

Over the same period, the annual mean discharge Q and the mean discharge in the recession period Q_re increased significantly, except for Q_re in the LS subbasin (Fig. 2d an e). The Q increased at a rate of 1.16–1.68 mm a⁻¹, which is higher than the rate of 0.22–0.47 mm a⁻¹ for Q_re. In contrast, Q_re in LS decreased insignificantly (Fig. 2d), possibly because of initial water storage when the ZK reservoir (Fig. 1a) began operating in around 2007.

In Wang et al. (2021), a point of dramatic change was detected around 1997 by the Pettitt test for T, G, and Q and around 1995 for annual P in the subbasins of NGS, YC, and NX. In this study, the same point of change in 1997 was detected in the annual series of T_re and Q_re in the NGS, YC, and NX subbasins and in the annual series of T, G, and Q in the YBJ and LS subbasins. The annual time series of Q in LS was detected to have two additional points of change in 1995 and 2005, which is possibly attributable to the increased impact of human (reservoir operation) activities (Cai et al., 2021). The point of change in 1997 has also been identified by the dramatic changes in surface conditions, i.e., the reversed NDVI trend (Fig. 2g), increased ALT (Fig. 2h), and accelerated thawing of permafrost after 1997 (Fig. 2i). Accordingly, we separate the study years from 1980–2015 into two periods, i.e., the early period from 1980 to 1996 and the recent period from 1997 to 2015.

Figure 3 (a–e) P, T, and Q in a hydrological year (from 1 March to 28 February of the following year) for the two periods in the five subbasins. The red dashed rectangle in (a) shows the hydrograph recession from 1 October to 15 February of the following year, and the shading shows the range of the daily variation in P, T, and Q in each period.

Climate in the recent period has changed to be markedly warmer and wetter. The mean annual P after 1997 increased by 27–46 mm or 7.9 %–10.7 % compared to that in the early period in the five subbasins. The mean annual T in the recent period is 0.75–1.52 ^∘C warmer than that in the early period, and a larger rise of 1.40–1.78 ^∘C was found for the mean T_re after 1997. These changes concurred at 8.6–55.9 mm or a 23.8 %–81.1 % increase in mean annual G after 1997 in the five subbasins. The mean MTD_a and MTD_re in the recent period is 8–18 and 2–7 d greater than that in the early period of 1980–1996, respectively. As a result, Q increased by 29.6–50.2 mm or 12.7 %–31.5 % in the recent period compared to that before 1997. This increase in Q and Q_re after 1997 is much larger in the upstream subbasins NGS and YC.

Figure 4 Discharge recession for selected years with approximately the same initial discharge Q₀ in each subbasin.

4.1.3 Annual recession characteristics

As shown in Fig. 3, the annual hydrographs in the five subbasins are consistent, delineating a single peak response to maximum precipitation and temperature in July–August. The statistic values of the daily discharge series in the two periods show that the mean value in the recent period exceeds that in the early period in all subbasins. Meanwhile, daily discharge variability in the recession of the annual hydrograph in the recent period is also greater than that of the early period, shown by larger CV after 1997 for most subbasins, except subbasin YBJ (Table 2).

Figure 4 further exhibits that changes in recession rates are different in the early phase and the later phase of the recession for the selected hydrographs with approximately the same initial discharge Q₀ in each subbasin. In the recent period, the streamflow recedes faster in the early phase of the recession and slows down in the later phase, except for LS.

Figure 5 (a–e) Plot of $-\mathrm{d}Q/\mathrm{d}t$ vs. Q in log–log diagram for each recession hydrograph during 1980–2015 and the fitting lines [$\text{log}(-\mathrm{d}Q/\mathrm{d}t)=b\text{log}\left(Q\right)+\text{log}\left(a\right)$] for the data points in the two periods for the five subbasins. (f) Differences in mean recession rates between the two periods (Δv_Q) estimated from the non-overlapping moving averages of the 5 d series.

The faster recessions in the recent periods are also illustrated by the regression on data points of $-\mathrm{d}Q/\mathrm{d}t$ vs. Q in the log–log diagram (Fig. 5). The fitted line of $-\mathrm{d}Q/\mathrm{d}t$ vs. Q after 1997 has a steeper slope (Fig. 5a–e) and a more negative intercept, except for LS, indicating a higher value of b and a smaller value of a for the subbasins of NGS, YC, NX, and YBJ after 1997. The increased b value and decreased a value after 1997 suggest that the recession curves tend to be more concave for these subbasins when the climate is warmer and wetter, as indicated in Fig. 4. According to the non-overlapping moving averages of the 5 d series of the recession discharge, the estimated average recession rate after 1997 [$(\mathrm{\Delta }Q/\mathrm{\Delta }t{)}_{\mathrm{2}}$] is larger than that before 1997 [$(\mathrm{\Delta }Q/\mathrm{\Delta }t{)}_{\mathrm{1}}$], as indicated by the positive values of $\mathrm{\Delta }{v}_Q=(\mathrm{\Delta }Q/\mathrm{\Delta }t{)}_{\mathrm{2}}-(\mathrm{\Delta }Q/\mathrm{\Delta }t{)}_{\mathrm{1}}$ in Fig. 5f.

Table 3 Mean values of parameters a (${\mathrm{mm}}^{\mathrm{1}-\mathrm{b}}{\mathrm{d}}^{\mathrm{b}-\mathrm{2}}$), a^′ (d⁻¹) and b (dimensionless), recession coefficient K (${\mathrm{mm}}^{\mathrm{b}-\mathrm{1}}{\mathrm{d}}^{\mathrm{2}-\mathrm{b}}$), and recession timescale τ (d).

The asterisk^∗ indicates the significance by the TFPW-MK test (p < 0.05). Values in parentheses refer to the range of annual value.

Figure 6 The box plot of (a) the RMSLE of $-\mathrm{d}Q/\mathrm{d}t$ between observed and estimated discharge and (b) E_MAP between annual mean observed and estimated discharge for individual recession hydrographs in each subbasin.

4.2 Estimation of the recession parameters

We used the observed hydrograph in each year and obtained the recession parameters a and b. Afterwards, a^′ is calculated using the scaling factor k (0.527, 0.602, 0.740, 0.594, and 0.611 mm d⁻¹ for NGS, YC, NX, YBJ, and LS, respectively). The results are summarized in Table 3. The mean values of annual a and a^′ during 1980–2015 range from 0.022 to 0.042 ${\mathrm{mm}}^{\mathrm{1}-\mathrm{b}}{\mathrm{d}}^{\mathrm{b}-\mathrm{2}}$ and from 0.015 to 0.025 d⁻¹, respectively, and the mean annual value of b ranges from 1.36 to 1.85 for the five subbasins. Figure 6 shows the errors of the estimated recession (RMSLE and E_MAP) for the recession curve of each subbasin. Mean annual RMSLE is less than 0.15 in all subbasins. The mean annual E_MAP is lower than 10 %, except in the subbasins of YBJ and LS, where E_MAP is 0.15 and 0.14, respectively.

Figure 7 (a–e) The exponential function of a^′ with T_re for each subbasin. (f–j) Same as panels (a–e) but for recession slope b. The solid and open circles represent the 4-year average and annual value, respectively.

Figure 7 shows the relationship of the annual value and 4-year moving average value of the recession parameters a^′ and b with mean surface air temperature in the recession period (T_re) in each subbasin. The exponential function between a^′ or b and T_re is fitted with a high determination coefficient (R² ranges 0.51–0.67 for a^′ and 0.58–0.87 for b). These results show that a^′ decreases exponentially with increasing T_re in the subbasins, except for LS, while b increases exponentially with increasing T_re in all subbasins.

For the multiyear mean values of parameters a, a^′, and b in the two periods (Table 3), the mean value of a and a^′ in the recent period ranges between 0.021–0.039 ${\mathrm{mm}}^{\mathrm{1}-\mathrm{b}}{\mathrm{d}}^{\mathrm{b}-\mathrm{2}}$ and 0.014–0.023 d⁻¹, respectively, which is smaller than that in the early period when they range between 0.022–0.046 ${\mathrm{mm}}^{\mathrm{1}-\mathrm{b}}{\mathrm{d}}^{\mathrm{b}-\mathrm{2}}$ and 0.015–0.027 d⁻¹, respectively, in the subbasins, except for LS. On the other hand, b ranges between 1.47–1.89 in the recent period, which is larger than that in the early period (1.25–1.81) for all subbasins. These results indicate that a warming and wetting climate increases the nonlinearity of the recession (b) and enhances stability of low streamflow [log(a)] for most subbasins in YRB.

Figure 8 Relationship of S and Q at different values of K (a) and m (b) and combinations of K and m (c) in the two periods for subbasin YC. The different set of numbers in the blue box with a dashed line in each panel shows the K and m values and therefore the different S–Q relationships in the recent warmer period.

4.3 Change in storage–discharge relationship in warming climate

The strong sensitivity of the recession parameters of a (a^′) and b to T means that a warming climate can change the nonlinear relationship of the water storage (S) and discharge/streamflow (Q) (Eq. 2). This change is shown by the increase in the recession coefficient K in the recent period (11.5–73.2 ${\mathrm{mm}}^{\mathrm{b}-\mathrm{1}}{\mathrm{d}}^{\mathrm{2}-\mathrm{b}}$ for ΔK in Table 3) and decrease in m ($=\mathrm{2}-b$) in Eq. (2) because b increases with the temperature (Table 3 and Fig. 7). The increase in K and decrease in m mean a lower discharge for a specific storage or a higher storage for a specific discharge in the recent warmer period. As an example, we show, in Fig. 8, the relationship of S with Q for different K and/or m between the two periods in subbasin YC. Q decreases significantly with the increase in K or decrease in m for storage S (Fig. 8a and b). The combination of the increase in K and decrease in m leads to a marked decrease in Q for storage S (Fig. 8c). Correspondingly, the recession timescale (τ = $\mathrm{d}S/\mathrm{d}Q$) increases by 5.2–20.7 d in the recent warmer period in all subbasins (Table 3), especially the glaciated subbasin YBJ. The increase in τ also means an increase in storage for Q because $S\left(Q\right)S_{\mathrm{0}}={\int }_{S_{\mathrm{0}}}^S\mathrm{d}S=-{\int }_{Q_{\mathrm{0}}}^Q\mathit{\tau }\left(Q\right)\mathrm{d}Q$ > 0 in the recession.

Table 4 The storage sensitivity of discharge (λ_S; mm⁻¹) and sensitivity coefficients of recession parameters of a^′ (${\mathit{\lambda }}_{a^{\prime }}$) and b (λ_b), T (λ_T), and Q (λ_Q) to storage changes (ΔS) during different periods for each subbasin. The values in parentheses refer to the range of annual value.

The asterisk ^* indicates the significance by the TFPW-MK test (p<0.05).

The lower discharge for a specific storage or higher storage for any specific discharge are further shown by the results of the storage sensitivity of the discharge (λ_S in Eq. 4) in Table 4. The mean value of λ_S during 1980–2015 ranges between 0.036 and 0.059 mm⁻¹ for the five subbasins. These values mean that 1 mm decrease in storage only results in 3.6 %–5.9 % decrease in discharge. Thus, in a warmer climate, a unit decrease in the storage releases less water to discharge in the recession period and is supported by the decrease in λ_S in the recent warmer period (see the negative values of Δλ_S in Table 4) and towards the warmer and wetter downstream for the main stem of YRB (from NGS to NX; see the mean annual value of λ_S during 1980–2015 in Table 4). This is especially so in glaciated basins, e.g., YBJ, where the decrease in λ_S (Δλ_S in Table 4) is largest in the recent period, corresponding to the largest increase in τ. We also note that the decrease in λ_S (Δλ_S) in LS is relatively small in the recent period, primarily because of the regulation of reservoirs on discharge.

4.4 Sensitivity of the recession parameters to storage change under climate warming

The change in the sensitivity of discharge to storage (λ_S) should affect the recession processes. This effect is described by the sensitivity of the recession parameters to storage change ΔS (${\mathit{\lambda }}_{a^{\prime }}$ and λ_b for a^′ and b, respectively). As listed in Table 4, λ_b is positive, while ${\mathit{\lambda }}_{a^{\prime }}$ is negative in all five subbasins. The larger the negative ${\mathit{\lambda }}_{a^{\prime }}$ is, the larger the positive λ_b would be in the recent warmer period, particularly in the glaciated subbasin of YBJ ($\mathrm{\Delta }{\mathit{\lambda }}_{a^{\prime }}$ and Δλ_b in Table 4). Thus, the increase in storage in a warmer climate will correspondingly enhance the nonlinearity of recession (b) and streamflow stability [log(a)] in the YRB, as indicated in Eq. (13). This sensitivity can be weakened, however, by an anthropogenic effect (reservoir regulations), as suggested by the less sensitive result of the recession parameters of a^′ and b to ΔS in LS subbasin.

Effects of climate warming on the change in the storage and streamflow recession characteristics can be further shown by the storage change ΔS in response to changes in T and Q (ΔT and ΔQ) using Eq. (16). The values of λ_T in Table 4 show an increase in the recent period in all subbasins, especially in YBJ and YC. Thus, the enlarged storage is largely attributed to climate warming. As expected, λ_T is smaller in the subbasin LS with its reservoir regulation. The values of λ_Q in Table 4 also become bigger in the recent warmer period in all subbasins. These changes indicate that climate warming increases storage and discharge.

Figure 9 Changes in storage ΔS in relation to changes in discharge ΔQ under different changes in temperature ΔT for each subbasin. The changes in temperature ΔT refer to the annual values relative to mean annual temperature in the recession period during 1980–2015. The solid circle refers the point of ΔS in response to 0.2 mm of ΔQ.

However, the increase in discharge ΔQ in response to the increase in storage ΔS can be quite different in response to the different rate of change in temperature ΔT in the five subbasins. According to Eq. (16), the relationship between ΔS and ΔQ for different ΔT in the five subbasins is shown in Fig. 9. As the temperature rises, ΔS becomes greater, while the greater increase in storage volume (in thawing soil layers) allows a smaller amount of water to be released as baseflow. For example, as the mean annual temperature rises by 1.2 ^∘C (ΔT = 1.2 ^∘C), the increase in discharge (ΔQ =  0.2 mm d⁻¹) corresponds to a storage increase of about 92 mm in NGS, 160 mm in YC, 63 mm in NX, 42 mm in LS, and a huge increase of 478 mm in the YBJ (see the value of the ΔS vs. ΔQ in Fig. 9). These results suggest that a larger increase in water storage caused a smaller increase in baseflow in the glaciated subbasins, reflecting a buffering effect of freezing on streamflow dynamics.

When the T_re increased from the early to the recent period, from −7.17 to −5.39 ^∘C in NGS, −6.84 to −5.16 ^∘C in YC, −6.54 to −5.02 ^∘C in NX, −9.20 to −7.69 ^∘C in YBJ, and −7.60 to −6.20 ^∘C in LS, the estimated water storage in terms of λ_TΔT in Eq. (16) increased between 15.2 and 132.6 mm for the five subbasins. These increases contribute to about 86.4 %–99.9 % of the total increase in storage (${\mathit{\lambda }}_T\mathrm{\Delta }T/\mathrm{\Delta }S)$ in those subbasins. They only cause a 0.1 %–13.6 % increase in discharge in terms of ${\mathit{\lambda }}_Q\mathrm{\Delta }Q/\mathrm{\Delta }S$ in Eq. (16). This relationship varies among the five subbasins with different glaciate conditions. In the warm and wet subbasin of NX with low glacial coverage, the increase in storage (λ_TΔT) is relatively small (35.4 mm and 86.4 % of the total increase in storage), and the increase in discharge (λ_QΔQ) is relatively large (5.56 mm and 13.6 % of the total increase in storage). In the cold subbasins with high glacial coverage, climate warming causes a large increase in storage but a small increase in discharge. As an example, the YBJ subbasin has 97.3 % of the total increase in storage vs. only 2.7 % of the total increase in discharge. Again, this relationship is distorted in basins with human regulatory actions in water management. In the subbasin of LS, with strong reservoir regulations, changes in both the storage and discharge are small (e.g., 0.1 % in the total increase in storage).

Figure 10 A schematic illustration of the climate warming effect on surface conditions and the subsurface profile and hydrological variables. The larger sizes of the arrows indicate a large increase in the hydrological variables, e.g., glacier melting, precipitation, discharge, and evaporation.

5 Discussion

Observations have shown that climate warming has accelerated glacier melting and permafrost thawing in cold climate and high-altitude regions. Subsequent changes are found in vegetation growth and the thickening of talik and active soil layer thickness. These changes have altered land surface conditions and unconsolidated soil profiles and subsurface permafrost and thereby redefined the surface and groundwater exchange and balance in those regions (Fig. 10). Our case study of YRB in southern TP shows that accelerated glacier melting and permafrost thawing during 1980–2015 have substantially increased its dynamic groundwater storage, which is defined as $S\left(Q\right)-S_{\mathrm{0}}=-{\int }_{Q_{\mathrm{0}}}^Q\mathit{\tau }\left(Q\right)\mathrm{d}Q$. These results, with the decrease in terrestrial water storage (TWS) in the southern TP, including YRB (Wang et al., 2020) in recent decades, indicate a transformation of water storage in the region from a solid form (glacier and permafrost) to a liquid volume (soil moisture, surface water in rivers/lakes, and groundwater; Fig. 10b). According to the water balance in a catchment, i.e., $\mathrm{d}S/\mathrm{d}t$ = $P_{\mathrm{r}}-E-Q$, where S is regarded as the liquid volume (here, the change in S is equal to the sum of changes in soil moisture and groundwater), E is evapotranspiration, and P_r is the recharge from glacier melting, permafrost thawing, and precipitation, the increase in S infers that P_r is larger than the sum of E and Q in a study region. Because cold regions tend to have a greater coverage percentage of glacier and permafrost, glacier melting and permafrost thawing could substantially increase water storage under climate warming. Higher water storage could extend the recession period and sustain healthy annual streamflow.

Our study also shows that the increase in water storage and its effect on the annual recession of streamflow weakened towards the warmer downstream areas of YRB (with diminishing glacier melting and permafrost thawing effect). Accordingly, if the climate warming continues, the shrinking of the glacier and permafrost volume could eventually reach a point when there is not enough melting to recharge the liquid volume of water in YRB. From that point onward, steady streamflow in YRB would be in danger.

While the processes initiated by the accelerated glacier melting and permafrost thawing extend subsurface flow paths (Hinzman et al., 2020) and the streamflow recession time (τ), the increase in surface temperature and E can also increase surface water loss. According to the discharge relation $-\mathrm{d}Q/\mathrm{d}t$ = $-\frac{\mathrm{d}Q}{\mathrm{d}S}\frac{\mathrm{d}S}{\mathrm{d}t}$ = $(-P_{\mathrm{r}}+E+Q)/\mathit{\tau }$ (Kirchner, 2009), where P_r can be neglected in the recession period, a faster recession ($-\mathrm{d}Q/\mathrm{d}t$) could occur under climate warming from a faster decrease in storage ($\mathrm{d}S/\mathrm{d}t$) due to increasing E (Tashie et al., 2020), and a faster response of discharge to storage ($\mathrm{d}Q/\mathrm{d}S$) could occur due to the increase in the effective hydraulic properties (Lamontagne-Hallé et al., 2018). This phenomenon has happened for the initial recession period of 1–2 months when temperature is relatively higher (see Fig. 4).

Meanwhile, the accelerated glacier melting and permafrost thawing have shortened the frost period (i.e., the prolonged MTD_a and MTD_re in Fig. 2j and k) and thereby increased the soil active layer thickness (ALT) for groundwater storage (Lamontagne-Hallé et al., 2018) and lengthened subsurface flow paths (Hinzman et al., 2020) and the streamflow recession time (τ). These changes weaken the sensitivity of the discharge to the storage ($\mathrm{d}Q/\mathrm{d}S$ = $\mathrm{1}/\mathit{\tau }$) and slow down the recession rate ($-\mathrm{d}Q/\mathrm{d}t$), which is found in the later recession period when temperature is relatively lower (see Fig. 4). Therefore, the competing effects from the warming climate on $\mathrm{d}S/\mathrm{d}t$ and $\mathrm{d}Q/\mathrm{d}S$ in different recession phases would lead to the recession curve becoming more concave, which indicates an increase in the nonlinearity of the recession (b) and streamflow stability [log(a)] in YRB. In comparison, in the warm climate area, the effect of the storage decrease ($\mathrm{d}S/\mathrm{d}t$) on the recession ($-\mathrm{d}Q/\mathrm{d}t$) strengthens, and the effect of the recession timescale ($\mathrm{1}/\mathit{\tau }$ or $\mathrm{d}Q/\mathrm{d}S$) on the recession weakens. As shown in Fig. 9, when temperature is higher (e.g., large and positive ΔT), hydrograph recession (negative ΔQ) is faster, along with the faster decline of storage (ΔS).

Additionally, deep circulating groundwater through macro structures, such as north–south-oriented active tensile faults (Fig. 1b), could also affect the baseflow and its recharge and discharge (Tan et al., 2021). According to studies using multi-tracer data (e.g., ²H, ³H, ¹⁸O, and Sr), modern meltwater is found to primarily maintain the rapid recharge of phreatic groundwater in alpine regions through faults and fissures (Shi et al., 2021). In the middle of the YRB (i.e., NGS–YC subbasins), changes in storage sensitivity to temperature (λ_Tin Table 4) and recession timescale (τ) are greater than those in the upstream and downstream areas (NGS and NX subbasin, respectively). Rising temperature can greatly increase storage (Fig. 10b).

Finally, the anthropogenic effects from reservoir regulation can reduce the climate warming effect on these storage–discharge responses in the YRB. For example, in the LS subbasin, operations of two reservoirs significantly reduced the sensitivity of the recession parameters a (a^′) and b to climate warming and increased streamflow stability [log(a)] (Fig. 4e). It remains questionable, however, as to how this human effort in water management in YRB would be practical/beneficial after the point when the increase in water storage from glacier and permafrost melt has exhausted the solid volume of water resources in the basin following climate warming.

6 Concluding remarks

Climate warming accelerated after 1997 in the Tibetan Plateau, especially in its cold and high-altitude upstream areas. Since 1997, the mean annual temperature has risen by 0.75–1.52 ^∘C, and the mean temperature in the annual recession period (1 October–15 February of the following year) has risen by 1.40–1.78 ^∘C in the five subbasins of the YRB. The largest rise in temperature occurred in the drier and colder subbasins in the upstream YRB. The recent strong warming has accelerated glacier melting and permafrost thawing, and thereby increased annual streamflow (12.7 %–31.5 % larger than the mean value in the early period before 1997) and streamflow in the recession period (20.9 %–25.8 % larger than before 1997) for the five subbasins, except LS where reservoir operations are active and heavily affecting the streamflow. These processes initiated by climate warming have changed the hydrological properties of subbasins considerably and altered the recession characteristics and the storage–discharge relationships.

We have found that the recession parameter a (a^′) that characterizes the stability of streamflow has decreased exponentially in the subbasins, except for LS. Meanwhile, the parameter b that describes the nonlinearity of the recession to discharge has increased exponentially in all the subbasins. These results indicate that climate warming increases the nonlinearity of the recessions and enhances streamflow stability in most of the subbasins in YRB. Our sensitivity analysis further shows the decrease in the sensitivity of discharge/streamflow to storage under the warming climate. Currently, the accelerated glacier melting and permafrost thawing have recharged the system, deepening the active subsurface zone and increasing groundwater storage. Only a small fraction of the enlarged storage is released in surface streams because the increase in active water layer lengthens subsurface flow paths. These changes have also increased the recession timescale and retarded the recession in the later period, particularly in high-altitude cold climate areas. In the relatively warm climate areas downstream of the YRB, the effect of these changes is minor.

As the liquid water storage has increased greatly from melting glaciers and thawing permafrost in the YRB in the recent warming climate, the fast erosion of the solid water storage has weakened its buffering effect of the streamflow, which is becoming less stable and more vulnerable to individual intense precipitation events, such as increase in the recession rate in the initial period found in this study. There are the following two potential consequences from these changes: one is the increase in flash flooding in the trend of rising precipitation in the high-altitude subbasins where more land is exposed after the retreat of glaciers, and the other is the extreme scenario of exhaustion of the water resources in the upstream of the YRB after the buffering effect of glacier and permafrost is lost following the continued warming of the climate.

While human interference with these processes, via reservoirs and regulations, can reduce and curb these impacts of climate warming on storage–discharge relationships, recession characteristics, and streamflow in short term, as shown in the subbasin of LS, long-term strategies need to be developed to not only cope with the short-term needs but also the sustainability of water resources in the Tibetan Plateau under the threat of the continued warming that could change the entire hydrological system in this critical water source region for the one of world's most populated nations.