Characterization of the Recharge-Storage-Runo ﬀ Process of the Yangtze River Source Region under Climate Change

: Storage and runo ﬀ are the two fundamental surface hydrological variables of a catchment. Research studies have been focused on the storage-runo ﬀ (S-R) hysteretic relationship of a catchment and its explanation very recently, thanks to satellite gravimetry. However, a complete analysis of a hydrological process starting from recharge to runo ﬀ has not been investigated. The S-R hysteretic relationship of Yangtze River Source Region (YRSR) situated in the northeast Tibetan Plateau is also unexplored. This study aims to investigate the Recharge-Storage-Runo ﬀ relationship of this catchment using gravimetrically-derived terrestrial water storage (TWS), satellite-derived and gauged precipitation, land surface modeled and gauged evapotranspiration, and runo ﬀ data measured during 2003–2012. We found that the Pearson correlation coe ﬃ cient (PCC) of S-R relationship is 0.7070, in addition to the fact that the peak of runo ﬀ every year comes earlier than that of the storage. This ﬁnding enables us to further calculate equivalent runo ﬀ based on water balance equation using the above data, while comparing to measured runo ﬀ time series. The comparison of Global Land Data Assimilation System (GLDAS)-derived (gauge-derived) equivalent runo ﬀ against measured runo ﬀ reveals a PCC of 0.8992 (0.9402), respectively, indicating both storage and runo ﬀ are largely controlled by the recharge derived from precipitation and evapotranspiration. This implies the storage is not coupled with runo ﬀ prominently due to steep topography in YRSR unable to hold the water in the form of storage. Exceptional anomalous water storage time series in 2006 has also been investigated. We speculate that the low rainfall might partly be related to an El Niño Southern Oscillation event. The low rainfall and abrupt groundwater transfer are likely to be the causes of the anomaly in 2006.


Introduction
A comprehensive process of a simple catchment should include input (e.g., precipitation), storage (e.g., soil moisture, groundwater, lakes, and glaciers), and release (runoff) [1]. Traditional hydrology attaches great importance to runoff simulation at a catchment scale based on precipitation as a collection input, referring to rainfall-runoff response (e.g., [2]). In this process, land water storage serves as an intermediate buffer relating rainfall to runoff, in which pore media are substantial hydrogeologic features that govern the storage behavior [3]. Recently, storage has been designed as a metric for catchments comparison [4]. The storage and runoff (S-R) relationship has also been explored to in the past few decades mainly due to the precipitation changes. Tang et al. [39] supplemented ice and snow melting caused by rising temperature as an additional reasoning.
The aforementioned studies focused on the storage and runoff trend and its potential underlying causes. Nonetheless, the seasonal and inter-annual characteristics, and potential hysteretic properties of the S-R relationship in the YRSR are still unexplored. In addition, as for the recharge element of the catchment, it is widely accepted to be the cause of the storage or runoff change. Long et al. [40] studied the relationship between recharge and storage in different parts of the Yangtze River Basin. The recharge is also regarded as a factor when analysing the runoff change in YRSR [37,38,41]. Chen et al. [42] also used satellite-based recharge data to assess the runoff for the entire Yangtze River Basin. However, no published materials have connected recharge with S-R relationship, not to mention analyzing their inter-relationships as a whole system. This paper aims to characterize the Recharge-Storage-Runoff Process as a whole in the YRSR, with a particular focus on the seasonal and inter-annual characteristics, and potential hysteretic properties of the S-R relationship. Time lags are quantified among storage, runoff, and recharge. Subsequently, quantitative explanation of anomalous water storage of a particular year will be illustrated and linked to El Niño Southern Oscillation (ENSO) [43]. It is structured as follows: Section 2 introduces the geography of YRSR and datasets; Section 3 presents the methodology after a pre-analysis of YRSR's S-R relationship, followed by calculating results in Section 4 that the equivalent runoff derived from water balance equation shows better correlation with observed runoff than storage; Section 5 concludes the resulting Recharge-Storage-Runoff relationship in YRSR, and discusses the causes of 2006 anomalous water storage event; Section 6 summarizes the conclusions.

