Study on the Variation of Terrestrial Water Storage and the Identification of Its Relationship with Hydrological Cycle Factors in the Tarim River Basin , China

The terrestrial water storage anomalies (TWSAs) in the Tarim River Basin (TRB) were investigated and the related factors of water variations in the mountain areas were analyzed based on Gravity Recovery and Climate Experiment (GRACE) data, in situ river discharge, and precipitation during the period of 2002–2015. The results showed that three obvious flood events in 2005, 2006, and 2010 resulted in significant water surplus, although TWSA decreased in the TRB during 2002–2015. However, while the significant water deficits in 2004, 2009, and 2011 were associated with obvious negative river discharge anomalies at the hydrological stations, the significant water deficits were not well consistent with the negative anomalies of precipitation.While the river discharge behaved with low correlations with TWSA, linear relationships between TWSA and climate indices were insignificant in the TRB from 2002 to 2015. The closest relationship was found between TWSA and Pacific Decadal Oscillation (PDO), with correlations of −0.56 and 0.58 during January 2010–December 2015 and during January 2006–December 2009, respectively. Meanwhile, the correlation coefficient between TWSA and El Niño-SouthernOscillation (ENSO) index in the period of April 2002–December 2005 was −0.25, which reached the significant level (p < 0.05).


Introduction
The TRB is located in the heart of China's Silk Road economic belt and has abundant natural resources and extremely vulnerable ecological environment [1].The main water resource in the TRB is sourced from the surrounding mountains for the glacier snowmelt in alpine regions and precipitation in mid-mountains [2].The famous Silk Road with lots of fertile oases located at the edges of the TRB is nurtured by the glacier-snowmelt water.These oases play important roles and are the inward and outward kernel for economic and social development in the region [3,4].Because of the arid and semiarid climate, water is the main limiting factor of the development of oasis in the TRB [1].
Therefore, climate variations and human activities would impact on the terrestrial water resources which are important to human society [5].Thus, it is important to detect variations of terrestrial water storage (TWS), as well as the spatial and temporal variation of its processes [6].Meanwhile, the atmospheric circulation in large scale is connected with the hydrological cycle in the TRB.For instance, Wang et al. [2] found that the drought evolution in the TRB was affected by the northern hemisphere polar vortex, North Atlantic Oscillation (NAO), and Arctic Oscillation (AO); Li et al. [7] reported that the increasing precipitation in the northwest China, including TRB, owned to the South China Sea Subtropical High (SCSSH), West Pacific Subtropical High (WPSH), and North America Subtropical High (NASH); Liu et al. [8] concluded that the El Niño-Southern Oscillation (ENSO) index could be another factor for the variations of the dryness/wetness conditions in the northwest of China.Worldwide, the hydrological cycle can be identified by precipitation [9], river discharge [10], and other hydrological variables [11].Generally, while the drought indices are valuable tools for dryness/wetness conditions assessment, the Standardized Precipitation Index (SPI) [12], the Palmer Drought Severity Index (PDSI) [13], the stream flow drought index (SDI) [14], and the surface water supply index (SWSI) [15] are commonly and widely used to qualify the dryness/wetness level.There are multiple drought indices and methods to assess the dryness/wetness level in the TRB [16,17].However, there are not enough observation data in highly temporal and spatial variability for dryness/wetness conditions analysis.Additionally, the local precipitation and nonclimatic factors commonly impact on the river discharge variations [13,18].Because the different dryness/wetness processes make the hydrological variables have different characteristics, the work of dryness/wetness identification is incomprehensive [19,20] .
GRACE has provided an essential and valuable tool for analyzing variation in terrestrial water storage (TWS) [21,22], enriching the data for detecting dryness/wetness conditions at large scale.Many studies have detected the TWSA in many regions of China by using GRACE and other hydrometeorological datasets [23][24][25][26][27] such as the Tian Shan Mountain [28], the Yangtze River Basin [18], and North China [29].It can be noted that the usefulness of GRACE data can be prolonged by combining with other datasets.Moreover, the detection of dryness/wetness conditions can be solved through GRACE data [30][31][32].
To examine the dryness/wetness conditions and their potential factors in the TRB, the GRACE and Global Land Data Assimilation System (GLDAS) data were applied, combined with the meteorological and hydrological data.Furthermore, the river discharge and precipitation data in the subbasins were investigated, which might directly reflect the changes of TWSA in the TRB during the period of 2002-2015.In particular, the derived factors of TWSA during 2008 and 2009 were investigated.Additionally, potential teleconnections between TWSA and NAO, AO, ENSO, and Pacific Decadal Oscillation (PDO) were also assessed in the TRB.

