A time series assessment of terrestrial water storage and its relationship with hydro-meteorological factors in Gilgit-Baltistan region using GRACE observation and GLDAS-Noah model

Mountains regions like Gilgit-Baltistan (GB) province of Pakistan are solely dependent on seasonal snow and glacier melt. In Indus basin which forms in GB, there is a need to manage water in a sustainable way for the livelihood and economic activities of the downstream population. It is important to monitor water resources that include glaciers, snow-covered area, lakes, etc., besides traditional hydrological (point-based measurements by using the gauging station) and remote sensing-based studies (traditional satellite-based observations provide terrestrial water storage (TWS) change within few centimeters from the earth’s surface); the TWS anomalies (TWSA) for the GB region are not investigated. In this study, the TWSA in GB region is considered for the period of 13 years (from January 2003 to December 2016). Gravity Recovery and Climate Experiment (GRACE) level 2 monthly data from three processing centers, namely Centre for Space Research (CSR), German Research Center for Geosciences (GFZ), and Jet Propulsion Laboratory (JPL), System Global Land Data Assimilation System (GLDAS)-driven Noah model, and in situ precipitation data from weather stations, were used for the study investigation. GRACE can help to forecast the possible trends of increasing or decreasing TWS with high accuracy as compared to the past studies, which do not use satellite gravity data. Our results indicate that TWS shows a decreasing trend estimated by GRACE (CSR, GFZ, and JPL) and GLDAS-Noah model, but the trend is not significant statistically. The annual amplitude of GLDAS-Noah is greater than GRACE signal. Mean monthly analysis of TWSA indicates that TWS reaches its maximum in April, while it reaches its minimum in October. Furthermore, Spearman’s rank correlation is determined between GRACE estimated TWS with precipitation, soil moisture (SM) and snow water equivalent (SWE). We also assess the factors, SM and SWE which are the most efficient parameters producing GRACE TWS signal in the study area. In future, our results with the support of more in situ data can be helpful for conservation of natural resources and to manage flood hazards, droughts, and water distribution for the mountain regions.


Introduction
The Gilgit-Baltistan (GB), Pakistan, is located in High-Mountain Asia (HMA), which is the highest glacierized territory outside the Arctic and Antarctic, the so-called Third Pole, covering an area of over 100,000 km 2 and containing over 40,000 km 2 of ice bodies (snow glaciers) [21]. GB holds Indus River's catchment, a key source of water for agriculture and hydroelectricity production in Pakistan. GB also famous for its natural forests, mineral reserves, and some of the world's highest mountains makes it one of the most attractive and visited tourist destinations of local and foreign tourist in the country [36]. These large glaciers of GB contribute 70% to the Indus water flow [1,3,33]. Any small variation in temperature can be disastrous for the whole region and downstream population in context of climate change as it is projected that 1 °C increase in temperature would result in 15-16% increase in river runoff [44]. As snow and glaciers are very sensitive to climatic conditions, the rate of snow and glaciers melts will ultimately affect livelihood and safety situation of large population. There are contrasting research conclusions regarding glaciers of GB as authors [6,43] have concluded that snow cover is increasing in the area, whereas [18,22] conclusions are in contrast. Considering these inconclusive results about water resources (snow and glaciers), which will make difficult the management of water resources and its related events like flood or drought, there is an urgent need to use independent and newer technologies to assess the status of water storage in the region [5].
To study surface and sub-surface water resources through remotely sensed data, a new satellite of Gravity Recovery and Climate Experiment (GRACE) observation by NASA have provided data since 2002 for observing variations in terrestrial water storage (TWS) [2,14,30,31], ground water changes [9], and flood hazard assessments [50] at continental or regional level with few hundred kilometers resolution and uniform data coverage. The data provided from satellite gravity measures are complement to the field ground stations' data. However, to validate and increase the accuracy, the results needs to be evaluated based on parallel data sets provided by other remotely sensed sources [8,35,42]. The Global Land Data Assimilation System (GLDAS) [37], launched by NASA and NOAA, provides four core surface models which include VIC [29], Noah [26], Mosaic [27], and CLM [11] to simulate TWS. GRACE is a complement and can overcome the limitations of the other remotely sensed products and gauge station measurements with the main advantage that it provides water stored at all levels, including surface and sub-surface water, whereas the GLDAS-Noah model does not simulate the ground water storage [13]. However, according to the study conduct by [23], the ground water storage is stable (no apparent trend) in Upper Indus basin (UIB) as there is no any human intervention such as ground water pluming has been reported in the region. Therefore, ground water storage estimation has been neglected in this study.
In GB, Pakistan, except traditional hydrological studies [4,6,7,15], TWS is not considered at a regional scale with the GRACE observation data and technologies. In GB, the previous studies using traditional remote sensing technologies about water resources, i.e., glaciers and snow have contrasting results, few authors have concluded expansion and increase in glaciers and snow cover areas, whereas other contradicted their results [6,15,18]. In this research work, our main objective is to assess the TWSA in GB and to identify the factors that may play major role in variation. For this purpose, the water storage change is considered for GB region of Pakistan at a seasonal and monthly basis for the duration from January 2003 to December 2016. The monthly GRACE grids data (Level 2-RL05) from CSR, GFZ, and JPL processing centers were used. Further, to examine the influencing factors of TWS variability and to validate the outcomes, monthly grids of GLDAS-Noah extracted data with a 1° × 1° resolution are also investigated.

