Characterization of Basin-Scale Dynamic Storage–Discharge Relationship Using Daily GRACE Based Storage Anomaly Data

Despite the fact that streamflow occurs mainly due to depletion of storage, our knowledge on how a drainage basin stores and releases water is very limited due to measurement limitations. A window of opportunity, however, is provided to us by GRACE (Gravity Recovery and Climate Experiment) satellite mission that provides storage anomaly (San) data. Many studies have explored a range of potential applications of San data such as flood forecasting. Here we argue that the capability of GRACE satellite mission has not been fully explored as most of the studies in the past have performed analysis using monthly San data for large river basins. In this study, we use daily San data for several mid-sized catchments to perform storage–discharge analysis. Our results support the earlier notion that storage–discharge relationship is highly dynamic. Furthermore, we demonstrate that San data can be exploited for prediction of k of the Brutsaet–Nieber equation −dQ/dt=kQα (Q is discharge at time t). For comparison we also use storage information provided by Catchment Land Surface Model (CLSM) as well as past discharge information to predict k. Our results suggest that GRACE based storage information can be used to predict k reasonably well in gauged as well as ungauged basins.


Introduction
Prediction of streamflow is crucial for planning and management of many water resources infrastructural projects. Unfortunately, it is not easy to predict streamflow as we do not know very well how rainfall produces streamflow. Hydrological processes responsible for streamflow generation are complex and our ability to observe them is quite limited. Thus, it is practically impossible to predict streamflow by employing fundamental equations of water flow that requires detailed information on the surface and subsurface characteristics of a basin. Hydrologists generally predict streamflow using simple conceptual models capable of performing reasonably when limited data is available. However, after the advent of satellite based remote sensing technologies our ability to obtain hydrologically relevant information has improved dramatically. Hydrological models too have evolved over the past few decades to properly utilize information obtained through satellite remote sensing. For example, many hydrological models are now well equipped to handle fine resolution topographic and land-cover information obtained through satellite remote sensing [1][2][3][4]. Similarly, soil moisture data products obtained through microwave remote sensing are widely used these days to improve streamflow prediction [5][6][7]. One limitation of these data products is that they provide moisture content information of the top soil layer only, whereas for streamflow prediction purposes it is important to obtain information on the hydrologically-active total water content in a basin at a given time [8][9][10][11][12][13]. In this regard the Gravity Recovery and Climate Experiment (GRACE) satellite mission launched on 17 March 2002, by the National Aeronautics and Space Administration (NASA) and the German Aerospace Centre (DLR), has made a remarkable achievement. GRACE mission has been providing total water storage anomaly information of the earth surface by monitoring the spatiotemporal variation of Earth's gravity field since 2002 [10,[14][15][16][17][18][19]. The fundamental measurement is derived from micron level tracking of the satellite-to-satellite distance, which varies due to individual gravitational attractions on the satellites as they pass over the Earth's surface. Further, mathematically, storage anomaly (S an ) for a basin at time t is given as: S an (t) = S(t) − S av (1) where S(t) is storage at time t, which includes water stored in all the surface and sub-surface storage units of the basin. S av is the average storage of the data set from period 2002 April to 2009 December. Since the main aim of a typical hydrological model is to explain how a basin stores precipitated water and releases it as streamflow or discharge, it is important to characterize the relationship between storage and discharge of the basin. GRACE based S an data is found to be useful by many for hydrological modelling studies such as detecting trends of anthropogenic groundwater depletion, hydrological flux estimation and climate model improvement, sea-level change and ocean dynamics, operational drought monitoring [20]. Further, Total Water Storage Anomaly (TWSA) data was used to predict flood in large river basins [10,21]. Further, Fang and Shen exploited GRACE based S an data for providing information on high flows as well as low flows in rivers [22]. GRACE based data was also used for measuring drainable water storage [23]. It should be emphasized here that the above studies have largely focused on use of monthly S an data for performing storage-discharge analysis for large river basins. A few studies have been carried out using daily GRACE data as compared to monthly data. A recent study used daily GRACE data to track major flood changes in the Ganga-Brahmaputra delta [24]. As hydrological complexity decreases with scale [25,26], it is expected that storage-discharge analysis at smaller spatial and temporal scales will reveal finer details about the hydrological processes of a river basin. The main aim of this study is to use daily S an data for providing information on basin-scale storage-discharge relationship. Daily S an data is the ITSG-Grace2016 gravity field model, computed in Graz University of Technology, providing unconstrained monthly and Kalman smoothed daily solutions. It has been calculated by assuming that the gravity field does not change arbitrarily from one-time step to the next. They have used geophysical models the WaterGAP global hydrology model (WGHM), the atmospheric model ECMWF, and the ocean circulation model OMCT to find the information about the temporal correlation patterns. Further utilizing this knowledge, the temporal resolution has been enhanced without losing spatial information within the framework of a Kalman smoother estimation procedure [27]. The stochastic prior information derived from geophysical models and the daily GRACE observations are included in the Kalman smoother to deliver an updated state of the gravity field for each day. The stochastic information is introduced in terms of the process model, which is constructed from spatial and temporal covariance matrices derived from the output of the geophysical models. Further, they used model output of the years 1976-2000 (i.e., outside the GRACE time span) to guarantee that the GRACE solutions are not biased towards the model. In particular, we attempt to utilize daily S an data to quantify the dynamic relationship between storage and discharge, characterized by the power-law coefficient of −dQ/dt vs. Q curve during streamflow recession. Estimation of the power law coefficient is important because it can be used to characterize the interaction of subsurface water and surface water systems [28] and to undertake prediction of recession flow which can be used for water resource management during no flow period [29]. It can be used for estimating drainable storage of the catchment [13]. Additionally, the coefficient can be used to obtain information on active drainage network [30] and baseflow index [31]. To assess the relative usefulness of S an data, we also consider storage information provided by GLDAS (Global Land Data Assimilation System) [32][33][34][35][36] and antecedent discharge, a proxy of past storage [12,37,38] in our analysis. Section 2 provides information about the data products used in this study, and Section 3 introduces the methods. In Section 4, the results and discussion about this study are presented. Finally concluding remarks are given in Section 5.