Study Area
The largest and longest inland river in China, named Tarim River, is located in the TRB, with an area of 1.029 × 10 6 km 2 , ranging from 73 ∘ E to 97 ∘ E and from 34 ∘ N to 45 ∘ N, including high mountains and low plain regions (Figure 1) [33].The 114 rivers in the TRB connected with the mainstream reduced throughout history to just 3 rivers in the present day [33].The mean annual runoff of 3.989 × 10 10 m 3 is exhibited in the TRB, where the water component from ice and snow and precipitation in the mountains accounts for 48.2% [33,34].The three largest headstreams, Aksu River, Hotan River, and Yarkand River, account for 73.2%, 23.2%, and 3.6% of the total runoff, respectively [35].For the annual precipitation, while it is more than 300 mm in the mountainous region, it differs from 60 to 200 mm for the plain regions.Additionally, the annual air temperature in the TRB exhibits intra-annual variations [33].The barren and sparsely vegetated land plays major role in the land use in the TRB, which accounts for 56.2% [34].

Materials and Methods
The GRACE and GLDAS data were applied to detect the spatiotemporal distribution of TWSA.In addition, the precipitation, river discharge, and climate indices were used to investigate the potential factors for the changes of TWSA.  2.

Climate Indices Data.
Because the precipitation and river discharge were included in the hydrological cycle and lots of researchers found that the precipitation [7], stream flow and its extremes [2], and dryness/wetness conditions [36] were associated with the climate indices, the North Atlantic Oscillation (NAO), the Arctic Oscillation (AO), ENSO/NINO3.4 (i.e., the sea surface temperature locates in 170 ∘ W-120 ∘ W and 5 ∘ S-5 ∘ N), and the PDO during the period of 2002-2015 were selected to analyze the potential factors of the variation of TWSA from Earth System Research Laboratory (ESRL), http://www.esrl.noaa.gov.

GRACE Data.
The couple satellites in GRACE survey the spatiotemporal variations of earth's gravity field based on detecting the distance between the satellites.The variations of earth's gravity field could be converted to TWSA in the form of equivalent water height (EWH) after removing the influence of atmospheric/oceanic circulations and glacial isostatic adjustment (GIA).The methods for optimizing the accuracy and resolution of this conversion are still a subject of active research [25].
The GRACE data can be obtained beginning at April 2002 from http://grace.jpl.nasa.gov/.Based on the Level-2 RL05 data of GRACE in forms of spherical harmonic coefficients, the TWSAs were computed.Firstly, the atmospheric mass variations were removed from the gravity filed.Secondly, the C 2,0 coefficients were replaced with Satellite Laser Ranging (SLR).Thirdly, the degree-1 was estimated by the empirical formula.Then, the corrected GIA, destriping filter, and Gaussian filter with 300 km smooth radius were included in the data processing.Lastly, the monthly averaged TWS during January 2004-December 2009 was removed from the TWS time series during April 2002-December 2015.Although there were leakage error and measure error in the GRACE data processing, they were not taken into consideration in this article.The missing data in 2003-2015 were directly remedied by the linear interpolation.

GLDAS Data
. GLDAS aims at producing the precise condition of the land surface [37].There are four models in the system (e.g., Mosaic, Community Land Model (CLM), and Variable Infiltration Capacity (VIC)), http://mirador.gsfc.nasa.gov/.The soil moisture, canopy water, and snow water equivalent which are intimately connected with TWSA from these four models were recombined for comparison with the TWSA by averaging these four models.

