A Climatic Perspective on the Impacts of Global Warming on Water Cycle of Cold Mountainous Catchments in the Tibetan Plateau: A Case Study in Yarlung Zangbo River Basin

Abstract

Global warming has a profound influence on global and regional water cycles, especially in the cold mountainous area. However, detecting and quantifying such changes are still difficult because noise and variability in observed streamflow are relatively larger than the long-term trends. In this study, the impacts of global warming on the catchment water cycles in the Yarlung Zangbo River Basin (YZRB), one of most important catchments in south of the Tibetan Plateau, are quantified using a climatic approach based on the relationship between basin-scale groundwater storage and low flow at the annual time scale. By using a quantile regression method and flow recession analysis, changes in low flow regimes and basin-scale groundwater storage at the Nuxia hydrological station are quantified at the annual time scale during 1961–2000. Results show annual low flows (10th and 25th annual flows) of the YZRB have decreased significantly, while long-term annual precipitation, total streamflow, and high flows are statistically unchanged. Annual lowest seven-day flow shows a significantly downward trend (2.2 m³/s/a, p < 0.05) and its timing has advanced about 12 days (2.8 day/10a, p < 0.1) during the study period. Estimated annual basin-scale groundwater storage also shows a significant decreasing trend at a rate of 0.079 mm/a (p < 0.05) over the study period. Further analysis suggests that evaporation increase, decreased snow-fraction, and increased annual precipitation intensity induced by the rising temperature possibly are the drivers causing a significant decline in catchment low flow regimes and groundwater storage in the study area. This highlights that an increase in temperature has likely already caused significant changes in regional flow regimes in the high and cold mountainous regions, which has alarming consequences in regional ecological protection and sustainable water resources management.

1. Introduction

Steady rising global temperature has already altered regional water cycles in many places around the world and has posed great challenges to local water resources management [1,2,3]. This issue is relatively more alarming for cold mountainous catchments because of faster warming rates [4], especially in the Asian water tower that is the Tibetan Plateau (hereafter denoted as the TP) [5,6], from where meltwater feeds several large rivers of Asia (e.g., Indus River, Brahmaputra River, Ganges River, Yellow River and Yangtze River) and supports almost half of the world’s population [7]. Many studies have examined the impacts of climate warming on the hydrological process and improved our knowledge of water cycle in this area [8,9,10,11]. However, these investigations basically examined trends in the mean annual flow using various statistical methods and numerical models and provided limited information on the impacts of climate warming on local water cycle due to large natural climate variability comparing to impacts signals and various compensating changes amongst different water cycle components. Previous studies suggested that impact signals can be detected more reliably and effectively from extreme hydrological events and low frequency hydrological variables [12,13], such as the low flow regimes [14,15,16]. Because low flow exhibits catchment long-term signal of climate variability through groundwater storage and thus is ideal for detecting variability and trends of climate at decadal time scale [17,18,19]. However, changes in low flow regimes have rarely been addressed and/or quantified to understand the impacts of global warming in the TP region.

Changes in low flow regimes of a pristine catchment are largely determined by evolution of catchment-scale groundwater storage. Typically, there are three types of data sources for quantifying groundwater storage information including well observation, the Gravity Recovery and Climate Experience (GRACE), and streamflow observations. One can also estimate groundwater storage by numerical modelling. Groundwater monitoring wells can provide most accurate information of regional groundwater levels [20]. However, such observations are usually expensive and sparse outside of industrialized regions [20] and represent only point scale information [19]. Hydrological models can provide dynamic variation of regional groundwater storage at much cheaper expanses and over various spatial and temporal scales. But the accuracy of modelled results remains a challenge due to the uncertainty in model structure, parameterization, and model calibration [21,22]. Recent development in satellite gravity measurement technology (e.g., GRACE), provides a new approach to monitor terrestrial water storage and GRACE data has been widely used to derive regional groundwater storage changes [23,24,25]. Unfortunately, the record of satellite monitoring is relatively short (only from 2002) for detecting changes in regional groundwater dynamics or water cycle as such studies typically require long-term observations (e.g., >30a) to get statistically meaningful conclusions [26]. Based on the catchment storage–discharge relationship or low flow recession theory, Brutsaert [27] proposed a method to infer annual change of catchment groundwater storage from annual lowest seven-day flow, which enables us to estimate decadal trends in catchment groundwater storage from long-term daily streamflow observations that is widely available around the world.