Study Catchments and Preliminary Data Analysis
A total of 51 basins with size ranging from 60 to 8500 km 2 were selected from the USGS database for this study (Table A1). The location of gauging stations is given in the map provided in Figure 1. In Figure 1 grey grid lines represent the GRACE data footprint used in the analysis, yellow circles correspond to location of gauging stations and both of them were plotted over a precipitation map of USA. The background of the map displays mean precipitation in mm/day from the period 1948 to 2010 (data collected from: https://psl.noaa.gov/data/gridded/data.unified.daily.conus.html) to provide an idea on the hydrological conditions of the study basins. GRACE based gridded S an dataset was obtained from TU Graz (Graz University of Technology), which is available at https://www.tugraz.at/institute/ifg/downloads/gravity-field-models/itsg-grace2016/ [39][40][41].
Geosciences 2020, 10, x FOR PEER REVIEW  3 of 16 recession flow which can be used for water resource management during no flow period [29]. It can be used for estimating drainable storage of the catchment [13]. Additionally, the coefficient can be used to obtain information on active drainage network [30] and baseflow index [31]. To assess the relative usefulness of S an data, we also consider storage information provided by GLDAS (Global Land Data Assimilation System) [32][33][34][35][36] and antecedent discharge, a proxy of past storage [12,37,38] in our analysis. Section 2 provides information about the data products used in this study, and Section 3 introduces the methods. In Section 4, the results and discussion about this study are presented. Finally concluding remarks are given in Section 5.

Study Catchments and Preliminary Data Analysis
A total of 51 basins with size ranging from 60 to 8500 km 2 were selected from the USGS database for this study (Table A1). The location of gauging stations is given in the map provided in Figure 1. In Figure 1 grey grid lines represent the GRACE data footprint used in the analysis, yellow circles correspond to location of gauging stations and both of them were plotted over a precipitation map of USA. The background of the map displays mean precipitation in mm/day from the period 1948 to 2010 (data collected from: https://psl.noaa.gov/data/gridded/data.unified.daily.conus.html) to provide an idea on the hydrological conditions of the study basins. GRACE based gridded S an dataset was obtained from TU Graz (Graz University of Technology), which is available at https://www.tugraz.at/institute/ifg/downloads/gravity-field-models/itsg-grace2016/ [39][40][41]. The dataset provides daily S an data from April 2002 to December 2016 at one degree by one degree spatial resolution. However, it should be noted that the actual spatial resolutions of the gridded daily TWSA data is lower with approximately 3° × 3° and 5° × 5° for monthly and daily solution, respectively. Although it is recommended to use GRACE derived TWSA data for river basins with drainage area greater than 200,000 km 2 [42], the data is useful even for small basins in USA. For example, Scanlon et al. 2016 used GRACE derived TWSA data for hydrological applications in 53 basins with drainage area ranging from 40,000 to 100,000 km 2 and found the data to be useful in smaller basins [43].
We obtained daily S an time series for each study basin from the gridded dataset. Further, first, a 30 m resolution digital elevation model (DEM) data was used to delineate drainage basins following D8 flow direction algorithm which was introduced by O'Catlaghan and Mark [44]. The D8 algorithm divides the plane into eight directions: (1) east, (2) north-east, (3) north, (4) north-west, (5) west, (6) south-west, (7) south, (8) south-east. The D8 flow directions are assigned to each grid following the concept of steepest gradient, i.e., water flows in the direction in which the gradient is the steepest and then for each outlet point (outlet points are gauging stations and their details are given in Table A1), we obtained boundary shape-files. The entire processing of getting boundary shapefiles was using The dataset provides daily S an data from April 2002 to December 2016 at one degree by one degree spatial resolution. However, it should be noted that the actual spatial resolutions of the gridded daily TWSA data is lower with approximately 3 • × 3 • and 5 • × 5 • for monthly and daily solution, respectively. Although it is recommended to use GRACE derived TWSA data for river basins with drainage area greater than 200,000 km 2 [42], the data is useful even for small basins in USA. For example, Scanlon et al. 2016 used GRACE derived TWSA data for hydrological applications in 53 basins with drainage area ranging from 40,000 to 100,000 km 2 and found the data to be useful in smaller basins [43].
We obtained daily S an time series for each study basin from the gridded dataset. Further, first, a 30 m resolution digital elevation model (DEM) data was used to delineate drainage basins following D8 flow direction algorithm which was introduced by O'Catlaghan and Mark [44]. The D8 algorithm divides the plane into eight directions: (1) east, (2) north-east, (3) north, (4) north-west, (5) west, (6) south-west, (7) south, (8) south-east. The D8 flow directions are assigned to each grid following the concept of steepest gradient, i.e., water flows in the direction in which the gradient is the steepest and then for each outlet point (outlet points are gauging stations and their details are given in Table A1), we obtained boundary shape-files. The entire processing of getting boundary shapefiles was using ArcGis 2010. For each basin spatial weighted average rule was followed to obtain S an each day. Weightage of a grid pixel for a basin is the area of the pixel falling within the basin boundary.
For comparison purposes, we also obtained storage data provided by GLDAS [32]. We particularly used storage data provided by Catchment Land Surface Model (CLSM) [45]. Unlike other land surface models that assume uniform topographic and hydrologic characteristics at the grid scale, CLSM divides the land surface into topographically defined catchments and models hydrologic processes based on each catchment's topographical statistics. It provides storage information at daily [45,46] timescale from 1948 to present. CLSM model based total water storage data (S cl ) is obtained for the study basins following the same spatial average technique discussed earlier.
Available daily discharge (Q) time series data for each study basin was obtained from the USGS website https://waterwatch.usgs.gov/. Since we were interested in analyzing recession flows, we first delineated the recession curves for each basin. A recession curve is defined as a continuously decreasing streamflow time series lasting at least 5 days [47]. We also identified the recession curves from each basin for which −dQ/dt continuously decreases with time to ensure minimization of the effect of errors in our analysis [48,49]. Note that we selected basins that are relatively free from human interventions using satellite maps (courtesy, Google Earth) to ensure that the streamflow time series were not considerably altered. The basins are located either in rural areas, forest area or sub-urban areas of USA states. Further we avoided the basins where dams have been constructed. The basins used in our analysis are basically from a list of basins used by Swagat et al. [38].

Theoretical Backgrounds
Storage (S) change in a basin at any time (t) has to be equal to inflow rate (I) minus outflow rate (O): dS/dt = I − O. Inflow is either rain or snowfall, whereas streamflow (Q), evapotranspiration (ET), and groundwater loss (GL) constitute outflow (O). GL corresponds to outflow from the basin through subsurface pathways. Although the mass balance equation provides a foundation to determine the relationship between S and Q for a basin, the main challenge in solving the equation is that some of the quantities are unobservable. In particular, direct measurement of storage S and GL is practically infeasible. Large-scale ET measurements are generally associated with significant errors, particularly when ground measurements are not available for calibration purposes [50]. An easy step towards simplification of the mass balance equation is to focus on recession periods when inflow I is zero, i.e., dS/dt = −O = −(Q + ET + GL). If all the three outflow variables (Q, ET and GL) decrease together in time during a recession event, we can expect to see a strong functional relationship between Q and dS/dt: Q = f(dS/dt). This assumption is particularly useful because dS/dt = dS an /dt, which means we can directly use GRACE based storage anomaly information for determining the storage-discharge relationship of a basin. Such a possibility, however, is mired with the fact that dS an /dt observations are expected to be associated with large errors with respect to change in discharge measurement per day (see Figure 2a). Large observational errors are likely involved because we are extracting S an for each basin from the GRACE dataset which is having a very poor resolution with respect to the average size of our study basins [51]. One alternative to the obstacle above is to express discharge as a function of storage itself: Q = f(S an ), since compared to dS an /dt, S an is not expected to be influenced a lot by observational errors. This is because as S an (t) = S(t) − S av , when we take dS an , i.e., S an (t) − S an (t − 1), we will be left with S(t) − S(t − 1), as S av will get cancelled. The term S(t) − S(t − 1) will have lot of noise and is not suitable for storage discharge analysis. (a) There is almost zero correlation between dSan/dt and Q due to high observational errors associated with dSan/dt, suggesting that dSan/dt cannot be used to predict Q. Although Q vs. San plot exhibits less scatter, the correlation is not very strong, likely because of the fact that the relationship is not one-to-one as shown in (b). Note that the figure was prepared using data from the sample basin of Bluestone river near Pipestem, West Virginia with USGS id 03179000 having drainage area of 395 square miles. Figure 2b shows the Q vs. S an plot for the sample basin. Although the correlation of this plot is better than that of the dS an /dt-Q plot, there is still a large amount of scatter, which is because the relationship between S an and Q can change from one recession event to another [33,48,49,[52][53][54][55]. Figure 3 essentially shows that although the relationship S an and Q can be strong for a recession event, the relationship may vary considerably across events. It essentially suggests that we cannot construct a single mathematical relationship between S an and Q for a basin to predict streamflow. We need to first characterize the dynamic relationship nature of storage-discharge relationship. One avenue in this regard is proper characterization of the relationship between dQ/dt and Q during recession periods following the method proposed by Brutsaert and Nieber [56]. Generally, the relationship between dQ/dt and Q can be expressed as −dQ/dt = kQ α , from which we can easily determine the storage-discharge relationship of the basin [48,57,58]. Biswal and Marani noticed that the value of α generally remains fairly constant for a basin, although the coefficient k varies by orders of magnitude across recession events [47]. The median of the α values observed across the recession events can be considered as the representative α (α r ) of the basin [11,48]. Once α r is determined for a basin, k can be computed for each recession event separately by fixing α at α r . This is essential because α and k are typically highly correlated as the unit of k is dependent on α, and thus, α has to be fixed before computing [54,55]. Although the value of α r is generally close to 2 [13], large variations are observed sometimes [59]. However, for the sake of simplicity we assume here that α r = 2 for all the study basins. (a) There is almost zero correlation between dS an /dt and Q due to high observational errors associated with dS an /dt, suggesting that dS an /dt cannot be used to predict Q. Although Q vs. S an plot exhibits less scatter, the correlation is not very strong, likely because of the fact that the relationship is not one-to-one as shown in (b). Note that the figure was prepared using data from the sample basin of Bluestone river near Pipestem, West Virginia with USGS id 03179000 having drainage area of 395 square miles. Figure 2b shows the Q vs. S an plot for the sample basin. Although the correlation of this plot is better than that of the dS an /dt-Q plot, there is still a large amount of scatter, which is because the relationship between S an and Q can change from one recession event to another [33,48,49,[52][53][54][55]. Figure 3 essentially shows that although the relationship S an and Q can be strong for a recession event, the relationship may vary considerably across events. It essentially suggests that we cannot construct a single mathematical relationship between S an and Q for a basin to predict streamflow. We need to first characterize the dynamic relationship nature of storage-discharge relationship. One avenue in this regard is proper characterization of the relationship between dQ/dt and Q during recession periods following the method proposed by Brutsaert and Nieber [56]. Generally, the relationship between dQ/dt and Q can be expressed as −dQ/dt = kQ α , from which we can easily determine the storage-discharge relationship of the basin [48,57,58]. Biswal and Marani noticed that the value of α generally remains fairly constant for a basin, although the coefficient k varies by orders of magnitude across recession events [47]. The median of the α values observed across the recession events can be considered as the representative α (α r ) of the basin [11,48]. Once α r is determined for a basin, k can be computed for each recession event separately by fixing α at α r . This is essential because α and k are typically highly correlated as the unit of k is dependent on α, and thus, α has to be fixed before computing [54,55]. Although the value of α r is generally close to 2 [13], large variations are observed sometimes [59]. However, for the sake of simplicity we assume here that α r = 2 for all the study basins.
It can be reasoned that k is primarily influenced by the initial storage in the basin [12,13]. Following this, we attempt here to relate k with GRACE based storage anomaly, S an . Since S an has negative values too, we modified S an by subtracting maximum negative S an from the storage time series of each basin and added 1 to the time series. Since storage depletes gradually, we can expect k to be influenced by past storage states of the basin [12,38]. In other words, k can be expressed as a function of past S an , and S anPN , the mean S an from N to 2 days before the hydrograph peak. The analysis was carried out following the least square linear regression method. Regression analysis was carried out between dQ/dt and Q to get k and α of the Brutsaet-Nieber equation −dQ/dt = kQ α . The correlation between k and S anPN (indicated in terms of the coefficient of determination R 2 anPN ) is expected to weaken with N. This logic was first proposed by Biswal and Nagesh Kumar [12], who related k with Q PN , mean discharge from N to 2 days before the hydrograph peak, a proxy for past storage. The correlation between k and Q PN is indicated in terms of R 2 QPN .
Geosciences 2020, 10, 404 6 of 15 recession events can be considered as the representative α (α ) of the basin [11,48]. Once α is determined for a basin, k can be computed for each recession event separately by fixing α at α . This is essential because α and k are typically highly correlated as the unit of k is dependent on α, and thus, α has to be fixed before computing [54,55]. Although the value of α is generally close to 2 [13], large variations are observed sometimes [59]. However, for the sake of simplicity we assume here that α = 2 for all the study basins. Although the correlations are strong (R 2 in each case is greater than 0.96), the relation between the two variables is clearly not unique (d), which supports the notion that Q-S an relationship is dynamic, not static. Note that to limit the influence of observational errors in the analysis we selected only the recession curves for which both S an and Q continuously decreased over time.
Further, we also obtained mean past storage from N to 2 days from CLSM (S clPN ) for each recession event and computed the coefficient of determination between k and S clPN , denoted as R 2 clPN . In this study, we considered the following values of N: 10, 30, 60, and 120 days ( Figure 4). Finally, we compared R 2 anPN with R 2 QPN and R 2 clPN to assess the relative importance of S an in predicting k.

Results and Discussion
Most studies assume that the relationship between basin-scale storage and discharge, which is characterized by the power-law coefficient k, to be unique for a basin [29,58]. Typically, k is related to some observable basin properties for hydrological modelling purposes [38]. However, our study supports the hypothesis that k can vary significantly across recession events as done by Biswal and Nagesh Kumar [12], which is revealed by the fact that for the sample basin 08153500, where, k varies by several orders of magnitude ( Figure 4).
We showed that k has a power-law relationship with basin storage represented by S anPN . R 2 anPN decreased with N in most cases, which was expected. Figure 4 demonstrates k vs. S anPN plots for a sample basin. This phenomenon essentially suggests that the effect of storage on streamflow diminishes with time [12,58]. We therefore consider only N = 10 for the subsequent

Results and Discussion
Most studies assume that the relationship between basin-scale storage and discharge, which is characterized by the power-law coefficient k, to be unique for a basin [29,58]. Typically, k is related to some observable basin properties for hydrological modelling purposes [38]. However, our study supports the hypothesis that k can vary significantly across recession events as done by Biswal and Nagesh Kumar [12], which is revealed by the fact that for the sample basin 08153500, where, k varies by several orders of magnitude ( Figure 4).
We showed that k has a power-law relationship with basin storage represented by S anPN . R 2 anPN decreased with N in most cases, which was expected. Figure 4 demonstrates k vs. S anPN plots for a sample basin. This phenomenon essentially suggests that the effect of storage on streamflow diminishes with time [12,58]. We therefore consider only N = 10 for the subsequent analysis. Figure 5 shows k vs. S anPN plots along with k vs. Q PN plots and k vs. S clPN plots for three sample basins for N = 10. All the plots exhibit power-law relationship. It essentially suggests any data that provides information on basin storage can be used to predict recession coefficient k. Considering all the study basins, we found that the 25th, 50th, and 75th percentiles are (0.46, 0.56, and 0.66), (0.28, 0.39, and 0.50), and (0.44, 0.50, and 0.60), respectively, for R 2 QP10 , R 2 anP10 , and R 2 clP10 (see Figure 6a-c). Performance of S anP10 (Figure 6b) seems to suggest that GRACE is useful in providing information on basin storage for streamflow modelling even though the data is available at such a coarse (1 • × 1 • ) resolution. Moreover, both Q P10 and S clP10 appear to be more useful than S anP10 in predicting recession coefficient k (Figure 6a-c), indicating that GRACE data is more erroneous. In other words, better estimation of k can be done if we are provided with more accurate S an data. To test this hypothesis, we continued our analysis for the recession events during which both Q and S an continuously declined with time. Note that for 18 study basins only we found more than seven recession events satisfying the above criterion. Figure 6d,e shows the boxplots considering R 2 anP10 values for the 18 study basins. The 25th, 50th, and 75th percentiles of R 2 anP10 improved from (0.27, 0.43, and 0.54) (Figure 6d) to (0.36, 0.54, and 0.78) (Figure 6e), respectively, which seems to be supporting the notion that better S an data can lead to better prediction of k. However, even according to Figure 6e R 2 anP10 is close to zero for some of the study basins, implying that GRACE data is not useful for estimation of recession coefficients in every basin. It might be also possible that the relationship between past storage and recession coefficient is not strong for those basins because of quick depletion of storage from groundwater reservoirs [12].  Note that all the plots were obtained for N = 10. For each sample basin, R 2 anPN was found to be greater than R 2 clPN but less than R 2 QPN , suggesting that although S anPN data is capable of explaining storage dynamics, it is associated with significant observational errors for smaller basin. Note that all the plots were obtained for N = 10. For each sample basin, R 2 anPN was found to be greater than R 2 clPN but less than R 2 QPN , suggesting that although S anPN data is capable of explaining storage dynamics, it is associated with significant observational errors for smaller basin. Geosciences 2020, 10, x FOR PEER REVIEW 9 of 16 It should be highlighted that GRACE data is known for being associated with large observational uncertainties [60], which is the reason why most studies use monthly average S an . That is, daily GRACE based S an information is expected to be associated with even larger observational errors. Furthermore, S an data with 1° × 1° spatial resolution (a single pixel covers an area 12,321 sq. km) was obtained from the original S an dataset with 500 km spatial resolution [24], which means its use for small to medium sized basins is not very appropriate. This is why it is recommended to use GRACE based S an data for large river basins with drainage area greater than 200,000 km 2 [32,61]. GRACE based S an data was also used for relatively smaller basins and observed that the agreement between S an and discharge declines with decreasing drainage area [62]. Nevertheless, despite the above limitations related to GRACE based S an data quality, our results highlight the potential of gravimetric satellite missions to provide crucial information on how the power-law coefficient varies across recession events. One possible reason behind this could be that storage does not fluctuate a lot in space and time. This is perhaps the reason why prediction of k in a basin can be performed even using S an data from a nearby basin. It should be kept in mind that recession flow is a result of complex interactions between various hydrological phenomena such as rainfall, evapotranspiration, overland channel flow, and groundwater flow. Nevertheless, despite the complexities, it seems basin scale hydrology can be described quite well by simply understanding the relationship between storage and discharge. In connection with this point, our results further strengthen the earlier notion that past storage information can be exploited to predict k of a basin. This is because a basin is expected to drain groundwater slowly over a long period of time, and hence past storage states can provide crucial information about the coefficient k of an event [13,29,37,38].
The advantage of GRACE based storage anomaly data is that it can be used to predict k, an important variable in hydrological modelling, without using streamflow data once the relationship between S anPN and k is known. One may thus wonder if there is a unique or one-to-one relationship between S anPN and k in which case we can predict k even for ungauged basins, i.e., basins that do not have streamflow gauging stations. Figure 7 shows k vs. S anP10 plot considering all the recession events from all the study basins. It is quite apparent that k vs. S anP10 data points from individual basins show appreciable correlation, although the combined k vs. S anP10 plot is characterized by very high scatter. It can be clearly noticed that (S anPN , k) data points from different basins show very different relationships, meaning there is no unique relationship between S anPN and k. This is perhaps because the value of S av (see Equation (1)) is not constant everywhere, and, therefore, S an cannot independently provide information on S, which actually controls the value of k. However, S av is not It should be highlighted that GRACE data is known for being associated with large observational uncertainties [60], which is the reason why most studies use monthly average S an . That is, daily GRACE based S an information is expected to be associated with even larger observational errors. Furthermore, S an data with 1 • × 1 • spatial resolution (a single pixel covers an area 12,321 sq. km) was obtained from the original S an dataset with 500 km spatial resolution [24], which means its use for small to medium sized basins is not very appropriate. This is why it is recommended to use GRACE based S an data for large river basins with drainage area greater than 200,000 km 2 [32,61]. GRACE based S an data was also used for relatively smaller basins and observed that the agreement between S an and discharge declines with decreasing drainage area [62]. Nevertheless, despite the above limitations related to GRACE based S an data quality, our results highlight the potential of gravimetric satellite missions to provide crucial information on how the power-law coefficient varies across recession events. One possible reason behind this could be that storage does not fluctuate a lot in space and time. This is perhaps the reason why prediction of k in a basin can be performed even using S an data from a nearby basin. It should be kept in mind that recession flow is a result of complex interactions between various hydrological phenomena such as rainfall, evapotranspiration, overland channel flow, and groundwater flow. Nevertheless, despite the complexities, it seems basin scale hydrology can be described quite well by simply understanding the relationship between storage and discharge. In connection with this point, our results further strengthen the earlier notion that past storage information can be exploited to predict k of a basin. This is because a basin is expected to drain groundwater slowly over a long period of time, and hence past storage states can provide crucial information about the coefficient k of an event [13,29,37,38].
The advantage of GRACE based storage anomaly data is that it can be used to predict k, an important variable in hydrological modelling, without using streamflow data once the relationship between S anPN and k is known. One may thus wonder if there is a unique or one-to-one relationship between S anPN and k in which case we can predict k even for ungauged basins, i.e., basins that do not have streamflow gauging stations. Figure 7 shows k vs. S anP10 plot considering all the recession events from all the study basins. It is quite apparent that k vs. S anP10 data points from individual basins show appreciable correlation, although the combined k vs. S anP10 plot is characterized by very high scatter. It can be clearly noticed that (S anPN , k) data points from different basins show very different relationships, meaning there is no unique relationship between S anPN and k. This is perhaps because the value of S av (see Equation (1)) is not constant everywhere, and, therefore, S an cannot independently provide information on S, which actually controls the value of k. However, S av is not expected to vary considerably within a small region, which means if two basins are not situated far apart, we can safely assume them to exhibit the same S anPN -k relationship. The basis for this assumption is that hydrological characteristics do not fluctuate that much in space. For example, Figure 1 shows that mean precipitation does not fluctuate so much within a GRACE pixel. assumption is that hydrological characteristics do not fluctuate that much in space. For example, Figure  1 shows that mean precipitation does not fluctuate so much within a GRACE pixel.
Following the reasoning above if we assume that S anPN values from both basins have appreciable correlation, we can use S anPN time series of one basin to predict k for the other basin. The consequence of this assumption is that we can predict k for an ungauged basin using S anP10 from a nearby gauged basin. The idea is implemented here to predict k of a recession event for an ungauged basin using its S anP10 value for the event in the k vs. S anP10 regression equation borrowed from the nearby gauged basin. Note that a similar hypothesis was tested by Varaprasad et al. who used past discharge data from nearby gauged basins to predict recession coefficients [29]. The main premise is that storage does not fluctuate a lot in space. Figure 8a shows the predicted k vs. observed k plot for a sample basin considered as pseudo-ungauged with the correlation coefficient (R anP10 2 * ) which compares well with its R anP10 2 (i.e., when the basin is considered as a gauged basin). The boxplot for the R anP10 2 * considering results from all the study basins is shown in Figure 8b . Overall, the results here suggest that GRACE based storage anomaly data can be used to predict recession coefficients in completely ungauged basins. Figure 7. k vs. S an10 scatter plots for four sample basins (Basin ids 0317000, 07268000, 03463300, and 03237500) shown together. The idea is to imply that k-S anPN relationship is not universal, and, therefore, it is not possible for us to use S anPN data for predicting streamflow in ungauged basins. Figure 7. k vs. S an10 scatter plots for four sample basins (Basin ids 0317000, 07268000, 03463300, and 03237500) shown together. The idea is to imply that k-S anPN relationship is not universal, and, therefore, it is not possible for us to use S anPN data for predicting streamflow in ungauged basins.
Following the reasoning above if we assume that S anPN values from both basins have appreciable correlation, we can use S anPN time series of one basin to predict k for the other basin. The consequence of this assumption is that we can predict k for an ungauged basin using S anP10 from a nearby gauged basin. The idea is implemented here to predict k of a recession event for an ungauged basin using its S anP10 value for the event in the k vs. S anP10 regression equation borrowed from the nearby gauged basin. Note that a similar hypothesis was tested by Varaprasad et al. who used past discharge data from nearby gauged basins to predict recession coefficients [29]. The main premise is that storage does not fluctuate a lot in space. Figure 8a shows the predicted k vs. observed k plot for a sample basin considered as pseudo-ungauged with the correlation coefficient (R 2 * anP10 ) which compares well with its R 2 anP10 (i.e., when the basin is considered as a gauged basin). The boxplot for the R 2 * anP10 considering results from all the study basins is shown in Figure 8b. The 25th, 50th, and 75th percentile values of R 2 * anP10 are, respectively, 0.31, 0.40, and 0.44. Figure 8c plots R 2 anP10 − R 2 * anP10 values for all the basin to show how R 2 * anP10 compares with R 2 anP10 . Overall, the results here suggest that GRACE based storage anomaly data can be used to predict recession coefficients in completely ungauged basins. values from all the study basin.

Summary and Concluding Remarks
Basin storage is an elusive entity of which direct observation is practically infeasible at present. For hydrological modelling purposes, therefore, information on basin storage is commonly obtained following a suitable analytical method. Once storage is estimated, it can be used to predict discharge

Summary and Concluding Remarks
Basin storage is an elusive entity of which direct observation is practically infeasible at present. For hydrological modelling purposes, therefore, information on basin storage is commonly obtained following a suitable analytical method. Once storage is estimated, it can be used to predict discharge by constructing the relationship between storage and discharge. However, storage-discharge analysis is typically associated with several uncertainties, and, therefore, any information on storage obtained through measurement is expected to be very helpful in hydrological modelling. In this regard, storage anomaly information (S an ) provided GRACE satellite mission has been proven itself to be quite useful. The main aim of this study was to use GRACE based S an data to predict the power law recession coefficient k that effectively characterizes basin-scale storage-discharge relationship. In particular, we attempted to predict k using past S an information, S anPN , i.e., mean S an from N to 2 days before the recession event. Results here seem to suggest that GRACE based S an information can explain variation of k across recession events. The relationship between S anPN and k was observed to be declining with N, which supports the notion that the effect of storage k declines with time. Furthermore, we observed that S anP10 is reasonably good at explaining variation of k. The correlation between S anP10 and k further improved when we considered only the recession curves for which both Q and S an declined continuously, which means there is scope to improve k prediction if more accurate satellite based S an information is provided. Many studies have used different land surface models to enhance data quality derived from GRACE and have successfully used them for hydrological analysis [63][64][65]. Our study can further be improved by assimilating GRACE derived storage data with CLSM storage data. Another limitation of GRACE based S an data is that we cannot obtain the value of average storage S av for a region, and thus it is not possible to construct a universal relationship between S anPN and k. However, results here imply we can safely assume that both S an and S av do not vary so much in space, in which case k can be estimated for a recession event in an ungauged basin by putting the S anP10 value in the k-S an regression equation of a neighboring gauged basin. Overall, our study demonstrated the potential of GRACE based S an information to explain the dynamic nature of the storage-discharge relationship.

Acknowledgments:
We would like to acknowledge that a portion of this work was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. We thank three anonymous reviewers for providing us very valuable comments and suggestions, leading to a significant improvement of the paper.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Coefficient of correlation between k (recession constant) and past discharge and between k and past total water storage from GRACE. Column 1 gives basin ID. Columns 2 and 3 are the latitude and longitude of the gauging station of the basin. The 4th column is the drainage area of basin in square miles. The 5th column is the correlation between k and past average discharge between 2 and 10 days of discharge from peak discharge. The 6th column is the correlation between k and corresponding past average total water storage from GRACE.