Terrestrial Water Storage Calculations.
Reference [38] proposed that the gravity spherical harmonic coefficients can be processed into TWS based on the following equation: where  is terrestrial water storage in form of EWH and , 0, ,  ave ,   ,   , and   are the latitudes, longitudes, equatorial radius, mean density of Earth, love number, and coefficients of the spherical harmonics (Stokes' coefficients), respectively.Additionally,   (sin()) is the thdegree and mth-order fully normalized Legendre function, with maximum degree  and order m, expanded to 60.
Based on the method, the spherical harmonic coefficients of GRACE were converted into the terrestrial water storage in form of EWH.

The Gauss Filter Weight Function.
To verify the TWSA of GRACE, the monthly average value during the period of January 2004-December 2009 of equivalent TWS of GLDAS was removed from monthly equivalent TWS to obtain the GLDAS-based TWSA.Firstly, the null values area of GLDAS in the south polar was reassigned zero and then expanded to harmonic coefficients at 60 degrees.Subsequently, the harmonic coefficients were smoothed with Gaussian filter at 300 km smooth radius.Lastly, the processed harmonic coefficients were recombined to TWSA in form of EWH.The Gauss filter weight function could be expressed as follows: where  is equal to ln 2/(1 − cos( 1/2 /)), in which  1/2 and  are the Gauss filter smooth radius and the radius of Earth.The combined TWSA of GLDAS were filtered by the Gauss filter for verifying the TWSA of GRACE.

The Fitness Function.
Due to the significant seasonal variation and linear trend in the time series of TWSA in form of EWH, precipitation and river discharge would impact on the real signal; the seasonal variation and linear trend were removed from the time series.The seasonal changes could be obtained from the least-square method to fit the time series, as shown in the following equation: where  1 and  3 are the annual and semiannual amplitudes, respectively,  2 and  4 are the phase in the annual and semiannual,  5 indicates linear trend, and  6 is the constant. 0 is the first year of reference, that is, 1 January 2002.

Spatial Distribution of Terrestrial Water Storage.
Figure 2 shows the averaged spatial distribution of TWSA from three different organizations (e.g., Center for Space Research at University of Texas, Austin (CSR), Jet Propulsion Laboratory (JPL), and GeoForschungsZentrum, Potsdam (GFZ)) in the TRB from April 2002 to December 2015.Significant negative TWSA was detected in the north of TRB, while the positive TWSA was monitored in the south of TRB, especially in the southeast of TRB.In addition, the least TWSA was less than −15 mm, while the most TWSA was more than 10 mm.Pondering on the difference of the dataset from the three organizations, it was noted that there were no obvious differences among the datasets in the spatial distributions of TWSA.For the regional different variations of TWS in the TRB, glacial retreat and material accumulation were defined as the vital causes [1].Because of the increase of summer 0 ∘ C level height [39], accelerated glaciers retreat and large-scale land use expansion emerged in the north of TRB [40], but the decrease of glacial melt and material accumulation happened in the south of TRB for the decrease of summer 0 ∘ C level height [39,41].

Temporal Distribution of Terrestrial Water Storage.
Meanwhile, the distribution of TWSA from GRACE and GLDAS in the TRB from April 2002 to December 2015 was shown in Figure 3.It could be seen that the decreasing trend was monitored in all of them (i.e., CSR, GFZ, JPL, and GLDAS), and the amplitude and phase were in line with each other.The correlation coefficients of TWSAs between GRACE and GLDAS were 0.42 (i.e., CSR and GLDAS), 0.27 (i.e., JPL and GLDAS), and 0.42 (i.e., GFZ and GLDAS), respectively, as shown in Figure 3(b).The trends, at the rates of −1.6 ± 1.1 mm/a, −1.6 ± 1.1 mm/a, and −1.6 ± 1.2 mm/a, respectively, agree well among CSR, GFZ, and JPL, which were less than the rate of −4.2 ± 1.7 mm/a for GLDAS.This indicated that the other water resource factors During these three special periods, while the annual precipitations in 2005 (e.g., 135.7 mm) and 2010 (e.g., 163 mm) were 20.6% and 44.9% higher than the mean value (e.g., 87 mm) during the period of 2002-2015, respectively, the annual precipitation in 2009 was 19.2% lower than the mean value.Moreover, the second highest mean annual temperature of 9.7 ∘ C occurred in 2009, as shown in Figure 3(d).Thus, the significant low TWSA in the last month of 2008 and the early month of 2009 originated from the increasing evapotranspiration and subsequently caused an extra water loss from soil and surface water bodies in the TRB.Also, this proved that there was a severe drought event in the last month of 2008 and the early month of 2009.In addition, the previous study has found that the summer snow coverage rate reached a minimum amount of 2.46% in 2008 for the period of 2002-2012 based on the MODIS snow-ice product.Subsequently, the height of summer 0 ∘ C level reached the lowest value and the snow coverage rate was at less level in 2009 [39,42,43].As a result, a drought event happened in 2009 [42].

Relationship between TWSA and River Discharge Anomalies.
To research the significant drought and flood events, the seasonal and linear trend signals were deleted from the original time series of TWSA, river discharge, and precipitation in the TRB during the period of 2002-2015.Then, the nonseasonal time series were smoothed at thirteenmonth moving windows and normalized during the study period.The variations of TWSA in the TRB exhibited in Figures 4-6 are the mean value of TWSA from CSR, GFZ, and JPL.The river discharge of the whole basin was the average of the measurements at six hydrological gauging stations (i.e., Shaliguilanke, Xidaqiao, Alaer, Tongziguluoke, Kaqun, and Dashankou) in the TRB.
Figure 4 shows the normalized variations of TWSA, river discharge, and precipitation with nonseasonal signal and linear trend signal in the TRB.While the water depletions were detected during the period of 2002-2004, the drought event was reported in 2009.Although the significant negative river discharge anomalies in 2003 and 2004 were in line with the water depletion, there were no significant negative precipitation anomalies in the related period.Before the 2008 and 2009 droughts, there was a long-term positive TWSA in 2005-2007.Additionally, severe water depletions and drying up in 2009 resulted from the significant precipitation deficits and low height of 0 ∘ C in summer in the TRB [39].Subsequently, while TWSA and river discharge anomalies increased to be positive in a short period in 2010 and 2011, the abundant precipitation and river discharge happened.Additionally, there was another period with negative TWSA during 2013-2015 but without significant variations in precipitation.It was proven that there were complex relationships, which might have consisted of the human activities [40], the height of 0 ∘ C in summer [39], and atmospheric circulation [2,7,36] in water resource system in the TRB during 2002-2015.Also, the uncertainties in TWSA of GRACE were evaluated, which resulted from its relative coarse spatial resolution in the early period [18,44].In addition, the negative TWSA in 2003 and 2004 can be attributed to the negative river discharge and less snow cover rate in the same periods [42].More potential factors in hydrological variables were needed to verify the TWSA of GRACE.
As the inland basin region, TRB is surrounded by mountains, which are the source of the water for the glacier and snow ice melting.Although the correlation coefficients were not high, indicating that the impact of precipitation on TWSA was insignificant, the extreme events were associated with the factors among the TWSA, precipitation, and river discharge during 2002-2015.As shown in Figure 1, six specific hydrological stations, grids of GRACE in the mountain areas, and six weather stations were selected.The river discharge anomalies were normalized for further comparison.Then the correlation relationship between TWSA and river discharge anomalies and precipitation anomalies in the TRB at six hydrological stations was detected (Figure 5).The results showed that the river discharge at the six hydrological stations behaved with low correlations with TWSA.This indicated that the TWSA may be influenced by other factors.However, while the significant water deficits observed in 2004, 2009, and 2011 were associated with negative river discharge anomalies in all stations, the significant water deficits were not consistent with the anomalies of precipitation in the weather stations.
The influences of subbasins in the TRB for the 2009 and 2015 drought events were different.In 2009, the extremely low TWSAs over half of year, which were reported at all stations, and the extremely low TWSAs at Xidaqiao, Alaer, Tongziguluoke, Kaqun, and Dashankou stations were between the two extremely negative river discharge anomalies.This demonstrated that the river discharge from the mountains decreased in 2008 followed by the negative TWSA in 2009 under the impact of river discharge.Subsequently, the river discharge from the mountains increased in the last months of 2009 followed by the positive TWSA in 2010 under the impact of river discharge and precipitation.For the four stations (e.g., Shaliguilanke, Xidaqiao, Alaer, and Dashankou), both the negative TWSA and positive river discharge anomalies were observed in 2003.It can be deduced that the decrease of snow ice and glacier increased the river discharges downstream.Moreover, it decreased the TWSA in the mountain areas.Additionally, the positive TWSA and positive river discharge were found at Shaliguilanke, Alaer, Tongziguluoke, and Kaqun stations during 2005-2007.It implied that although the accumulation of TWSA happened less in the mountain areas for the snow and ice melt, the river discharge increased at the stations.
Despite the significant negative TWSA in the last months of 2008 and the early months of 2009, there were no significant positive anomalies of river discharge displayed at the hydrological stations.But both significant positive TWSA and obvious positive river discharge anomalies were detected at the hydrological stations in 2010 and the early months of 2011.As shown in Figure 5, this might be owed to the extreme precipitation in 2010 and the early months of 2011 in the TRB.In the last months of 2011, the river discharge anomalies and TWSA dropped down together.Also, there were, respectively, lowest TWSAs during 2014-2015.Additionally, there were obvious differences in the variation patterns of precipitation anomalies, river discharge anomalies, and TWSAs at different stations, especially at the hydrological station of Shaliguilanke.There were large high amplitudes in the nonseasonal precipitation anomalies at Shaliguilanke station, especially in the year of 2010.

Relationship between TWSA and Climate Indices. The relationships between TWSA and climate indices in the TRB
were also investigated in this study (Figures 6 and 7).In addition, based on the previous studies, the climate indices of AO, NAO, PDO, and NINO3.4/ENSO were selected to analyze the potential factors about variations of TWSA.Similarly, the seasonal signal and linear trend were removed from the hydrological factors, and the nonseasonal time series were smoothed at thirteen-month moving windows to better understand the relationships (Figure 6).As the results reported above, both the dry events around 2009 and 2015 and the wet periods in 2005 and 2010 in the TRB are exhibited in Figure 6 The significant surplus of TWSA and precipitation anomalies during 2010 owed to other factors, such as ENSO, and it was given that this water surplus was consentient with a negative ENSO phase (La Niña) (Figure 6(d)).The cool sea water in the eastern Pacific equatorial region would increase the precipitation on the western Pacific continent, thus increasing TWS.On the contrary, the warm condition of sea water in eastern Pacific equatorial region would decrease the precipitation on the western Pacific continent with decreasing the TWS [45,46].Similar cases were also investigated in the midstream of Yellow River and Qinghai province in China, which found that NINO3.4/ENSOimpacted on the regional precipitation and water conditions in 2008-2009 and 2010 [8,47].
Figure 7 reports the monopolized liner relationship between TWSA and climate indices.It can be noted that the linear relationships were not significant in the TRB during 2002-2015.In addition, the correlation coefficients between TWSA and climate indices in the TRB at different periods are detected (Table 3).Since the correlations of −0.56 and 0.58 between TWSA and PDO were, respectively, obtained during January 2010-December 2015 and January 2006-December 2009, the TWSA in the TRB was more closely related to PDO than NINO3.4/ENSO,AO, and NAO, as opposed to insignificant correlation between TWSA and NINO3.4/ENSO,AO, and NAO for the different periods (Table 3).Additionally, the correlation coefficient between TWSA and NINO3.extremes in many regions [48].Both the transport of water vapor and intensity of evaporation would be affected by the powerful atmospheric circulations; as a result, the conditions of TWS would be affected.

Conclusion
In the present study, the spatiotemporal variability of TWSA in the TRB during April 2002-December 2015 and the potential related factors of water sources in the mountain area were investigated using GRACE data, river discharge, and precipitation data.TWSA exhibited the decreasing trend in the TRB between April 2002 and December 2015, which was largely caused by the glacial retreat and melted snow ice.A few of drought events in 2002, 2004, 2009, and 2015 were observed from GRACE-derived TWSA, while the flood event in 2010 was detected from GRACE as well.Thus, it is also implied that while the quantification of hydrological extremes (e.g., floods and droughts) can be detected by GRACE, the combination between GRACE and other hydrological parameters would broaden the applications of GRACE in closed area.Comparing the six hydrological stations, it was found that the impacts from TWSA in the mountain areas were different at different stations.TWSA in the TRB was influenced by POD and NINO3.4/ENSO for the subperiods.The correlations of −0.56 and 0.58 between TWSA and PDO were obtained during January 2010-December 2015 and during January 2006-December 2009, respectively.Thus, the TWSA in the TRB was more closely related to PDO than NINO3.4/ENSO,AO, and NAO, as opposed to insignificant correlation between TWSA and NINO3.4/ENSO,AO, and NAO for the different periods.Additionally, the correlation coefficient between TWSA and NINO3.4/ENSOduring April 2002-December 2005 was −0.25, which reached the significant level.This research can provide useful information to detect the wetness/dryness conditions and their potential factors in the TRB.

Figure 1 :
Figure 1: The map of the study area.

Figure 3 :
Figure 3: The distribution of water resource factors.(a) The distribution of averaged TWS; (b) correlation relationship between GRACE and GLDAS; (c) time series of monthly TWS of GRACE and GLDAS; and (d) time series of monthly precipitation and temperature."EWH" stands for equivalent water height.

Figure 4 :
Figure 4: Time series of nonseasonal variations in TWSA, mean discharge, and precipitation in the TRB.The error bars indicate the uncertainties of GRACE data.
. However, there was no monopolized relationship between TWSA and climate indices for the different periods in the TRB during 2002-2015.As for AO, while the positive correlation relationship with TWSA was exhibited in the TRB during 2009-2012, the negative correlation relationship was displayed during 2002-2009 and 2012-2015 (Figure 6(a)).Pondering on NAO, while the positive correlation relationship with TWSA was exhibited in the TRB during 2009-2010, the negative correlation relationship was displayed during 2002-2009 and 2010-2015 (Figure 6(b)).Taking PDO into consideration, while the positive correlation relationship with TWSA was exhibited in the TRB during 2005-2010, the negative correlation relationship was displayed during 2002-2005 and 2010-2015 (Figure 6(c)).As for NINO3.4/ENSO,while the positive correlation relationship with TWSA was exhibited in the TRB during 2006-2008, the negative correlation relationship was displayed during 2002-2006 and 2009-2015 (Figure 6(d)).
4/ENSO in April 2002-December 2005 was −0.25, which reached the significant level.While the opposite correlations before and after December 2005 were found between TWSA and AO, NAO, PDO, and NINO3.4/ENSO, the converse correlations before and after December 2009 were discovered between TWSA and AO, NAO, PDO, and NINO3.4/ENSO.After December 2009, the positive TWSA was consistent with two continued La Niña events in 2011-2012 and the longterm significant cool PDO phases imposed various climate

Table 1 :
The summaries of weather stations in the TRB in this study.

Table 2 :
The summaries of hydrological stations in the TRB in this study.

Table 3 :
Correlation coefficients between TWS/precipitation and teleconnection indices.Note: the numbers in bold represent the correlations at 95% level of confidence.