Increasing annual streamflow and groundwater storage in response to climate warming in the Yangtze River source region

The area of the Yangtze River source region is about 13.77 × 10⁴ km², which is located in the central QTP (figure 1). It is one of the most representative alpine areas in China with the most concentrated biodiversity and it has undergone significant changes due to climate warming (Jiang et al 2015, Grosse et al 2016, Wang et al 2017). The main rivers in the region are the Tuotuo River, Dangqu River and Chumar River. The downstream of the Tuotuo River is referred to as the Tongtian River, which is the largest trunk stream and its change could affect the streamflow in the whole source region (Chorography Compile Committee of Qinghai Province 2000). The climate in the source region is a continental and transitional type, which is mainly controlled by the upper Westerly belt (Cao 2013). The mean annual air temperature ranges from −5.5 °C to −1.7 °C, the annual mean precipitation is 375.4 mm, and the annual evaporation is between 1340 mm and 1700 mm (Li et al 2020).

Zoom In Zoom Out Reset image size

Streamflow data were obtained from the most upstream hydrological station of the Yangtze River (Zhimenda Station). Temperature, precipitation, and evapotranspiration were obtained from the four meteorological stations in Tuotuohe, Wandaoliang, Qumalai, and Yushu (table 1, http://data.cma.cn). Weighted average of the meteorological data were calculated by the Thiessen polygon (Thiessen 1911, Bhavani 2013, Panno and Luman 2018). Permafrost data were simulated by the top temperature of permafrost (TTOP) model (appendix) (Smith and Riseborough 1996, Nan 2003) based on the MODIS LST dataset (http://ladsweb.nascom.nasa.gov/data) offered by the National Aeronautics and Space Administration. The five kinds of data (streamflow, temperature, precipitation, evapotranspiration, and permafrost) are available from 1962 to 2012. The latest terrestrial water storage data estimated by reconstructing the GRACE satellite data (resolution of 0.5° × 0.5°) during 1990–2012 (Li et al 2021) were used to compare with the groundwater storage estimated by our baseflow recession analysis.

2.3.1. Trend analysis and mutation detection

2.3.2. Wavelet analysis

2.3.3. Baseflow separation

The recursive digital filtering method is one of the most reliable methods for separating baseflow due to its high accuracy, since it eliminates the random interference and proves the detection accuracy (Bozic 1976, Eckhardt 2005, Stamenković et al 2019):

$$\begin{equation}{\text{B}}{{\text{F}}_t} = \frac{{\left( {1 - {\text{BF}}{{\text{I}}_{\max }}} \right)\omega {\text{BF}}{{\text{I}}_{\max }}{Q_t}}}{{1 - \omega {\text{BF}}{{\text{I}}_{\max }}}}\end{equation} \tag{ 1 }$$

where ${\text{B}}{{\text{F}}_t}$ is the daily baseflow (m³ d⁻¹) subject to ${\text{B}}{{\text{F}}_t}$ < ${Q_t}$, $t$ is the time step (day), ${Q_t}$ is the daily discharge (m³ d⁻¹), $\omega$ is a dimensionless recession constant and $\omega < 1$, ${\text{BF}}{{\text{I}}_{\max }}$ is the maximum value of the ratio of baseflow to total streamflow with ${\text{BF}}{{\text{I}}_{\max }} < 1$. The ${\text{BF}}{{\text{I}}_{\max }}$ is tightly related to the total baseflow, and $\omega$ controls the shape of baseflow hydrograph.

2.3.4. Estimation of groundwater storage

Groundwater storage is estimated by three main functions (Brutsaert and Nieber 1977, Malvicini et al 2005):

$$\begin{equation}S = K{y^\lambda },\end{equation} \tag{ 2 }$$

where $y$ is the volumetric flow rate in winter baseflow (m³ d⁻¹), $S$ is the total amount of groundwater storage in the catchment aquifers (m³), $K$ and $\lambda$ are the dimensionless constants depending on the hydrogeological conditions of the study area. $K$ is the baseflow recession coefficient, which is a parameter that could quantitatively represent the streamflow recession process of the catchment. $\lambda$ is a parameter that typically ranges from 0.5 to1.0.

During the dry season, when there is no precipitation and other inputs, the conservation of mass equation can be expressed as follows (Brutsaert and Nieber 1977, Parlange et al 2001, Brutsaert 2008, Lin et al 2020):

$$\begin{equation} - \frac{{{\text{d}}S}}{{{\text{d}}t}} = y.\end{equation} \tag{ 3 }$$

Substitution of equation (2) into equation (3) leads to (Brutsaert and Nieber 1977, Mendoza et al 2003, Rupp et al 2004, Rupp and Selker 2006):

$$\begin{equation}\frac{{{\text{d}}y}}{{{\text{d}}t}} = - a{y^b}\end{equation} \tag{ 4 }$$

where ${\text{d}}y/{\text{d}}t$ is the temporal change in the baseflow during recession, and constants $a$ and $b$ are the intercept and slope of the plots of $- {\text{d}}y/{\text{d}}t$ versus $y$ in log–log space. Parameters of $K$ and $\lambda$ in equation (2) can be expressed in terms of $a$ and $b$ as $K = 1{\textrm{ }}/\left[ {a\left( {2 - b} \right)} \right]$ and $\lambda = 2 - b$ (Lyon et al 2009, Gao et al 2017).

Figure 2 shows the 5% (Q_5%), 10% (Q_10%), 25% (Q_25%), 50% (Q_50%), 75% (Q_75%), and 90% (Q_90%) dates of the annual total streamflow between 1962 and 2012. For example, the 5 percentile date (Q_5%) is the day by which the 1st 5% of the annual streamflow has occurred in the catchment. The 10% (Q_10%), 25% (Q_25%), 50% (Q_5%), 75% (Q_75%), and 90% (Q_90%) dates are defined in the same way. The Q_50% date is similar to the 'date of center of mass of flow' (Stewart et al 2005, Hodgkins and Dudley 2006). According to these data, we can infer the relationship between streamflow and seasons. The occurrence time of the six typical annual streamflow (Q_5%, Q_10%, Q_25%, Q_50%, Q_75% and Q_90%) in the Yangtze River source region shows a significantly upward tendency except Q_25% and Q_50% from 1961 to 2012 (figure 2). The slopes of the dates of Q_5%, Q_10%, Q_75%, and Q_90% are increasing mean the dates of these four typical flows are delayed. Delaying dates indicates there are much more precipitation first becoming groundwater and then discharging into river instead of discharging into the river directly. The slopes of the dates of Q_25% and Q_50% are decreasing mean the dates of these two typical flows are advanced. The advance in dates is mainly caused by the melt of the solid water remaining in winter and the dates are advanced by rising temperature. The dates of Q_25% and Q_50% are in summer (from June to August) according to the results, so we estimate that the variations of streamflow are closely related to seasons.

Zoom In Zoom Out Reset image size

In order to determine the mechanisms of the streamflow variation, we first chose the recursive digital filtering method to separate baseflow from the streamflow. Figure 3 shows that the increasing trends of the streamflow and baseflow are close to each other, with the slopes of 1.68 m³ s⁻¹ and 1.35 m³ s⁻¹ per year, respectively, indicating that the increase of baseflow plays an important role in the increase of streamflow. ${\text{BFI}}$ in the graph ranges from 0.68 to 0.74, which means the annual mean ${\text{BFI}}$ is relatively stable.

Zoom In Zoom Out Reset image size

The Mann–Kendall test shows that annual streamflow and winter baseflow exhibit a significantly increasing trend at the 95% confidence level ($\alpha$ = 0.05) during 1962–2012 (figure 4). The values of ${\text{UF}}$ show the variation of corresponding time series, ${\text{UF}} < 0$ represents the decrease of time series, and the larger the absolute value of ${\text{UF}}$ is, the greater the decrease is. Figure 4 shows that there is a significant decrease of both streamflow and baseflow near the year 1966, and a shortly increase since 2010. The annual streamflow decreases more sharply in 1972–1981 and 1996–2000 ($\alpha$ = 0.05), while the winter baseflow is in 1979–1981 ($\alpha$ = 0.05). The intersection at the beginning is not representative, and it has been often ignored in the Mann–Kendall test (Li et al 2011). Therefore, annual streamflow's mutation point of 1965 is removed and the mutation years of the annual streamflow are 2008 and 2011 (figure 4). The mutation year of winter baseflow is 2008, which is the same as the result of previous studies (Li et al 2017).

Zoom In Zoom Out Reset image size

The contour map of the real part of wavelet coefficients depicts the periodic characteristics of streamflow time series in different time scales. When the real part of the wavelet coefficients is positive, it represents wet year, and negative represents dry year. The vertical coordinates corresponding to the centers of contour lines represents the main period scales and the horizontal coordinates represents the wet years or dry years. The characteristics of wetness–dryness alternation at main period scales are helpful to predict streamflow variations in the next period scales. Figure 5(a) shows that the 1st main period scale is 37 a for streamflow variation. There are three wet years (1966, 1989, and 2008) and two dry years (1976 and 1998) of the annual streamflow in the 37 a, which means it has five obvious wetness–dryness alternations (first increases and then decreases) at the period scale of 37 a.

Zoom In Zoom Out Reset image size

The wavelet variogram can reflect the distribution of fluctuation energy and determine the main period precisely in different time series, and it has been widely used in analyzing streamflow variations. There are five obvious peaks in figure 5(b), corresponding to the time scales of 11 a, 20 a, 26 a, 16 a, and 37 a from small to large. The maximum peak is 37 a, indicating that the period scale of 37 a is the strongest, and 37 a is the 1st main period of the streamflow variation. That is to say, there is a 37 years periodicity (five obvious wetness–dryness alternation at the period scale of 37 a) in streamflow variation, which is consistent to the results in figure 5(a). The periodic fluctuation energy of the remaining time scales of 16 a, 26 a, 20 a, and 11 a decreases successively. Time scales of 37 a, 16 a, 26 a, 20 a, and 11 a are the 1st to 5th main periods of annual streamflow. The above five periods control the characteristics of streamflow variations during 1962–2012.

Similarly, the winter baseflow wavelet variogram indicates that the winter baseflow also have five significant periods, from the 1st to the 5th main periods are 34 a, 15 a, 24 a, 18 a, and 10 a, respectively (figure 5(c)). The periodic oscillation energy in 34 a is the largest, controlling the variation in the winter baseflow (figure 5(d)). Continuous hydrologic years should be chosen during the wavelet analysis of winter baseflow. As there are no streamflow data in 1974–1976, 1983–1984, and 2010, winter baseflow cannot be calculated in 1962, 1977, 1985, and 2011. The winter baseflow's period scales are smaller than those of the annual streamflow. But the five main period scales of both winter baseflow (34 a, 15 a, 24 a, 18 a, and 10 a) and annual streamflow (37 a, 16 a, 26 a, 20 a, and 11 a) are similar. According to the obvious wetness–dryness alternations (first increase and then decrease) of streamflow and baseflow in their 1st main periods, both the streamflow and baseflow are possibly increasing significantly in the near future.

Daily baseflow records in winter (December–February) were used in calculating the groundwater storage. Figure 6(a) shows all data points of $- {\text{d}}y/{\text{d}}t$ versus $y$ in log–log space during 1962–2012. The fitted slope $b$ of the blue dashed line is equal to 1.02 via the nonlinear least squares fit from equation (4). It should be noted that the value of real intercept is its corresponding log value in figure 6(a), it is $\log _{10}^a = - 0.615$ according to equation (4). The values of $K$ and $\lambda$ then can be calculated as 42.25 and 0.98, respectively. The total amount of the groundwater storage $S$ for each year can be directly estimated by equation (2). Figure 6(b) compares the groundwater storage estimated by baseflow recession analysis during 1962–2012 and the terrestrial water storage estimated by reconstructing GRACE data (Li et al 2021) during 1990–2012. The height difference between two curves mainly includes the soil moisture storage and snow water equivalent storage. The two curves in figure 6(b) are well correlated with a Pearson's correlation coefficient (Pearson 1895, Fischer et al 2017) of 0.63 during 1990–2012, which is higher than the previous study (Liu et al 2020). The coefficient indicates a generally good correlation between the baseflow-derived groundwater storage and the GRACE-derived measurements. As groundwater storage is essentially directly derived from streamflow data, there are some uncertainties in the estimated groundwater storage. One of the main uncertainties is from neglecting glacier melt water. Human activities and vegetation were also not taken into account.

Zoom In Zoom Out Reset image size

Figure 7 depicts the variations in three kinds of parameters (slope $b$, intercept $a$, and coefficient $K$) in equation (4) of all the recession data points (from December to February), they are processed by moving average according to the 1st to 5th main periods of streamflow (sliding windows at 37 a, 16 a, 26 a, 20 a, and 11 a, respectively). Figure 7 shows that the variations in these three parameters (slope $b$, intercept $a$, and coefficient $K$) are related to the five main periods of the annual streamflow. The stronger the periodicity is, the larger the absolute values of slopes in three pictures are (figures 7(a)–(c)). This indicates there is a positive correlation between the increase of groundwater recession coefficient $K$ in equation (4) and the 1st to 5th main periods of the streamflow in the Yangtze River source region from 1962 to 2012. As shown in figure 7(a), the variations of slope $b$ of the fitted recession curves in equation (4) during each period gradually move upward, which are opposite to the results of the intercept's variation (figure 7(b)). As the slope $b$ increases, the intercept $a$ gradually decreases, and the recession coefficient $K$ increases (figure 7(c)). Furthermore, the three figures are at a confidence level of 95% (α = 0.05), and this long-term variation in the recession coefficient $K$ from December to February indicates that baseflow recession gradually speeds up in the catchment during winter.

Zoom In Zoom Out Reset image size

Analyzing variation of streamflow and its influencing factors in the Yangtze River source region is vital to the variation of water resources in the middle and lower reaches of the Yangtze River and the entire QTP. Streamflow and baseflow show decreasing trends during 1960–1979, while ${\text{BFI}}$ shows an increasing trend. The streamflow and baseflow reached a minimum in 1979 (figure 3). The streamflow and baseflow fluctuate significantly after 1980. The reason is that the annual mean precipitation in the region increased rapidly in the 1980s, 1990s, and 21st century (Wu et al 2020). So the period from 1980 to 2012 was selected as the period for focused analysis.

The streamflow in the source region is correlated with climatic factors including temperature, precipitation, and evapotranspiration. In order to investigate the periodic characteristics and multi-year variation of the climatic factors, daily data from meteorological stations of Tuotuo River, Wudaoliang, Qumalai, and Yushu between 1962 and 2012 (table 1) were analyzed according to the Thiessen polygon method. The wavelet periodic analysis and the Mann–Kendall test were also used to investigate the variation characteristics of the various climatic factors (table 3). The result of Morlet wavelet analysis shows the 1st to 6th periods of the annual mean temperature are 26 a, 36 a, 16 a, 8 a, 11 a, and 5 a, respectively (figures 8(a) and (b)). The 26 a period is the most significant period scale that can reflect the variation characteristics during 1962–2012. According to the obvious wetness–dryness alternations (first increase and then decrease) of temperature in figure 8(a), we predict that the annual mean temperature would be a rising tendency in the near future (Michael 2019). The periods of the annual mean precipitation are 37 a, 26 a, 20 a, 17 a, 7 a, and 10 a, respectively (figures 8(c) and (d)), which shows a similar increasing trend as temperature. For evapotranspiration, the periods are 36 a, 26 a, 21 a, 17 a, and 6 a, respectively (figures 8(e) and (f)). The wavelet variogram of evapotranspiration implies that it will also be increasing in the following years.

Zoom In Zoom Out Reset image size

On the basis of the Pearson correlation coefficient between the three climatic factors and streamflow, the correlation from high to low is annual average temperature (0.68) > rainfall (0.43) > evapotranspiration (0.34). Such correlations suggesting that the increase of streamflow in the Yangtze River source region during 1962–2012, especially after 1980, is closely related to climate warming.

Permafrost thawing under climate warming can significantly influence the mechanism of streamflow (Walvoord and Kurylyk 2016). A reasonable hypothesis is that increased groundwater storage in winter is highly related to permafrost thawing, and the permafrost thawing may enlarge the groundwater storage capacity (Lyon and Destouni 2010, Niu et al 2016). The TTOP model has been validated to be of high adaptability in simulating the permafrost areal extent and distribution in the Yangtze River source region (Smith and Riseborough 1996, Nan 2003, Reich et al 2017). The results show that the permafrost area in the catchment has undergone great changes. Figure 9 clearly shows the degradation of permafrost from 1962 to 2012. The distribution of seasonally frozen soil or no permafrost area in 1962 is mainly concentrated in the southeast of the study area, and then extended to the northwest area along the river net in 2012. The permafrost area decreased from 13.39 × 10⁴ km² (97% of the study area) in 1962 to 12.57 × 10⁴ km² (91% of the study area) in 2012, which has decreased by 8200 km² (figure 10). The minimum permafrost area of 9.69 × 10⁴ km² (70% of the study area) occurs in 2006. There is a significant downward trend after 1995. Figure 10 shows the annual mean permafrost area is decreasing at a rate of 401 km² a⁻¹ during 1962–2012, and total groundwater storage $S$ shows an increase trend at a rate of 1.62 km³ a⁻¹. Since the variations of the permafrost area and $S$ are highly relative to streamflow and streamflow has a significant downward trend after 1980, we choose the typical period of 1980–2012 for a detailed analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 shows that the annual mean permafrost area decreased at 577 km² a⁻¹ between 1980 and 2012. In order to eliminate the effect of extreme values, all the data in figure 11 were processed by moving average (a sliding window at 10 years) based on the data in 1980–2012 from figure 10. The variation of permafrost area during 1997–2012 is much more fluctuant than that of 1980–1997. As shown by the cumulative anomaly curves (figure 12), the anomaly of permafrost area is 1997, so we divided the trend analysis into two sections before and after 1997. The annual mean permafrost area decreasing rates are −334 km² a⁻¹ in 1980–1997 and −425 km² a⁻¹ in 1997–2012, respectively (figure 11). The fluctuation in 1997–2012 is comparatively larger than that in 1980–1997 and its slope is closer to the overall trend, revealing that the variation of permafrost area is mainly influenced by the period after 1997. The year of 1997 is close to the mutation point of annual mean temperature (1998), indicating that the variation of permafrost area is closely related to the variation of temperature. Meanwhile, the permafrost area shows an increasing trend and the groundwater storage shows a decreasing trend after 2008 (figure 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The groundwater storage increases from 279 km³ in 1980 to 309 km³ in 2012 and its increasing rate is 3.43 km³ a⁻¹ in this period (figure 11). The cumulative anomaly of groundwater storage is 2004 as shown in figure 12. The weighted average of the mutations of all the time series can be calculated out is 2000 (table 4). So we divided groundwater storage into two periods (before and after 2000) for a better understanding of its trends (figure 6(a)). The blue solid line in figure 6(a) depicts the fitted linear trend of $\log ( - {\text{d}}y/{\text{d}}t)$ in 1962–2012, with a slope of 1.02. The linear trends in 1962–1999 (the black dash line) and 2000–2012 (the red dash line) are also shown in figure 6(a). The change of intercept $a$ due to slope $b$ in equation (4) can be observed to be an increasing trend from 0.0241 in 1962–1999 to 0.0245 in 2000–2012. The result means the absolute values of both the coefficient $K$ and the groundwater storage are increasing about 1.8% after 2000.

The correlation coefficients between groundwater storage and permafrost area and between groundwater storage and temperature in 1980–2012 can be calculated to be 0.74 and 0.67 by the Pearson correlation coefficient method. The correlation coefficients indicate that they are negatively correlated, and the variation of groundwater storage is more closely correlated with the permafrost area than temperature. Therefore, melt water in the permafrost is released as permafrost thaws, influencing the variation of water resources in the whole Yangtze River source region. The groundwater storage in the near future is predicted to be increasing according to the periodicity analysis of the baseflow, as baseflow first increases and then decreases in its 1st main period (figure 5).