By making use of the methods of Brutsaert [27], the impact of climatic changes in precipitation, temperature, radiation, etc., i.e., forcing of the catchment water cycle, can be detected and quantified from long-term daily streamflow observations via the catchment low flow recession processes. As catchment low frequency hydrological variable, essentially, annual tendency of low flow is an indirect yet reliable indicator of long-term climatic change impact signals in regional water cycle substantiated by the evolution of basin-scale groundwater storage–discharge relationship [19]. In this way, Brutsaert [27] defined it as a climatic perspective for understanding the climatic impacts signals in catchment water cycle. This climatic perspective method for detecting climate impact signal in regional water cycle has been successfully applied to detect variability and long-term trend in catchment groundwater storage under different climatic and land cover conditions around the world including the U.S [28], Australia [29], China [30,31], Japan [32], eastern Siberia [33] and northern Sweden [34]. In the TP region, streamflow observation is possibly the best and the most reliable data for estimating long-term changes in groundwater dynamics as both well and meteorological gauges are scarce and hydrological cycle and groundwater dynamic in such region are very complex for establishing numerical model as well [35]. Therefore, the method of Brutsaert [27] is adopted in this study to detect and quantify the impacts of climate change on the catchment water cycle in the TP region from a climatic perspective based the relationship between basin-scale groundwater storage and low flow.

The Yarlung Zangb River basin (YZRB) is the most important river to understand the hydrological cycle in the TP region. It is not only the largest river system in the TP region taking up about 47% of water resource with less than 10% area in the whole TP region [36], but also is the moisture transportation channel that transferring more than 50% of total precipitation water of the inner region of TP from the Indian Ocean [37,38,39]. Meanwhile, the YZRB is also one of most sensitive regions for global warming as the elevation-dependent warming [4]. Previous studies indicated that the YZRB had experienced faster temperature rising than global average level [40,41], which had caused changes in evaporation [9], precipitation characteristics [42] and meltwater [43]. Therefore, the YZRB is chosen as a case study area to detect the impacts of climate warming on catchment water cycle from climatic perspectives (i.e., long-term annual changes of low flow regimes and groundwater storage). Meanwhile, as one of the largest ecologically sensitive areas around the world [44], low flow regimes of the YZRB and its response to climate change have not been examined, but are indeed important for maintaining the ecosystem health of YZRB and sustainable development in the downstream. An understanding of groundwater storage changes will also improve local water management capacity for stream ecology protection [45], irrigation in the downstream (e.g., Bhutan, Northeast India, and Bangladesh) [46], society development, and adaptation in the YZRB.

To understand the impacts of climate waring on the water cycle in the YZRB, in this study, the method proposed by Brutsaert [27] also referred to as the climatic perspective, is adopted to quantify the impacts signals from changes in catchment long-term annual low flow regimes based on the basin-scale groundwater storage–discharge relationship. By making use of long-term daily streamflow records, climate forcing (precipitation and temperature), remotely sensed snow coverage, underlying mechanisms that are responsible for detecting changes in low flow regimes and groundwater storage dynamics are further investigated and discussed. The primary objectives of this study are: (1) to investigate the changes of the low flow regimes; (2) to estimate groundwater storage changes inferred from long-term daily streamflow records; (3) to detect the signals of climate warming in the local water cycle and to quantify its impacts on the basin-scale groundwater storage.

2. Study Area and Methodology

2.1. Description of Yarlung Zangbo River and Data Set

The Yarlung Zangbo River originates from the Gema Yangzong Glacier in the south-central TP and is one of the highest basin in the world with an average altitude of more than 4000 m (ranging from 132 m to 7258 m) [47]. The study area of this study is the upstream of Nuxia hydrological station, of which drainage area is about 1.8 × 10⁶ km² (Figure 1). Climate of the YZRB is characterized as cold plateau mountain climate with intense solar radiation and low air temperature with annual mean temperature of 0.1 °C. Winter (November to April) is the dry season of the YZRB with strong and dry westerly winds. Precipitation falls dominantly (over 90%) as rainfall in the warm season (May to October) when wet airflow from the Indian Ocean brings a substantial quantity of rainfall to the YZRB.