Study area
GB is located 72°-76° E and 34°-37° N. The total area is 68917.95 km 2 . The landscape of GB is unique in terms of glaciers and high mountains as it contains one of the highest mountain ranges of the world, i.e., Karakuram range and most of the area is glacierized (i.e., one of the largest glaciers in the world is located here, i.e., Baltoro Glacier). The geography of the GB is also unique as the Karakoram to the east, the Hindu Kush to the west, the Pamir to the north and the Himalaya to the south. The GB land area is about 87% situated 3000 m above sea level shown in Fig. 1. Moreover, it has many highest peaks of the world including world's second highest mountain K-2 (Mt. Godwin Austin) and Gasherbrum-I. The total glacier area of GB is about 13,082.94 km 2 . These glaciers and snow cover area have major contribution to the runoff of the Indus River. The climatic conditions of GB have large variation ranging semi-arid cold desert in the northern Karakoram and moist temperature zone in the western Himalayas where temperature varies from less than zero degrees to 40 °C based on season and location.

Data from GRACE
In this study, we used GRACE spherical harmonic coefficient data of level 2 released by three processing centers: the Center for Space Research (CSR) of the University of Texas at Austin, the German Research Center for Geosciences (GFZ) and the Jet Propulsion Laboratory (JPL) with maximum degree and order of 60 to computer terrestrial water storage (TWS) variation in the study area. Problem related to C20 coefficient would impact that the outcome of TWS estimated by GRACE is replaced with satellite laser ranging [10] and degree 1 coefficients are estimated based on method [40]. Glacial isostatic adjustment (GIA) has been done based on the model from [17]. A de-striping filter has been applied to the data to minimize the effect of correlated errors, and a Gaussian filter has been applied with smoothing radius of 300 km to smooth the data [12].

GLDAS model
GLDAS is a joint project of National Aeronautics and Space Administration (NASA)-Goddard Space Flight Center (GSFC) and the National Oceanic and Atmospheric Administration (NOAA). GLDAS is an international offline terrestrial modeling system that combines satellite and field station measurements to generate uniform land surface models or states (e.g., soil moisture (SM) and snow water equivalent (SWE) and changes) or predict in real time (e.g., evapotranspiration, rainfall) [37].
GLDAS drives four surface models: MOS, VIC, CLM, and Noah. In this research work, GLDAS monthly data of Noah model are taken from [32] with spatial resolution of 1 degree by 1 degree, for the period from 2003 to 2016. GLDAS-Noah dataset contains a time series of states related to land surface and changes, which we will use to investigate water storage. Besides, the variances related to the main portion of the signal to TWS will be observed which may be expected to increase from the variation in canopy water storage, SWE and SM. Consequently, the GLDAS-Noah model's data were used to partition the TWS changes into SM, SWE, and CWS components to understand how the variability of these hydrological fluxes plays a role in enhancing or dissipating the TWS. Additionally, groundwater water storage components are missing in GLDAS model, which can be calculated from the GRACE signal. It has been observed that integrated total water content estimated by GRACE is comparable with the GLDAS model's TWS [1,23]. Therefore, these land surface state variables were derived, and then, the TWS from the GLDAS model is calculated. In the first step, these time series of land surface and changes' states related to the GB are extracted, and in the next step, the TWS from GLDAS-Noah model is computed, as given in Eq. 1: In Eq. 1, ΔTWS GLDAS is the variation in TWS from GLDAS-Noah model, ΔSM is the variation in SM, ΔSWE is the variation in SWE, and ΔCWS is the variation in canopy water (1) ΔTWS GLDAS = ΔSM + ΔSWE + ΔCWS