Study Area
The Yangtze River Source Region (YRSR), located in northeast Tibetan Plateau, that is widely regarded as "the Third Pole", covers the northwest part of Sichuan Province, south part of Qinghai Province and the boundary area of Qinghai and Xizang Provinces ( Figure 1). Bounded by Kunlun Mountains in the north and Tanggula Mountains in the south, the region contains a number of sub-streams in the upper Yangtze River Basin [44], making it a significant ecological zone. Any abrupt hydrological changes in the YRSR would exert a great influence on most parts of China situated in the downstream [38]. Due to global warming, research studies in YRSR have been gaining concern since the last decade. Along with rapid development of remote sensing and data assimilation techniques nowadays, water storage of this area has been inferred from various sources, including GRACE, Global Land Data Assimilation System (GLDAS), WaterGAP Global Hydrology Model (WGHM) and water budget estimates [32,34,40]. Long et al. [40] made an attempt to restore the signal losses in GRACE data using scale factors from other model data sources. Data and models have also been used to assess the discharge in the study region [30,37]. It is a common consensus that both water storage and runoff show an upward trend validated by various datasets and modeling outputs, due to the rising temperature and precipitation [30,32,34,[37][38][39]. For instance, the increase in temperature causes a significant melting of permafrost and ice cover, complicating the hydrological regime of this region [45]. Meanwhile, because of the complex terrain, the meteorological and hydrological elements of the catchment manifest themselves considerable spatial heterogeneities [41]. With regard to the groundwater storage, Xiang et al. [35] has estimated its trend rate (i.e., 13.23 ± 12.03 mm/year).

In-Situ Meteorological and Hydrological Data
Considering that the evapotranspiration data from GLDAS-Noah model is the result of data assimilation rather than observation [46], in-situ evapotranspiration data from China National Meteorology Information Center (CNMIC), which were available at http://data.cma.cn/site/index.html, are used. Among the CNMIC datasets, Tuotuo river, Qumalai, and Yushu stations, which are located along the upper Yangtze River Basin, were chosen ( Figure 1).
The hydrological station at Gangtuo near the border of Sichuan Province is the only exit of the main stream discharge in the study area. Monthly runoff data during 2003-2012 was obtained by request from the Yangtze (Changjiang) Water Resources Commission, Ministry of Water Resources. Therefore, despite longer data time span of the below remotely-sensed and modeled data, the data analyses were restricted from 2003 to 2012. Additionally, Gangtuo daily discharge dataset (in m 3 /s) is converted into the runoff unit (in mm) by accumulating the daily data into monthly one divided by the area of study region (i.e., 149,600 km 2 ).

GRACE Data
Launched by NASA and German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt, abbreviated as DLR), Gravity Recovery and Climate Experiment (GRACE) satellite observes time-variable gravity field enabling the calculation of terrestrial water storage on a global scale [47]. As GRACE monthly data products are offered by different organizations around the world with various versions, the University of Texas Center of Space Research (CSR) Level-2 Release 05 (RL05), NASA Jet Propulsion Laboratory (JPL) RL05, GeoForschungsZentrum (GFZ) RL05, and the Institute of Theoretical Geodesy and Satellite Geodesy (ITSG2016) are chosen. Same processing steps were applied as follows: (1) degree-1 term generated from Satellite Laser Ranging (SLR) measured geocenter time series [48] was added, and C 20 term was replaced in spherical harmonic coefficients, respectively [49]; (2) a de-striping operation and a Gaussian filter with 350-km radius were applied to reduce the spatial-correlated noise at higher degrees [50][51][52]. The processed data products were then computed into equivalent water height (EWH) with a spatial resolution of 1 • [53]. By calculating the mean time series among the chosen data products, large anomalies for each individual data product become less apparent (Figure 2), which was employed for our investigation ( Figure 2). These data can be downloaded at ISDC (http://isdc.gfz-potsdam.de/grace).
To account for the missing monthly data during 2003-2012, equivalent water height, EWH(t), varying with time t (in months) with an initial epoch t 0 (set at January 2003), is represented in the form of a sinusoidal function accounting for yearly (i.e., c 1 and c 2 ) and half-yearly (i.e., c 3 and c 4 ) periodic parameters with an offset (i.e., a) and a trend (i.e., b) expressed as: After determining for the above six parameters, the missing monthly data can be interpolated via Equation (1).

TRMM Data
The Tropical Rainfall Measuring Mission, jointly set up by NASA and the Japanese Aerospace Exploration Agency (JAXA) in 1997, aims at monitoring the rainfall variation in the tropical and the subtropical regions. In order to retrieval the data from the TRMM satellite, algorithms have been developed including the TRMM Multi-satellite Precipitation Analysis (TMPA), which combined ground gauging data with satellite observations while outputting a data time series at a monthly interval covering the latitude 50 • N-50 • S [54], with a spatial resolution of 0.25 • [55,56]. We chose the monthly TRMM 3B43 Version 7 products, re-sampling its spatial resolution into 1 • × 1 • grid scale. The dataset can be found from the Goddard Earth Sciences Data and Information Services Center (GES DISC), available at https://disc.gsfc.nasa.gov/.

GLDAS-Noah Model
Global Land Data Assimilation System (GLDAS), widely used in hydrological study fields [46], was built by NASA with the goal to ingest satellite-and ground-based observational datainto the LSMvia data assimilation techniques [22,46]. There are four land surface models driven by GLDAS: Noah, Mosaic, the Community Land Model (CLM), and the Variable Infiltration Capacity (VIC) [57]. With a near global coverage at 0.25 • or 1 • spatial resolution, set of meteorological, hydrological and physical parameters are available at 3h, daily, and monthly temporal scales. Here 1 • × 1 • monthly datasets are used. Among them, the derived evapotranspiration, snow depth water equivalent, and soil moisture (with four soil layers: 0-10, 10-40, 40-100 and 100-200cm) from the GLDAS-Noah land surface model was employed for our investigation.

Storage-Runoff Relationship Analysis
To investigate the storage-runoff relationship, a phase diagram of runoff versus storage was generated. Without anomalous events, both runoff and storage time series should be periodic in nature.
The phase diagram should manifest as an oblique line when runoff and storage increase or decrease simultaneously, or an ellipse when time lag between runoff and storage is present. The presence of time lag is called the S-R hysteretic process.
In our study region, despite regular time series of runoff, TWSA time series in some years presents distinct pattern (Figure 3a). Another aspect is that the runoff's peak of each year comes earlier than that of the storage, opposite to our perception of regular hydrological systems. This results in a complex S-R phase diagram with a low correlation coefficient between runoff and storage (Figure 3b), attributable to internal (i.e., uncoupled storages [5]) or external factors (e.g., precipitation, evapotranspiration, monsoon and climate variability, etc) with a long period that disturbs the system [16]. Therefore, the above causes are required to be dealt with in order to increase the consistency between storage and runoff.  To investigate whether long-periodic components exist, spectral analysis has been conducted for these two time series (Figure 3c). Both time series exhibit strong periodic signal at yearly and half-yearly periods, whereas relatively weaker but obvious signals arise from the low frequency range (e.g., 40 months per cycle) in the TWSA power spectrum. This indicates the long-periodic components play a substantial role in the hydrological process, in particular TWSA.

Method Based on the Water Balance Equation
The Water Balance Equation is a generic formula for describing a hydrological cycle within the catchment, which is formulated as: where P, ET, S, and R are the precipitation, evapotranspiration, storage, and runoff, respectively. In the equation, the difference between P and ET can be defined as the recharge N to the system [58]. From Figure 4a,b, ET derived from the GLDAS-Noah model matches the P time series well with almost zero time lag, with strong periodic signals at yearly and half-yearly periods. In contrast, ET from CNMIC has a broader range of periodic signals when compared to that from the GLDAS-Noah model. The difference should lie in the measurement and model assimilation method, in addition to spatial scale. Apparently, the recharge from CNMIC is larger than that of GLDAS-Noah one, revealing that GLDAS-Noah modeled ET might be overestimated (Figure 4c). Through Equation (2), a runoff time series can be calculated by recharge N and derivative of storage dS/dt; we call this the equivalent runoff, Reqv, in this study. The derivative here is replaced similarly by a difference quotient through Equation (3) where ∆t is set to one month. During the calculation, it is important that the time lag should be accounted for. Previous hydrological and hydraulic research studies solved the problems related to time lag by importing some parameters, including hydraulic time constant τ, relative time lag ω, followed by a least-squares fit [5,58]. Here, the time lag removal procedure is done by interpolating the monthly time series into the daily one via a cubic spline while assuming 30 days each month, and phase shift is made for each year's data in order to determine the time lag that maximizes the correlation. Same procedure is applied to estimate the time lag between equivalent runoff and in-situ runoff.

Results
Before calculating the equivalent runoff, Reqv, time lag analysis is necessary for understanding of the hydrological environment in our study region. Time lag between storage change and recharge for each year has been determined ( Table 1) by maximizing the correlation (i.e., Pearson Correlation Coefficient (PCC)). After removing the time lag, it is apparent that the correlation between storage and recharge (i.e., PCC > 0.836) has increased significantly almost every year, when compared to that of Figure 3b (i.e., 0.707). However, regardless of time lag or not, the PCC between the time-lagged storage change and recharge for year 2006 is exceptionally low. The reason will be discussed in the next section. Note that the time lag here means the time required for the recharge to be partly absorbed into the storage. Overall, the time lag varies between one and two months each year. The time series of Reqv are displayed in Figure 5.  (Table 2). We found that the Reqv leads behind the in-situ runoff ranging from one day to about one month. Because the Reqv derived from the water balance equation is a macroscopic description of the catchment, it can be interpreted as the total surface discharge which includes the runoff in the river, and thus, the time lag here might be because of the recharge which delays the change of total surface discharge.  Modifying the Reqv time series by removing the time lag and detrending, the phase diagram of R-Reqv, shown in Figure 6, indicates approximately linear relationship between R and Reqv for the combination of 10 years' figures, despite fluctuation around a slope line. The Reqv using CNMIC gauging data (i.e., PCC of 0.9402) performs better than that of Reqv using GLDAS-Noah model (i.e., PCC of 0.8992). Probable reasons are that the GLDAS-Noah model's sum data does not contain new CNMIC gauging stations covering the study region for data assimilation process.

Characteristics of Storage-Runoff Relationship in YRSR
From Tables 1 and 3, it can be found that both water storage changes and runoff time series lag behind the recharge series. As mentioned in Section 3.1, the fluctuation of runoff has phases leading ahead the storage. The above analysis indicates that the recharge derived from precipitation and evapotranspiration is an important factor to both runoff and storage in YRSR, leading the changes with different time lags presented each year. Distinguishing from other typical catchments in the tropics, the water storage of our study region is significantly uncoupled from runoff. In other words, it manifests as a bypass flow that the significant portion of the recharge directly discharges away from our study region in the form of the runoff, rather than first residing in the form of water storage before discharging in the form of runoff. This should be attributed to the steep topography of the YRSR, where the water directly flows downstream along the sloped surface due to gravity effect. Thus, by eliminating the storage part of recharge via water balance equation, the study results perform well. This mechanism also makes sense to the result that the time lags between storage and recharge range from 1 to 2 months, much longer than the lag between runoff and recharge. In essence, recent studies have been conducted on the runoff variation in YRSR. Nevertheless, different interpretations and explanations are present. For instance, though precipitation has played the main role in runoff variation in the past decades [38,59,60], using the same gauging station data in YRSR, Li et al. [61] argued that the suprapermafrost dominated the runoff process in YRSR from June 2016 to May 2018 through isotopic analysis.
In view of the above, the analysis of runoff during 2003-2012 should synthesize both viewpoints. Though the precipitation (or recharge) represents a significant contribution to the runoff, the melting contribution of permafrost cannot be underestimated. Given the background of global warming, YRSR has experienced permafrost degradations [45], as a consequence of temperature increases for the past few decades [34,35,41]. Gao et al. [62] and Li et al. [63] suggest the increase of the runoff in YRSR might be due to the lateral melting contribution of permafrost and glacier. Besides, the baseflow of YRSR acts as the drainage of groundwater aquifers that represents part of the indirect recharge to the runoff. Research study also indicates that the baseflow variation relates to the precipitation and temperature orderly [44]. In summary, the runoff in YRSR is not only affected by the precipitation (or recharge) but also permafrost and glacier change caused by rising temperature. Similar conclusions have been found in Zhu et al. [64] and Tang et al. [39]. The above discussions justify our further investigation into the cause of the TWS anomalous change in year 2006.

Analysis of TWS Anomalous Change in 2006
As mentioned above, the storage in YRSR is controlled by the recharge. However, the storage time series in year 2006 seemingly does not obey this rule under normal circumstances, showing large inconsistency with the recharge series. Therefore, explanation of the cause is required for the anomalous TWS change in 2006. TWSA time series are plotted year-by-year, among which a distinct temporal pattern for the entire data time series in year 2006 is observed (Figure 7). In general, the water storage has two peaks within a year: the biggest peak comes in August or September, and the smallest one usually occurs during spring (i.e., March to May). However, the autumn peak of year 2006 vanishes, which results in its anomalously low PCC (Tables 1 and 2  Our result in the previous subsection suggested that the contribution of the recharge is more dominant than that of water storage. From Figure 4c, the most significant recharge every year is the rainy season between June and September, in which the total recharge of the rainy season is calculated (Table 4). We found that both GLDAS-derived and CNMIC-derived total rainy season recharge reaches the minimum value in 2006, significantly lower than that of other years. Considering the precipitation comes earlier than the TWS biggest peak, the TWS biggest peak in autumn should be related to the previous rainy season, as manifested from Table 1. Therefore, the disappearance of the TWS biggest peak in 2006 should be partly caused by low rainy season recharge in the catchment. We speculate that this anomaly might be partly related to ENSO. According to the conclusions of Zhang et al. [43], an El Niño (La Niña) event affects the TWS in Yangtze River Basin in terms of precipitation with a phase lag of 7-8 months, and there is an El Niño event continuing till the second half of 2005 [52]. As mentioned in Section 2.4, GLDAS-Noah model offers initial insight into different components of TWS from the internal aspect. According to Sproles et al. [16], TWSA change consists of ground water (GW) change, snow water equivalent (SWE) change, soil moisture (SM) change and reservoir change. In case of YRSR, the reservoir change should be replaced by the sum of the changes of lake (LA) and glacier (GA), such that the equation should be rewritten as: So far, the groundwater change has been investigated via removing other components from TWS using LSM models [35,36,[65][66][67][68][69] or integrated with geographic information system [70]. This is because of the lack of the groundwater measurements. Xiang et al. [35] and Zou et al. [34] stated that the water storage of lakes and glaciers in YRSR have been decreasing with a stable trend based on ICESat data. Nonetheless, this interpretation fails to explain the anomaly in 2006. Therefore, we further investigate the SWE and SM components.
Based on GLDAS-Noah model, the SWE and SM in YRSR are displayed in Figure 8a. Compared to TWSA time series, SWE peaks in February or March every year, earlier than TWSA's spring smallest peak. As for the autumn when TWSA shows the big peak, the SWE series is at its minimum value. Being one of the TWSA components, SWE displays a substantial effect on the winter season serving for the TWS spring smallest peak, while seemingly irrelevant to the autumn one. Thus, SWE can provide little information about the anomalous TWS in 2006. In contrast, SM data time series yield that the peak matches the TWS autumn biggest peak very well every year. More importantly, a significant decrease in TWS can be found exactly in Autumn 2006.
In the GLDAS-Noah model, soil moisture is divided vertically into four layers based on depth: 0-10, 10-40, 40-100 and 100-200 (unit: cm). In Figure 8b, soil moisture in the four layers all show a significant decrease around 2006, with an earlier sharp decrease takes place in the two deeper layers. The earliest decrease started in 2005 autumn in the 100-200 cm layer and it took almost a year to transmit to the top layer, suggesting the slow moisture changing rate might be associated with groundwater in a larger region. Xiang et al. [35] showed a decreasing trend of the groundwater in YRSR in 2006, while adjacent regions including Nujiang-Lancangjiang Rivers Source Region, Upper Jinsha River Basin, Yellow River Source Basin, Qiadam Basin, and the central Qiangtang Nature Reserve yielded an increasing trend. This proposed part of groundwater from YRSR flows to adjacent regions through hydrogeological tunnels in deep aquifers [35]. Though using different models, Bibi et al. [71] also revealed a sharp increase of groundwater in Qiadam Basin in the second half of 2006. Besides the rising temperature, water seepage and other natural causes, Chao et al. [36] suggests that the groundwater change in this region is also influenced by human activities, such as dam construction. Overall, the origin of TWSA anomaly in 2006 can be arisen in two aspects. First, precipitation in 2006 was at a relatively low level that fails to recharge the TWS. On the other hand, due to both natural and human factors, groundwater in YRSR was transferred to adjacent regions which further led to the sharp decrease of soil moisture in 2006.

Conclusions
In this study, we analyze a recharge-storage-runoff relationship in YRSR which might be potentially suitable for catchments situated at high altitude, mountainous topography and frigid environments. Datasets from satellites (GRACE, TRMM), hydrology model (GLDAS-Noah) and in-situ stations are used to calculate the equivalent runoff via water balance equation, aiming at quantitative explanation of the storage-runoff (S-R) relationship in Yangtze River Source Region during 2003-2012. We found that equivalent runoff correlated well with in-situ runoff with Pearson correlation coefficient of 0.8992 (ET from GLDAS-Noah model) and 0.9402 (ET from in-situ data), while compared to that of S-R relationship with a correlation coefficient of 0.7070.
The successful usage of the method promotes the analysis of the relationship between water storage and runoff. Dominated by water recharge (i.e., precipitation minus evapotranspiration) simultaneously, the hysteretic process between storage and runoff in YRSR is fairly weak. In other words, the storage and runoff variations are nearly synchronized rather than in a sequence (i.e., storage before runoff), behaving as if a bypass flow. This distinct mechanism makes sense to the fact that the runoff peak comes earlier than that of storage every year during 2003-2012.
Further investigation has been conducted to examine the anomaly of YRSR TWSA in 2006. Through the storage data at different depth layers from the GLDA-Noah model, two peaks (small one in Spring and big one in Autumn) of each year's TWSA series are associated with snow water equivalent and soil moisture respectively. For the abrupt changes of the TWSA time series in 2006, the causes come from two aspects. One is the low precipitation rate in 2006 possibly due to El Niño event in 2005 [43], the other is the decrease in soil moisture might be due to groundwater seepages flow to adjacent regions through some faults or underground tunnels in deep aquifers. Nevertheless, the verification of the second reason should be subject to further study.