Grassland, forest, and bare land are the most common land-cover type in the YZRB and account for approximately 94% of the catchment area. Percentage of built-up is less than 0.05%. During the period 1985–2005, only about 1% land cover changed and the YZRB remains a natural watershed [47]. Vegetation cover within the area exhibits district vertical belts and varies from mountain forest, mountain broad-leaved forest, mountain coniferous forest, subalpine shrub meadow to alpine meadow along with raising elevation. Glaciers distribute at the higher elevation belt and its area accounts for about 1.6% of the basin based on the global glacier data set from the Randolph Glacier Inventory (RGI [48]). Permafrost distributes widely in approximately 52% area estimated by the widely used and validated altitude model proposed by Wu, et al. [49].

Daily streamflow observations at Nuxia hydrological station were available over the period 1961–2000 and it controls about 80% drainage areas of the whole YZRB (see Figure 1). Dry season flows were extracted within the drought year following [50], which is defined as from 1 August of last year to the 31 July of current year for the study area. Years with less than 95% availability of observations during the dry season were excluded to minimize the influence of missing data on sample of annual low flow and lowest 7-day flow used for estimating groundwater storage in next section.

Monthly gridded precipitation and temperature data with a spatial resolution of 0.5 degree during the same period as streamflow data were collected from China Meteorological Data Service Center (http://data.cma.cn) to investigate the changes in the meteorological conditions in the YZRB. The gridded precipitation and temperature data were generated from station observations using thin plate spline method. This data had minimized the influence of the undulating terrain on the accuracy of precipitation and temperature space based on the observed meteorological variables of 2472 national-level ground-based meteorological stations in China. Huang, et al. [51] reported that this gridded data set demonstrates higher precision than other global data sets (e.g., Global Land Data Assimilation System (GLDAS), Climatic Research Unit (CRU)) and are also more reliable than meteorological observation from sparse gauging station network within the YZRB.

Monthly snow depth with a spatial resolution of 0.25 degree accessed from the Cold and Arid Regions Sciences Data Center (http://westdc.westgis.ac.cn) [52,53,54] was also used to examine changes in long-term precipitation regimes and to infer the influence of changes in long-term snow fraction on recharge of groundwater storage.

2.2. Quantile Regression Analysis for Investigating Changes in Flow Regimes

The quantile regression analysis, introduced by Koenker and Hallock [55] as a location method to extend ordinary least square, was employed in this study to investigate the changes in the flow regimes in the YZRB. The method has been widely applied for investigating changes in different regimes of daily flow [15,16,27,56,57,58,59]. For example, trend in 25% percentile flows is a least-square regression of the flow that is exceeded 25% daily flows in a given year. In this study, linear trends in 10th and 25th percentile annual flows (i.e., 10% and 25% annual flows) were used to represent the changes in low flow regimes [16]. Meanwhile, linear trends in 50%, 75%, and 90% were also analyzed to compare with low flow regimes for a better understanding of the changes in flow regimes and distributions.

2.3. Groundwater Storage Trend Inferred from Daily Streamflow

In this study, a robust method proposed by Brutsaert [27] is adopted to derive inter-annual variability and trend of groundwater storage from the annual low flows. This method is primarily based on the assumption of linear reservoir in groundwater recession and water balance equations. For a pristine catchment (i.e., human interference is negligible), if rainfall or meltwater can be ignored, the flow during the dry season is primarily the discharge from catchment groundwater storage and can be expressed as an exponential decay form as below following Bussinesq equation [18]:

$$y=y_0\exp \left(-t/K\right)$$

(1)

where $y=Q/A$ is the rate of low flow in the stream per unit of basin area, $Q$ is the volumetric rate of flow (L³/T), and A is the basin area (L²). The purpose of replacement of Q by y is for convenience to compare with other components of hydrological cycle such as evaporation, precipitation, etc. y₀ is the outflow rate at the arbitrary given origin time (t = 0). K is called the characteristic time scale of the catchment drainage process (T) determined by the basin physical characteristics [18,60].

Furthermore, groundwater storage S, the volume of water over the area A of the catchment, is the low flow integration of time t under the assumption (i.e., absence of rainfall, evaporation, or meltwater):

$$S={\displaystyle {\int }_t^{+\infty }ydt}$$

(2)

Substitution of (1) in (2) integration could yield a linear relationship between groundwater storage and the rate of flow:

$$S=Ky$$

(3)

Note that the groundwater storage S here should be considered as the average thickness of water storage above the zero-flow level of the water table over the catchment. Based on Equation (3), the temporal trend of groundwater storage can be determined from the temporal trend in low flow, namely:

$$\frac{dS}{dt}=K\frac{dy}{dt}$$

(4)

To meet the assumption that recharge to groundwater storage can be neglected, the annual lowest seven-day flow (hereafter denoted as y_l7), which is a widely used streamflow index in drought studies [13], was chosen in this study to substitute y in Equation (4), then Equation (4) is rewritten as:

$$\frac{dS}{dt}=K\frac{d{y}_{l7}}{dt}.$$

(5)

Essentially the y_l7 was chosen for two reasons. One is that low flow depends primarily on the depletion of groundwater, which satisfies the assumption of Equation (5) as much as possible [29]. The other reason for this is that lowest flow reflects the lowest catchment groundwater level in that year and thus trend in long-term groundwater storage–discharge dynamics can be tracked.

In addition to y_l7, the implementation of Equation (5) also requires knowledge of the characteristic time scale K. Although Equation (1) provides a solution of K, it is impractical to be used for determining the value of K because this expression is sensitive to the beginning of the flow that is primarily sourced from groundwater storage, which is nearly impossible to define a time origin for each recession in a long-term flow hydrographs [29]. To obtain a stable value of K, the method proposed by Brutsaert and Nieber [61] (hereafter denote as BN77) was used in this study, which can eliminate the variance of time t by derivation of time t in both sides of the Equation (1) as:

$$\frac{dy}{dt}=-\frac{y}{K}$$

(6)

Therefore, the value of K, i.e., the slope of −dy/dt versus y in a log-log plot, can be estimated from the lower envelope line with unit slope by excluding 5% points below the line. Lower envelope line is used here to minimize influence of rainfall, evaporation, and other influences on baseflow recession. And 5% is suggested by Brutsaert and Nieber [61]. Low flow data points are selected from the daily streamflow data through an objective algorithm proposed by Cheng, et al. [62]. In this automatic method, a series of stringent criteria are considered to guarantee the chosen flow data represent pure low flow, such as flow data with negative values of dy/dt, more than two points after the last positive and zero dy/dt, not the suddenly anomalous flow data that are much larger or smaller than the previous one etc.

2.4. Trend Analysis

Monotonic linear trends in hydrological and meteorological variables in this study were estimated using the rank-based non-parametric Mann–Kendall test [63,64]. In the Mann–Kendall test, the standard normal statistics T is given as followings:

$$T=\sum _{i=1}^{n-1}\sum _{j=i+1}^n\mathrm{sgn}\left(x_j-x_i\right)$$

(7)

$$\mathrm{sgn}(\theta )=\{\begin{array}{c}1\\ 0\\ -1\end{array}\begin{array}{c}\\ \\ \end{array}\begin{array}{c} if \\ if \\ if \end{array}\begin{array}{c}\theta >0\\ \theta =0\\ \theta <0\end{array}$$

(8)

where n is the length of data records, x_i and x_j are sequential data values. Then the normal statistic Z is derived by T and compared with the standard normal deviate $Z_{\alpha /2}$, and the test statistic Z is only statistically significant if $Z<-Z_{\alpha /2}$ or $Z>Z_{\alpha /2}$.

$$Z=\{\begin{array}{c}\left(T-1\right)/\sqrt{Var(T)}\\ 0\\ \left(T+1\right)/\sqrt{Var(T)}\end{array}\begin{array}{c}\\ \\ \end{array}\begin{array}{c}if \\ if \\ if \end{array}\begin{array}{c}S>0\\ S=0\\ S<0\end{array}$$

(9)

$$Var\left(T\right)=\left[n\left(n-1\right)\left(2n+5\right)-{\displaystyle \sum _{p=1}^q{m}_p\left(m_p-1\right)\left(2m_p+5\right)}\right]/18$$

(10)

where q is the amounts of tie groups and m_p is an arbitrary given p-th tie.

If a linear trend was present, the magnitude of the trend was estimated by the slope estimator $\beta$ proposed by Sen [65] and extended by Hirsch, et al. [66]:

$$\beta =Median\left[\frac{x_j-x_i}{j-i}\right]where 1<i<j<n$$

(11)

A positive value of $\beta$ represents an increasing trend, and a negative value of $\beta$ represents a decreasing trend [67].

3. Results

3.1. Long-Term Climate Variability of the YZRB

Figure 2 shows the linear trends of annual precipitation (Figure 2a) and temperature (Figure 2b) at the grid level during the study period of 1961–2000. It indicates that there are both increase and decrease in precipitation with rate ranging from (−18 to 44.5 mm/10a) (Figure 2a). Increase in precipitation mainly appears in upstream and downstream while the decrease mainly locates in the midstream. However, linear trends in almost all cells (72/73) are insignificant (p > 0.1). Trends of temperature are mostly positive (65 of 73) with rate ranging from −0.23 to 0.91 °C/10a (Figure 2b). The larger temperature increase shows in the midstream with trend over 0.4 °C/10a. Almost all the grids show significant (p < 0.05) positive trend in air temperature, indicating the overall climate warming in the YZRB.

Figure 3 shows the inter-annual variability of precipitation (P, Figure 3a) and air temperature (T_a, Figure 3c) and trends in monthly P (Figure 3b) and T_a (Figure 3d). Annual mean precipitation during 1961–2000 is 321 mm (ranging from 218 to 395 mm) and fluctuation of inter-annual P is generally moderate (coefficient of variation Cv = 0.13). There are several relatively dry periods during the late 1960s and early 1980s and 1990s. The long-term trend in annual P estimated in this study is about 0.2 mm/a which indicates that climate of the YZRB became slightly wetter, but the trend is statistically insignificant (p > 0.1). For the trends in P of different months (Figure 3b), 10 out of 12 months show an increasing trend, except June and November. Monthly precipitation trends range from 0.8 to 1.8 mm/a (about −2.4%~5.5% of monthly precipitation), but only the trend in April is statistically significant (p < 0.05). Trends in annual and monthly P are relatively small and mostly insignificant. Therefore, it can be regarded that there is no significant trend in annual P during the period of 1961–2000 and the seasonality of P has not changed.

In contrast, significant trends in T_a at annual time scales (Figure 3c) and in most of months (Figure 3d) were detected over the study period. The trend in annual T_a is 0.28 °C/10a (p < 0.05) estimated using the Mann–Kendall method, which is about two times larger than that in other regions at the global scale reported by IPCC AR5 (0.08 to 0.14 °C/10a during 1951–2012) [68]. Roughly T_a has increased about 1.1 °C over the study period. Increasing trends in different month range from 0.01 to 0.54 °C/10a, of which eight are statistically significant (p < 0.05). Apparent increasing trends in annual and monthly T_a indicate that climate in the YZRB during 1961–2000 became warmer at a much faster rate than the global average.

3.2. Changes in Annual Flow Regimes

The trends in mean and different percentiles (i.e., 10%, 25%, 50%, 75%, and 90%) of annual flows are were found in both 10% and 25% annual flows (Figure 4b,c). The significant trends in shown in Figure 4. Mean annual flow is 1810 m³/s and it shows a statistically insignificant (p > 0.1) downward trend (−6.2 m³/s/a, c.a. 0.34% of mean annual flow rate) (Figure 4a). Statistically insignificant (p > 0.1) negative trends are also found in 75% (−5.0 m³/s/a) and 90% (−27.8 m³/s/a) annual flows (Figure 4e,f). A slight positive trend (0.4 m³/s/a) is detected in 50% annual flows but is not statistically significant (p > 0.1) (Figure 4d). While statistically significant (p < 0.05) negative trends 10% and 25% annual flows are—1.6 m³/s/a and—1.7 m³/s/a, about 16.3% and 14.3% of corresponding mean annual flow percentiles over the study period, respectively. In addition, it is worth noting that the linear correlation between annual mean flow and 90% annual flows is very high (R² = 0.96) and their changes are very similar. It suggests that annual mean or total streamflow is dominated by the extremely large flow events and the distribution annual flows is highly skewed. Therefore, changes in mean flow may not be a good indicator for representing the changes in flow regimes in the YZRB. Overall, quantile and trend analysis of long-term streamflow records indicate that total runoff of the YZRB is relatively unchanged but annual low flows are decreased significantly during the period of 1961–2000.

3.3. Trends in Catchment Groundwater Storage

Annual lowest seven-day flow (i.e., y_l7) and its timing are shown in the Figure 5a,b, respectively. Annual y_l7 shows a significant (p < 0.05) downward trend (Figure 5a) with the rate of 1.09 × 10⁻³ mm/d/a (or 2.2 m³/s/a). Figure 5b shows that the date of annual y_l7 usually occurred during the period of February to April and its timing has advanced markedly by about 2.8 day/10a (p < 0.1). A significant decreasing trend and advanced occurrence date of annual y_l7 indicate that the lowest flow in the YZRB had changed continuously during the study period.

Based on the method of BN77, estimated characteristic time scale of groundwater recession (i.e., K in Equation (1) in study area is about 75 days (Figure 6). Inter-annual variability of groundwater storage of the YZRB is then derived from y_l7 using Equation (5) and is shown in Figure 7 (the lower orange line). Results show that annual groundwater storage had experienced two relatively dry periods (i.e., 1983–1985 and 1995–1997) and several relatively wet years (e.g., 1961, 1991, and 1999). Annual groundwater storage shows a statistically significant declining trend (p < 0.05) at a rate of 0.079 mm/a. Groundwater storage depleted by about 3.16 mm (c.a. 7.57 × 10⁸ m³) over the period of 1961–2000 in the YZRB.

4. Discussion

4.1. Changes in Low Flow Regime under Climate Change Conditions

In this study, estimated insignificant trends in annual precipitation (Figure 3a) and total runoff (Figure 4a) suggest that long-term precipitation–runoff relationship and water availability in YZRB is unchanged over the period of 1961–2000. This is consistent with previous studies by You, et al. [41] and Yang, et al. [69], who reported that annual precipitation showed an insignificantly increasing trend and annual runoff decreased slightly but not significantly. However, significantly and remarkable (~15%) decrease in low flow regimes are observed (Figure 4b,c and Figure 5a) and the occurrence date of the annual lowest flow is also found to become progressively earlier (Figure 5b) over the study period, which suggests that flow regime, spring melting and hydrological cycle in the YZRB possibly have changed already. Similar changes in low flow regimes and timing (i.e., decline in low flow and advancing timing of annual lowest flow) found in this study were also reported in other cold mountain regions, such as the Pacific northwest of the U.S. (PNUS) [9,16,57,70], where shifts in streamflow timing resulting from earlier spring melting were caused by climate warming [71]. Further, Kormos, et al. [15] attributed quantitatively rising temperature as the dominant factor to the decrease in low flow regimes from 42 stream gauges in the PNUS using path analysis. These studies in the PNUS shed light on the possible drivers (e.g., climate warming) for the changes of low flow regimes and the hydrological cycle in YZRB.

4.2. Changes in Low Flow and Groundwater Storage

Low flow regimes are closely related to dynamics in groundwater storage at the catchment scale [18,72], which provides the basis for estimating groundwater storage from low flow regimes. Based on the hydraulic groundwater theory, exponential decay functions (Equation (5)) has been widely used to described the relationship between groundwater storage and the lowest flow, which has been validated with observation records in the U.S. [27,28], Australia [29] and the cold eastern Siberia [33]. In practice, successful application of this relationship requires to check carefully whether the assumption satisfies or not (i.e., absence of precipitation, evaporation, leakage, and other water users). By introducing an additional water flux c into Equation (2), Cheng, et al. [60] systematically discussed influences of potential water users on low flow recession.

$$S={\displaystyle {\int }_t^{+\infty }\left(y+c\right)dt}$$

(12)

If there are non-negligible water users (i.e., c ≠ 0), a nonlinear lower envelope curve (i.e., −dy/dt versus y) would be in the log-log space. For example, in a paired catchment, a convex shape of low flow recession data grouped by year was found in treaded catchment after afforestation [60]. Similar recession analysis has been carried out to investigate groundwater evaporation [73], bedrock leakage [74] and human interferences [75]. In this study, low flow recession data grouped by year are shown in the Figure 6. No obvious change of recession data represents possible water users is negligible for low flow recession analysis and annual groundwater storage estimation.

Change in dry season (November to April) flow can be considered as an indicator of change in annual groundwater storage. Because, during the dry season, precipitation in YZRB is less than 10% of annual total precipitation and meltwater recharge is negligible as average temperature −7.7 °C. Figure 7 shows the inter-annual variability of dry season flow is very close to that of estimated groundwater storage (r = 0.88). Basically, the derived change in annual groundwater storage from daily streamflow observations is closely related to dry season flows, which reaffirms the decreasing annual groundwater storage during the study period of 1961–2000.

Limited by available streamflow observations, the derived evolution of annual groundwater storage based on runoff only covers the period before 2000. Fortunately, the GRACE can provide some helpful information for the groundwater storage changes after its launch on 2002 [20,23,24,25], although it is not suitable for long-term analysis with only a length of 13 years. Based on the GRACE record, outputs of the land surface model (e.g., soil moisture, snow water equivalent storage, etc.) and water mass balance theory, Feng, et al. [20] estimated that groundwater storage of the YZRB had decreased significantly (more than 10 mm/a), which is much bigger than that derived in this study. Cao, et al. [28] also found that the magnitude of the trends in groundwater storage estimated from the flow recession method is smaller than that derived from GRACE in the eastern United States (EUS). Although apparent differences in the magnitude of derived trends in groundwater storage existed, the trends derived using the two methods agree with each other in terms of the direction (i.e., increase or decrease) and groundwater storage estimated from flow recession method had quite similar patterns with well record over the most EUS [28]. Unfortunately, it is still difficult to build a certain relationship between low flow and GRACE data because, essentially, they represent different waters on and/or below land surface, although both have been widely used to infer changes in regional groundwater storage. Basically, reported downward trends in groundwater storage in the YZRB after 2000 based on GRACE data suggests that possibly groundwater storage has decreased continuously over the past six decades in the YZRB.

4.3. Impacts of Climate Warming on Groundwater Storage Decline

Climate warming has been reported as one of key drivers for groundwater storage trend in many cold mountainous catchments, such as alpine snow-dominated catchments [76], the north Canada [77] and the cold mountainous catchments of the western U.S. [78]. However, the effects of climate warming on groundwater storage are complex and are difficult to quantify as it can have both a negative influence (i.e., reduce groundwater) and a positive influence (i.e., increase groundwater storage). On one hand, warming can reduce groundwater recharge by lengthening growing season, increasing regional evapotranspiration, and decreasing snowfall. Climate warming had lengthened vegetation growing season over the North Hemisphere and played a dominant role in vegetation greening in the YZRB [79], which possibly lead to more soil water has been exhausted by plants transpiration and thus reduce groundwater recharge [80]. Warming is also expected to enhance regional evapotranspiration (ET) because a rising temperature can result in higher potential ET and thus lead to more ET from soil surface and plants, which can deplete more soil water that may recharge groundwater storage, which was indeed reported in the YZRB by Yang, et al. [9]. Given the unchanged annual precipitation, increasing ET can enhance groundwater storage depletion. Additionally, meteorological record shows a decrease in fraction of precipitation falling as snow in the YZRB [42], which is expected to decrease the natural high-mountain storage reservoir in the form of snowpack [81]. Reduction in the annual snowpack had been confirmed by monthly snow depth record (Figure 8) based on remote sensing data. It shows that snow depth in summer is much close zero and annual mean snow depth has a decrease trend (p < 0.1) with the slope of 2 mm/10a over the whole catchment, which means recharge from meltwater to groundwater storage has possibly reduced. One the other hand, degradation of frozen soil and glacier could increase groundwater storage. The degradation of widely distributed frozen soil is expected to increase groundwater storage by increasing active layer thickness [82,83]. Meanwhile, glacier retreat can play an important role in sustaining drought flow and in recharging groundwater storage [6].

Increasing intensity of precipitation events caused by climate warming during 1961–2005 reported by You, et al. [84] is possibly another reason that resulted in decrease of the basin-scale annual groundwater storage of the YZRB. This causal relationship between precipitation events intensity and climate warming has been well recognized globally as warming can increase atmospheric water content and lead to more extreme precipitation events [85,86]. Higher rainfall intensity is widely reported to result in lower infiltration rates and less cumulative infiltration, i.e., negative correlation between rainfall intensity and infiltration [87,88,89,90]. This negative relationship in rainfall intensity and groundwater storage was documented in a simulated experiment Wang, et al. [90] and with well observations in India [91]. Increasing rainfall intensity in the YZRB over 1961–2000 is thus expected to decrease the recharge of groundwater storage from precipitation. Overall, more sophisticated methods, such as process-based modelling or multi-source regression analyses over similar catchments, are required to quantify specific effects of climate warming on groundwater and flow regime in the future.

Changes in groundwater dynamics is of critical importance not only for understanding regional water cycles and managing local water resources effectively under changing environment [92], but also for understanding and predicting its ecohydrological consequences [93,94]. A decrease in catchment low flows in the YZRB means a wider range and/or variability of streamflow, which could translate into greater uncertainty for water management [16]. Previous studies on the land cover change indicated the expansion of the desert in the YZRB [95,96], which could be caused by the depletion of groundwater storage. Basin desertification is closely related to groundwater storage decline and has been used as a robust indicator for examining desert expansion [19]. Therefore, long-term changes of low flow regimes and groundwater storage found in this study could be of crucial importance for water resource management and environment protection in the YZRB.

In this study, the flow recession method from climatic perspective is demonstrated to be a reliable for understanding the changes in low flow regimes and groundwater storage in a cold mountainous catchment. This method has several advantages, such as requiring only daily observations of streamflow and providing an integrated value of groundwater storage over the entire catchment rather than points affected by various factors (e.g., geography, precipitation, or meltwater from scattered snowpack and glacier etc.). Thus, this method can be also applied successfully for understanding changes in water cycle in cold mountainous catchments, such as other catchments in the TP, where water availability is likely to be affected substantially by steady rising air temperature [97] and well observations are too limited to quantify groundwater storage changes accurately.

5. Conclusions

In this study, the impacts of global warming on cold mountainous catchment are investigated from a climatic perspective in a climate-sensitive catchment, namely the Yarlung Zangbo River basin. This climatic perspective is used to quantify climate warming in the regional water cycle through the evaluation of the basin-scale groundwater storage trend inferred from the recession of daily streamflow at annual scale, because the tendency of low flow is an indirect yet reliable indicator of long-term climatic impact signal and is largely reliant on groundwater changes. Results show that the annual low flow has decreased significantly and the date of annual lowest seven-day flow also occurs earlier during the period of 1961–2000. The annual minimum seven-day flow has decreased at a rate of 2.2 m³/s/a and its timing has advanced significantly about 3 day/10a. Inter-annual groundwater storage has declined significantly at a rate of 0.079 mm/a over the study period. These changes in low flow regimes and groundwater storage are accompanied by insignificant changes in precipitation, total runoff, and high flow, which indicates changes in water cycle in the YZRB over the study period. Further analysis suggests that evaporation increase, decreased snow depth, and increased annual precipitation intensity induced by climate warming in the study are considered as drivers of the detected decline in low flow regimes and catchment groundwater storage. This study provides a new perspective of changes in low flow and groundwater storage to detect and quantify the response of the regional water cycle under global warming. It highlights that global warming has likely already caused significant changes in regional flow regimes in the high and cold mountainous regions, such as the YZRB, along the south boundary of the Tibetan Plateau.