In situ precipitation data
The daily precipitation data for the 12 weather stations for the period 2003 to 2016 by Pakistan Metrological Department (PMD), Government of Pakistan, were used for this study. Most of the weather stations are in elevation zone between 2000 and 4000 m.

Computation of trend, annual phase and amplitude
The long-term trends of the annual and semi-annual components of the above-defined data sets were estimated using least square regression as follows: where t represents time, 0 is the intercept, 1 is the slope, ( 2 , 3 ) represents the annual components, ( 4 , 5 ) are the semi-annual components, and ε is the random error.

TWS variations from time series
To determine the variability of TWS in Gilgit-Baltistan, TWS variations over the time drawn were examined from January 2003 to December 2016 using GRACE data and GLDAS-Noah model. Figure 2 shows residual (GRACE TWS-mean CSR, GFZ, and JPL) monthly TWS changes in Gilgit-Baltistan (mm) according to GRACE and GLDAS-Noah models. It can be observed that GRACE TWS time series from different data processing centers (CSR, GFZ, and JPL) and (2) GLDAS-Noah model are consistent. However, three data sets from GRACE (CSR, GFZ, JPL) observation and GLDAS-Noah model show a decreasing trend of (− 0.74 ± 0.31, − 0.37 ± 0.32, − 0.26 ± 0.41 mm/year) and − 0.79 ± 0.46 mm/ year, respectively, although the trend is not statistically significant. The annual amplitude, phase and trend of GRACE data sets and GLDAS-Noah model are shown in Table 1.
The highest annual amplitude measured for the periodic components by GLDAS-Noah model is 70.57 mm, followed by GFZ, CSR and JPL of GRACE observation being 57.94 mm, 41.26 mm, and 30.32 mm, respectively. The amplitude of GRACE TWS signal is lower than GLDAS-Noah model because GRACE estimated gravity field is composed of a set of spherical harmonic coefficient complete to degree and order 120. Spatial averaging and smoothing is required to reduce the contribution of noisy short wave length component in GRACE gravity field data as a result, degradation of actual signal may be possible [41]. Moreover, the SWE component time series (Fig. 8) shows that the minimum phases during summer-fall do not go lower than around − 40, while the peak variations match those of GRACE. It may be caused by permanent SWE building up or not melting completely during summer, which could be caused by the snow model deficiency or a lack of glacier model in GLDAS. Additionally, the annual phase shift estimated by GRACE (CSR, GFZ, and JPL data processing centers) and GLDAS-Noah model is almost the same with degree values of 255.500, 255.640, 256.320 and 254.590, respectively. The phase distortion describes how well the timing of positive and negative parts of TWSA time series matches. By comparing phase spectra of these signals, we get the delay spread information or dispersion, the information about the phase distortion.

Monthly mass change over the study period
The changes of TWS are determined based on GRACE (CSR, GFZ, and JPL) and GLDAS-Noah models in GB region over the 13 years period. Figure 3 shows residual mean monthly TWS change in GB (mm) from GRACE data and GLDAS-Noah model from 2003 to 2016. The particular  The TWS variation is negative from August to January, with a minimum of − 50 mm from October to November. Similarly, GLDAS-Noah models also show the same behavior with maximum of mean TWS greater than 100 mm in April and reached minimum in October with a value of almost − 100 mm as shown in Fig. 3. The monthly mean TWS variation in the region is corresponding to the variations of precipitation and snow fall in the region. The snow fall reaches its maximum in spring and its minimum in summer [22].

Seasonal mass change
The spatial changes in TWS estimated by GRACE (average of CSR, GFZ, and JPL) and GLDAS-Noah model for GB for the period January 2003 to December 2016 are shown in Fig. 4. The inverse distance-weighted (IDW) method of interpolation is applied to generate the maps and classified into seven classes to show variations in TWS. The TWS is low (negative) in winter (Fig. 4a, e) and autumn (Fig. 4d, h) seasons estimated by both datasets, but their amplitude of variation is varied in different locations. On the other hand, the TWS is high (positive) in whole study area for the spring (Fig. 4b, f) and summer (Fig. 4c, g) season estimated by GRACE and GLDAS-Noah models. This can be explained by the recent studies which emphasized that winter snow falling got shifted to spring and now there is more snowing in spring as compared to winters [5,22]. The maximum decrease in TWS is observed in western side of the study area for winter and autumn season estimated by GRACE solution, while GLDAS-Noah model shows decrease in TWS in southeast side for both seasons. For the spring season, TWS is low in eastern side estimated by GRACE solution, while GLDAS-Noah model shows decrease in TWS in southeastern side for the same season. Summer season also shows decrease in TWS in western side of the study area for GRACE solution, while GLDAS-Noah model estimates minimum TWS in southeastern region. Overall, the TWS is very low in winter and autumn, while it is high in spring and summer season as shown in Fig. 4. The spatial variation of TWS estimated by GRACE is converted into time series data, and a nonparametric Mann-Kendall's test (MKT) is applied to compute seasonal trends in the dataset as shown in Fig. 5. The MKT is a powerful technique generally used to compute trends in environmental variables [45]. A statistically significant (p < 0.05) decreasing trend of TWS is observed for the winter (DJF) and autumn (SON) season with a slope of − 3.382 and − 2.963 mm/year, respectively. However, the spring (MAM) and summer (JJA) seasons show an increasing trend of TWS with a slope of 0.295 and 0.414 mm/year but statistically not significant (p > 0.05).

Impacts on TWS changes
In the following section, we made scrutiny of the factors (precipitation, SM, SWE) which play major role in variability of TWS in the region.

Influence of precipitation
To determine the influence of precipitation on TWS changes, we utilized precipitation (mm/month) data from Pakistan Metrological Department (PMD). The total monthly precipitation for the period of 2003 to 2016 is represented in green histogram, while the TWS estimated by GLDAS-Noah model is represented in pink dash line as shown in Fig. 6. The three data sets of GRACE (CSR, GFZ, and JPL) which estimated TWS are represented in red, blue and black color. We compared the precipitation data to the GLDAS-Noah and GRACE data shown in Fig. 6. The TWS estimated by GRACE (CSR, GFZ, and JPL) and GLDAS-Noah model is increased with the increasing precipitation and vice versa, which indicates that precipitation has a direct positive correlation with TWS.
To make certain, the correlation between GRACE TWS time series and precipitation is computed using Spearman's rho method. The three datasets of GRACE (CSR, GFZ, and JPL) show significant positive correlation with precipitation with minimum value of r = 0.410 as presented in Table 2. Similarly, GLDAS-Noah model also has a statistically significant positive correlation with precipitation (r = 0.478).

Influence of soil moisture (SM)
The soil moisture data provided by the GLDAS-Noah data files from the global grids are concentrated to the GB region scale. Soil moisture measurements are summed for  four-layer model of Noah and transferred to mm and averaged. At the end, residual mean monthly changes are computed. Figure 7 shows the results of soil moisture changes (mm) and GRACE TWS of GB from 2003 to 2016. Soil moisture measurements range between − 75 and 60 mm, whereas GRACE TWS data range from 125 to − 102 mm. It can be determined that for GB, one of the parameters producing the GRACE TWS signal is the soil moisture. Table 3 shows the relationship between soil moisture and GRACE TWS using the correlation coefficients. The value of minimum correlation coefficient is 0.433 between CSR and soil moisture, which shows positive significant correlation. To summarize, soil moisture is one of the very important parameters in terms of TWS variation in the study area.

Influence of snow water equivalent
Likewise soil moisture data, the snow water equivalent (SWE) data provided by the GLDAS-Noah data files from the global grids are concentrated to the GB region scale using the same approach. Residual mean monthly changes are computed. Figure 8 shows the assessment results of snow water equivalent changes (mm) and GRACE TWS of GB from 2003 to 2016. Snow water equivalent measurements range between 150 and − 98 mm, whereas GRACE TWS data range from 125 to − 102 mm. It can be noticed that SWE results are in complete agreement with GRACE TWS time series variations. It can be concluded that for GB, one of another most important parameters developing the GRACE TWS signal is the SWE. According to the [20], major area of Upper Indus basin is covered with glaciers and snow. Table 3 shows the relationship between snow water equivalent and GRACE TWS using the correlation coefficients.
The minimum correlation coefficient is 0.45 which shows positive significant correlation between SWE and GRACE TWS time series. To summarize, SWE is an important parameter in order to monitor water storage variations and any unusual event related to water like flood or drought can be predicted and managed in advance from these results. The flood or water shortage management from GRACE data can be useful not only for economy and people of GB region but also for the whole country and millions of populations living in downstream country.

Discussion
The TWS storage estimated by GRACE (CSR, GFZ, and JPL) and GLDAS-Noah model for the study indicates that the TWS in GB has a slightly decreasing trend, while the trend is insignificant. The regional studies show contrasting results for similar elevations and uniformed temperature rise as the total mass in the Karakuram is stable or at gain, whereas in eastern and central Himalaya, the glaciers mass is losing [19,46]. Similarly, several studies found a negative trend using measured gravity signals for the Karakuram region [24,25,38], whereas many have opposite opinion or even they argue about mass gain [7,16,47]. Our results are in consistent with former researchers; however, the negative trend in GB region is insignificant. Furthermore, the annual amplitude and phase calculated from GRACE (CSR, GFZ, and JPL) and GLDAS-Noah model indicate that annual phased is almost the same for all data sets, but the amplitude is different. The annual amplitude of GRACE (CSR, GFZ, and JPL) is lower than GLDAS-Noah model, because spatial averaging and smoothing filters are necessary to apply on GRACE data to minimize the effect of noisy short wave length components of the gravity field solutions [41]. However, the GLDAS-Noah refers to a terrestrial modeling system that integrates ground-based observations and satellite data products to estimate land surface states and fluxes by using advance land surface models (LSMs) and data assimilation techniques. Long-term monthly analysis of TWS over the GB region indicate that increased pattern of TWS starts from January and reaches its maximum in the month of April. This corresponds to the maximum precipitation in this time period; according to the study of source [49], major precipitation occurred in winter and spring in Upper Indus basin originated from westerly disturbance. The TWS decreases from April and reached it minimum in October with the increase in temperature. Furthermore, seasonal analysis of TWS indicates that TWS is low in winter and autumn seasons, while it is high in spring and summer. The reason might be the decreased snowfall at high elevation in winter and increased evaporation [34,39]. Statistical analysis result indicates that the decreasing trend of TWS in DJF and SON is statistically significant with decreasing rate of − 3.38 and − 2.96 mm/year, respectively.
The precipitation is the major contributor for the variation of TWS. The TWS in the study region is increasing with the increase in precipitation and vice versa and have a positive statistically significant correlation. Previous studies also found high correlation between precipitation and TWS for the whole region of High-Mountain Asia which is in agreement with our results [48]. Furthermore, SWE and CWS also have a significant relationship with TWS, but the SWE plays major role in variations of TWS. Other regional studies argued that snowfall is the main contributor toward variability of glacier mass balance or TWS, which is in consistent with our results [28].

Conclusions
According to this research work, TWS studied for GB from GRACE and GLDAS-Noah derived TWS time series from 2003 to 2016 are consistent, which shows slightly decreasing trend. However, the trend is statistically insignificant. Our results are overall in agreement with few recent studies regarding stability of snow cover area and glaciers in GB. The annual amplitude of GLDAS-Noah model is higher than GRACE due to degradation of GRACE signal by applying decorrelation and Gaussian filter to reduce noise. The mean monthly TWS indicates that TWS reaches its maximum in April, while its reaches it minimum in October. This monthly variation of TWS is also in agreement with previous studies related to change in snow cover and glaciers in GB except in previous studies case the increasing snow cover and glacier starts from September. Our results help to conclude accurate results about the increasing or decreasing TWS trend in GB and having contradicting conclusions in earlier studies which do not use satellite gravity data. Seasonal analysis results indicate that TWS is decreasing in whole region in winter and autumn, while TWS is positive in spring and summer season. The TWS is increasing with the increase in precipitation which indicate that precipitation has high positive correlation with TWS estimated by GRACE and GLDAS-Noah models. Soil moisture and snow water equivalent are the most efficient parameters which play significant role in producing GRACE TWS signal and have a significant correlation with TWS. Further, to monitor and predict floods, drought, or any possible environmental conditions in the future, we suggest having an up-to-date data of following gravity observation of "GRACE-FO" based on satellite gravity. Additionally, increasing number of in situ measurements devices in GB region which in combination with GRACE-FO measurements will increase accuracy of forecasting and trend analysis to manage water and hazard for downstream population and bio-diversity.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.