the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Scaling methods of leakage correction in GRACE mass change estimates revisited for the complex hydroclimatic setting of the Indus Basin
Vasaw Tripathi
Andreas Groh
Martin Horwath
Raaj Ramsankaran
Total water storage change (TWSC) reflects the balance of all water fluxes in a hydrological system. The Gravity Recovery and Climate Experiment/FollowOn (GRACE/GRACEFO) monthly gravity field models, distributed as spherical harmonic (SH) coefficients, are the only means of observing this state variable. The wellknown correlated noise in these observations requires filtering, which scatters the actual mass changes from their true locations. This effect is known as leakage. This study explores the traditional basin and grid scaling approaches, and develops a novel frequencydependent scaling for leakage correction of GRACE TWSC in a unique, basinspecific assessment for the Indus Basin. We harness the characteristics of significant heterogeneity in the Indus Basin due to climate and humaninduced changes to study the physical nature of these scaling schemes. The most recent WaterGAP (Water Global Assessment and Prognosis) hydrology model (WGHM v2.2d) with its two variants, standard (without glacier mass changes) and Integrated (with glacier mass changes), is used to derive scaling factors. For the first time, we explicitly show the effect of inclusion or exclusion of glacier mass changes in the model on the gridded scaling factors. The inferences were validated in a detailed simulation environment designed using WGHM fields corrupted with GRACElike errors using full monthly error covariance matrices. We find that frequencydependent scaling outperforms both basin and grid scaling for the Indus Basin, where mass changes of different frequencies are localized. Grid scaling can resolve trends from glacier mass loss and groundwater loss but fails to recover the small seasonal signals in trunk Indus. Frequencydependent scaling can provide a robust estimate of the seasonal cycle of TWSC for practical applications such as regionalscale water availability assessments. Apart from these novel developments and insights into the traditional scaling approach, our study encourages the regional scale users to conduct specific assessments for their basin of interest.
 Article
(8548 KB)  Fulltext XML

Supplement
(372 KB)  BibTeX
 EndNote
The terrestrial water storage (TWS) includes all water components on and underneath the Earth's surface, i.e., soil moisture, surface water, groundwater, snowpack, and the water contained in biomass (Zhang et al., 2016). Regionalscale hydrological studies mainly deal with terrestrial water storage changes (TWSCs) over time. In situ measurement of TWSC from its components is practically impossible at basin scales. This is due to the inability of current observational methods to include all possible water storage compartments and to the pointscale nature of existing measurements, which does not capture the spatial variability of TWSC in a basin. With the launch of the Gravity Recovery and Climate Experiment (GRACE) satellite mission in 2002, globalscale observation of TWSC was made possible. These observations come at a monthly timescale, making GRACE an invaluable tool to study seasonal mass changes significant to hydrology (Jiang et al., 2014) and more extended timescales required for climate change studies (Tapley et al., 2019). However, the spatial resolution is of the order of several hundred kilometers, so GRACE is most accurate and valuable in global or continentalscale studies only. The limited spatial resolution is majorly contributed by the errors arising from the GRACE measurement process and instruments, modeling deficiencies, and background model errors used in estimating gravity field parameters (Flechtner et al., 2016). Measurement of intersatellite ranges along the orbit causes the range rate observations to be more sensitive to mass changes in this direction, resulting in error correlations manifesting as north–south stripes in maps. Satellite altitude (∼ 450 km) and intersatellite distance (∼ 220 km) cause these errors to be even more prominent at smaller spatial scales, limiting effective spatial resolution.
These errors are addressed by filtering the data, including destriping filters to remove striping and low pass filters to reduce random errors in small spatial scales. Destriping filters followed by low pass filtering are relatively simple, making them attractive for many users and performing well overall (Klees et al., 2008). However, the inevitability of filtering leads to additional uncertainties in TWS estimates arising from signal attenuation and leakage. Leakage errors occur from truncation of maximum degree (due to all inevitable limitations inherent to GRACE solutions) in the spherical harmonic model along with additional filtering, which leads to mass changes in the region of interest (ROI) to be affected by mass changes outside the ROI and vice versa.
Signal restoration and leakage correction have been an active area of research in hydrogeodetic communities, and have been traditionally carried out using three hydrological or land surface modeldependent approaches; the additive correction approach, the multiplicative correction approach, and the scaling factor approach (Long et al., 2015). In addition, the data driven correction (DDC) approaches (Vishwakarma et al., 2017) are also routinely applied now as a modelindependent method of leakage estimation and correction at basin scale (for basins with an area greater than ∼ 65,000 km^{2}). The scaling factor approach has been the most widely used for leakage correction. Its simplicity in application to the gridded TWS products (1^{0}×1^{0}) popularized and revolutionized the use of GRACE data in the hydrological community. The GRACE Tellus website (GRACE Tellus, 2019) provides gridded scaling factors calculated from the Community Land Model (CLM4.0) that must be applied to the GRACE grids as a regular procedure. The scaling factor approach (Landerer and Swenson, 2012) uses simulated TWSC from global hydrology models (GHMs) or land surface models (LSMs) and processes them in the same way as GRACE to obtain filtered simulations. A scaling factor is obtained through the least squares fit between original and filtered model simulations, which is multiplied to GRACE estimates to account for leakage.
Scaling factors have been explored in multiple schemes. Spatially, they can either be lumped, i.e., a single scaling factor for the basin, or distributed, i.e., gridded at a specific resolution. A timescaledependent scheme uses different scaling factors for mass changes occurring at different timescales. While basin scaling and gridded scaling are the most popular schemes due to their simplicity and globally acceptable performance, the timescaledependent approach is far less studied. The timescaledependent (or frequencydependent) scaling has been suggested on an application basis in previous studies (Rodell et al., 2009; Landerer and Swenson, 2012; Velicogna and Wahr, 2013), but has not been adequately assessed as a scaling method along with traditional methods. Hsu and Velicogna (2017) used a different scaling factor for seasonal and trend components at each grid cell to determine land water storage contribution to sealevel change. However, the properties of such a scaling were not discussed further. In any case, it is established that a single scheme cannot be ascertained to perform well for all regions (Long et al., 2015). In the absence of in situ observations of TWSC, no single criterion for comparison and performance assessment of these different schemes can be utilized. Therefore, we design a simulation environment using Terrestrial Water Storage Anomaly (TWSA) from a hydrological model corrupted with GRACElike errors to validate the inferences made in a regionspecific intercomparison to establish the benefits and drawbacks of each of these schemes.
The sensitivity of the scaling factor approach to the choice of LSM or GHM has been established in numerous studies by deriving and comparing scaling factors from different models (Huang et al., 2019). This sensitivity arises from the difference in underlying model physics, modeled water storage compartments, and the accuracy of forcing datasets used in these models. Most LSMs and GHMs do not model the effect of human intervention on water storage changes. These human interventions include irrigation water use, groundwater use, and artificial reservoir storage, and affect the distribution of TWS changes through complex feedbacks between climateinduced and humaninduced changes (Döll et al., 2003). The Water Global Assessment and Prognosis (WaterGAP) hydrology model (WGHM) is one such model that accounts for the human intervention in modeling TWS changes (Döll et al., 2003; Müller Schmied et al., 2016).
Moreover, an integrated version of the WGHM and global glacier model (GGM) (Marzeion et al., 2012) has been recently produced (Cáceres et al., 2020), which includes glacier mass changes in the modeled TWSC. Hence, models that include these complex interactions in deriving the scaling factors promise a more effective leakage correction in GRACE estimates. Comparing gridded scaling factors from two versions allows us to quantify the change in the spatial distribution of the factors when the model excludes a water storage compartment. Indus Basin is chosen as the study region due to its complex hydroclimatic nature, which will provide an opportunity to explore the different scaling schemes and two versions of the model with and without the glacier mass changes.
The objectives of this study are (1) to present simple processing of GRACE level 2 data to obtain TWS changes for the Indus Basin, (2) to quantify the amount of leakage error and develop scaling factors under three different schemes, and (3) to evaluate the schemes through simulation experiments and residual leakage and scaled GRACE noise levels. Through these objectives, we highlight the need for a basinspecific assessment, develop a novel frequencydependent scheme, and show the effect of including or excluding a crucial water storage component in the model used for deriving scaling factors.
2.1 Study area
The Indus River basin, shown in Fig. 1 (basin boundaries by ICIMOD, 2021), covers an area of 1.14 million km^{2}. The basin spans over four nations, India, Afghanistan, Pakistan, and China, and supports over 215 million people with an approximate water availability of 1329 m^{3} per head (Frenken, 2012). The Indus River is a perennial river originating from the Bokar Chu glacier in Mt. Kailash in Tibet. It flows through the Ladakh region in Kashmir, fed by the glaciers of the Himalayas, Karakoram, and the Hindu Kush ranges, entering the plains of Punjab in Kalabagh, Pakistan. In its journey, it is joined by the Kabul River from the west and Panjnad (Jhelum, Ravi, Chenab, Sutlej, and Beas) from the east to drain into the Arabian Sea. Almost 50 % of its discharge is due to snowmelt (Shrestha et al., 2019). Due to the complex terrain conditions, this basin is a datascarce basin which makes monitoring of hydrological variables sparse and inconsistent. Being a transboundary basin, cooperation between nations also aggravates this situation. Hence, given that the TWSC estimates from GRACE are the only source for the entire Indus Basin, their uncertainty due to leakage and possible correction approaches must be studied thoroughly.
The winter climate over the Indus Basin is dominated by western disturbances embedded with the Indian winter monsoon. During the summer, the Indian monsoon (ISM) brings precipitation to the southern parts of the basin (Dimri et al., 2019). The mean annual rainfall varies between 90 and 500 mm in the downstream and midstream segments, while there is more than 1000 mm in the upstream catchments. The climate in the Indus Basin varies from subtropical arid and semiarid to temperate subhumid in Sindh and Punjab, to alpine in the mountainous highlands in the north, with average temperatures ranging between 2 and 49 ^{∘}C (Dimri et al., 2019). The hydroclimatic parameters in the upper Indus Basin are primarily influenced by glacier melt and in the lower part by human intervention in the form of different irrigation schemes and extensive groundwater depletion (Rodell et al., 2009).
Hence, a complex hydrological regime exists in the entirety of the Indus Basin due to rapidly varying climatological and anthropogenic conditions across the basin. Figure 2, obtained from the Copernicus Global land cover service, shows the land cover distribution in the Indus Basin for 2016. Heavily irrigated areas lie along the river course and the central–east region of the basin (Cheema and Bastiaanssen, 2010). Due to the reduction in precipitation and increase in demand, groundwater has been the primary source of irrigation. The groundwater depletion in Punjab and the national capital region (just outside the Indus Basin boundaries) is a serious problem leading to severe localized land deformations (Garg et al., 2022). Due to rising mean annual temperatures, the evapotranspiration (ET) is also high in these regions. Hence, human intervention is the dominant shortscale driver of TWS changes in these regions, as opposed to the upper Indus Basin, where glaciers and snowmelt dominate the TWS changes.
2.2 Datasets
2.2.1 GRACE data
The GRACE Level 2 data consisting of monthly Stokes coefficients of Earth's geopotential was used to derive estimates of TWSC over the Indus Basin. The data from three Science Data Centers, Jet Propulsion Laboratory (JPL), University of Texas, Center for Science and Research (CSR), and GeoForschungsZentrum (GFZ) in release 6 (RL06) were utilized up to degree 90. Only the months common to all three solutions and WGHM output were used. This resulted in 147 monthly solutions from April 2002 to December 2016. The GSM (denoted as is) products were used, containing fully normalized geopotential coefficients representing the full magnitude of land hydrology, ice, and solid Earth processes. In addition, the RL06 MASCON product from CSR was also utilized to compare the overall results of level 2 processing. The ICE6G_D model was used to account for glacial isostatic adjustment (Peltier et al., 2018). This model is available up to degree and order 256 and provides secular rate of change of gravity field in mm yr^{−1}. The degree 1 coefficients from (Sun et al., 2016), distributed as Technical Note (TN13) on the Tellus website, were used. The replacement C_{20} coefficients from (Cheng and Ries, 2017), distributed as TN11 on the Tellus website, were used. To derive GRACElike errors for the design of simulation experiments, the monthly normal equations provided by TU Graz for the ITSG2018 solutions (Kvas et al., 2019) in SINEX format were used.
2.2.2 Model data
TWS anomalies from WGHM 2.2d, both standard and integrated versions, were used (Döll et al., 2003; Müller Schmied et al., 2016; Cáceres et al., 2020). The integrated version simulates glacier mass changes, unlike the standard version. Grids of 0.5^{∘} resolution from April 2002 to December 2016 were used and resampled to 1^{∘}. For each version, there were four variants available with two climate forcings; ERAInterim reanalysis (WFDEI) applied with WATCH Forcing Data (WFD) methodology and biascorrected using precipitation data sets derived from rain gauge observations of either GPCC v5/v6 (Global Precipitation Climatology Centre; Schneider et al., 2015) or CRU TS3.10/TS3.21 (Climate Research Unit; Harris et al., 2014) and two irrigation variants; irrigation water use at 70 % of consumptive use (CU) and optimal CU. No single irrigation scenario can be assigned for the Indus Basin due to the heterogeneity in irrigation patterns. Hence, an average of these four variants was taken under each version (standard and integrated) to obtain the respective grids. Table 1 summarizes the variants of the model used in this study.
WGHM simulates water flows at a daily timescale on a 0.5^{∘} by 0.5^{∘} global grid, excluding Antarctica. Human intervention and its effect on water flows are also considered by simulating human water use in five sectors: irrigation, manufacturing, livestock farming, domestic use, and thermal power plants. The model requires daily meteorologic inputs of nearsurface air temperature, precipitation (rainfall and snow), downward shortwave radiation, and downward longwave radiation. Reservoir data come from the Global Reservoir and Dam (GRanD) database, which includes 6862 reservoirs with a total storage capacity of 6197 km^{3} (Lehner et al., 2011).
Water balance is carried among three water storage compartments in the vertical direction: canopy, snow, and soil water storage. WaterGAP represents soil as a onelayer soil water storage compartment characterized by a land cover and soilspecific maximum storage capacity and soil texture. The groundwater is recharged from the soil, and storage is simulated after accounting for net abstractions from groundwater due to human use. A fraction of groundwater recharge is assumed to flow back to surface water bodies, and groundwater recharge under surface water bodies is neglected. The surface water storage comprises five subcompartments: local lakes, local wetlands, global lakes, reservoirs, and global wetlands. At each stage, the storage is simulated by accounting for ET losses and net abstractions from lakes and reservoirs. Finally, the output from surface water bodies and groundwater flows into rivers, where the river water storage is simulated after accounting for the streamflow.
Groundwater surface water use (GWSUSE) model accounts for net abstractions from surface water and groundwater. Irrigation is assumed to vary with countries (efficiency of irrigation water use) and source of irrigation water. In groundwaterdepleted areas, the CU from irrigation is considered 70 % of the optimal CU since farmers have less availability to satisfy the optimal requirements. In all other areas, it is assumed to be optimal. Countryspecific efficiency values are used for surface water irrigation, while in case of groundwater irrigation, water use efficiency is set to a relatively high value of 0.7 worldwide. These abstractions are accounted for in the water balance for each compartment as described above.
The glacier mass changes not included in the WGHM are obtained from global glacier model (GGM) (Marzeion et al., 2012). GGM consists of a surface mass balance model based on the temperature index approach and model that accounts for changes in glacier geometry in feedback to compute the changes in glacier mass. The model is forced by a mean ensemble of seven atmospheric datasets, reducing the uncertainty inherent in the individual forcing datasets. As an initial condition for glacier geometry, data like glacier areas and boundaries are taken from Randolph Glacier Inventory (RGI v6) (Pfeffer et al., 2014). The model is calibrated using observed surface mass balance from World Glacier Monitoring Service (WGMS).
2.3 Methods
2.3.1 GRACE data processing
The GRACE level 2 SH coefficients from the GSM product represent the geopotential of Earth containing the full magnitude of land hydrology, ice, and solid Earth processes during each month t as in Eq. (1),
where θλ are colatitude and longitude respectively, n,m are degree and order of the expansion, N is the maximum degree used which was 90 in this case, ${P}_{n}^{m}$ are the associated Legendre's functions, ${C}_{n\phantom{\rule{0.125em}{0ex}}m}^{t}$ and ${S}_{n\phantom{\rule{0.125em}{0ex}}m}^{t}$ are the spherical harmonic coefficients at month t, and R_{e} (6378 km) is Earth's equatorial radius.
The mean static gravity field is removed to obtain gravity field anomalies. In this study, the mean period was taken as the mean of the entire study period. Monthly solutions common to all three (JPL, CSR, and GFZ) centers till December 2016 were extracted with a tolerance of 15 d in the middle of each monthly solution epoch. Running means over three consecutive solutions are taken, which reduces the noise in the solutions. Glacial isostatic adjustment (GIA) is the ongoing secular response of Earth's surface to the mass changes that occurred due to the deglaciation process on Earth and is removed using the ICE6G model (Peltier et al., 2018). The same GIA model used in CSR MASCONS is used to maintain consistency. Assuming mass changes occur in a thin spherical shell (∼ 15 km) around Earth (Wahr et al., 1998), these coefficients are multiplied by degreedependent factors to obtain surface mass changes coefficients in terms of equivalent water height (EWH). The resulting change in surface mass density is obtained as Eq. (2),
where ρ_{avg} is the average density of Earth (5515 kg m^{−3}) and ${k}_{n}^{\prime}$ are the elastic load Love numbers of Earth.
Swenson filter (Swenson and Wahr, 2006) is used for destriping, which removes the correlated errors in even and odd degrees. The filter parameters are chosen after Sasgen et al. (2018), which minimizes the signal and noise contamination between midlatitude and polar regions during destriping. A Gaussian low pass filter of half width 300 km is employed to reduce the random errors in higher degree terms. Grids of filtered surface mass changes are obtained as Eq. (3),
where W_{n} are degreedependent factors of Gaussian filter in spectraldomain (Sasgen et al., 2018).
From the filtered GRACE surface density coefficients, mass changes for the Indus Basin are obtained by regional integration using the region function defined in Eq. (4),
where RF is the region function inside the region R of the Indus Basin and region Ω is the total surface area of Earth. This region function is transformed to spectraldomain (up to degree 90), and the corresponding coefficients are multiplied with filtered surface density coefficients to obtain mass changes for the Indus Basin, ΔM, as in Eq. (5),
where R_{c} and R_{s} are respectively the cosine and sine coefficients of the region function in the spectral domain.
2.3.2 Scaling factors determination
WGHM anomalies were centered around the mean of the entire observation period (April 2002–December 2016) by including only GRACE months. Unfiltered mass change time series for the Indus Basin from both versions were obtained using the region function in Eq. (4), by integrating the global grids over the surface of Earth as in Eq. (6),
where $\mathrm{d}\mathrm{\Omega}={R}_{\mathrm{e}}^{\mathrm{2}}\mathrm{sin}\mathit{\theta}\mathrm{d}\mathit{\theta}\mathrm{d}\mathit{\lambda}$ is an element of the area on the region surface. The global grids are then transformed to the spherical harmonic domain up to degree 90 and applied with Swenson destriping and Gaussian low pass filter in the same manner as described in Sect. 2.3.1. The filtered coefficients are transformed back to the spatial domain to obtain filtered model grids. Regional integration in the spectral domain is performed with filtered model coefficients to get filtered model mass change time series.
A single scaling factor (k) for the entire basin is derived by a least squares regression between unfiltered model time series and filtered model time series, that is, by minimizing L in Eq. (7),
over the entire period of study. In Eq. (7), $\mathrm{\Delta}{S}_{t}^{\mathrm{u}}$ is the unfiltered and $\mathrm{\Delta}{S}_{t}^{\mathrm{f}}$ is the filtered model mass for the tth month. A scaling factor for each grid cell in a 1^{0}×1^{0} grid is derived using the unfiltered model grid and filtered model grid and minimizing L in Eq. (7) but at every grid cell in the basin. This is done for both versions of WGHM to obtain a 1^{∘} resolution map of scaling factors.
To derive frequencydependent scaling factors, the total mass changes from unfiltered and filtered model versions are decomposed into longterm linear, seasonal, and residual components as Eq. (8),
For this purpose, a time series model is fit to the data using least squares as in Eq. (9),
where β_{1} is the linear trend, ${\mathit{\beta}}_{k}^{\mathrm{c}}$ and ${\mathit{\beta}}_{k}^{\mathrm{s}}$ are the amplitudes of the cosine and sine terms of the periodic component with a period T_{k}, and t represents the month with respect to the mean of the observation period. However, the periods (T_{k}) of the seasonal terms are unknown. For this study, the unknown periods are found from the data itself using a Lomb–Scargle (LS) periodogram analysis (Scargle, 1982) which allows detection of weak periodic signals in otherwise random, unevenly sampled data.
From this analysis, the peaks were found at annual and semiannual periods. The false alarm probabilities associated with each peak were minimal ($\sim {\mathrm{10}}^{\mathrm{5}}$), which shows that these periods are significant. These periods were used in Eq. (9) for time series decomposition. The trend, annual, and semiannual components are separated in both model versions' unfiltered and filtered time series. Then a least squares fit is carried out for all three components separately to obtain three scaling factors that minimize the expression in Eq. (10),
where k_{j} is the scaling factor for jth component.
Figure 3 presents the schematic of deriving scaling factors (in green boxes) from filtered and unfiltered model time series (in yellow boxes) and applying them to corresponding GRACE time series (in blue boxes) to obtain scaled GRACE estimates (in purple boxes).
2.3.3 Implementation of datadriven correction approach
The DDC approach or the method of deviation (Vishwakarma et al., 2017) uses twicefiltered GRACE fields to determine two correction terms for the regional averages of filtered GRACE fields, namely the leakage and deviation integrals for the entire basin. We used the method codes provided by the authors (and through personal communication with the lead author) to compute basinaveraged DDC corrected time series for the Indus Basin to compare with our scaled basin averaged time series.
2.3.4 Design of simulation experiment
Due to the unavailability of distributed in situ data, we designed a simulation experiment using WGHM and monthly GRACE TWS errors to assess the extent to which the scaling schemes can recover the signal lost to filtering. We use WGHM fields as a proxy for true TWS signals and corrupt them with GRACElike errors (Willen et al., 2022). For this, we first derived the error covariance matrices from the normal equations provided by TU Graz for the ITSG solutions (Kvas et al., 2019) from degree 2 onwards. For the degree 1 coefficients, the corresponding error covariances were computed from an ensemble of 12 different estimates derived following Swenson et al. (2008), arising from four GRACE products (JPL, CSR, and GFZ release 6, ITSG2018) and three GIA models (A et al., 2013; Caron et al., 2018; Peltier et al., 2018). Assuming the errors to follow multivariate random distribution, we derived their random realizations in the SH domain from the covariance matrices using Cholesky decomposition. We then added these errors to the WGHM fields in SH domain and filtered using the same filter used in the study for GRACE (Swenson destriping +300 km Gaussian). These filtered–corrupted WGHM fields now represent GRACElike observations. Finally, using the scaling factors derived in the study, we rescaled these filtered–corrupted WGHM fields as per the respective schemes (Fig. 3) to recover the lost signals as rescaled WGHM fields.
We use GRACE errors from the ITSG solutions for assessing scaled estimates of GRACE solutions in our study since the normal equations from three centers – JPL, CSR, and GFZ – are not publicly available. Moreover, using ITSGbased errors can be justified by the following: our study uses a monthly running average of the three solutions, which means the errors are already reduced. The ITSG solution which has been shown to have less noise than the JPL, CSR, and GFZ solutions (Ditmar, 2022) thus provides a reasonable estimate of the average error content.
2.3.5 Assessment of uncertainty components
Two quantities are computed using WGHM as a proxy to the true signal to assess the uncertainty specific to leakage and scaling. Equation (11) provides the initial leakage due to filtering, and Eq. (12) provides the amount of leakage not covered by the scaling factors, termed residual leakage.
where std represents the standard deviation, RMS represents the root mean square of the time series indexed by t, k represents any one of the scaling factors schemes, ΔG_{t} represents GRACE mass time series, and $\mathrm{\Delta}{S}_{t}^{\mathrm{fc}}$ represents the filtered–corrupted model time series. If the scaling factors work, the residual leakage should be less than the initial leakage. This quantification assumes all the variation in the difference between filtered–corrupted and unfiltered model time series to represent leakage error (making no distinction between individual effects on underlying signal and noise). It is thus meaningful for intercomparison of different schemes but not as a true estimate of leakage error. This standard deviation is multiplied by the ratio of RMS of GRACE and filtered–corrupted model time series to account for the amplitude difference in GRACE and model (Landerer and Swenson, 2012).
The uncertainty estimation of GRACE mass estimates is done following Groh et al. (2019), which uses only the time series of derived GRACE masses to determine noise. The unscaled GRACE time series contains errors from measurement, leakage, GIA, and C_{20}. Since this approach provides the temporally uncorrelated component of total GRACE errors, the effect of scaling on the random component of time series (thus a random component of leakage) can be compared. The uncertainties from GIA models (<1 mm yr^{−1} equivalent sea level) and C_{20} (<0.1 mm yr^{−1} equivalent sea level) coefficients are almost negligible for the Indus Basin and, hence, are not considered (Caron et al., 2018; Blazquez et al., 2018). A high pass filter is applied to the residuals of the GRACE time series (after removing trend and seasonal components) to remove the unmodeled interannual components in the residuals. The filter width is chosen to be 18 months (6 sigma width of Gaussian hat =18 months), which means all the signal with a larger than 18month period is removed, leaving temporally uncorrelated errors referred to as noise. This noise contains the GRACE measurement error and the random component of the initial leakage error. A scaling factor generated from high pass filtering random white noise signals is multiplied to account for the damping of any white noise component during high pass filtering. Similarly, this process is repeated for scaled GRACE time series. The corresponding scaled noise contains the scaled measurement and residual random leakage errors.
3.1 Mass change estimates from GRACE and WGHM
Figure 4 shows the result of the processing scheme followed to obtain mass change time series from the average of three GRACE level 2 SH solutions and MASCON time series from CSR. The shaded region depicts the 3σ interval of uncertainties from formal errors of the SH coefficients provided by the centers. These uncertainties are propagated to the mass estimates and averaged for three centers by the law of error propagation over Eq. (5) in Sect. 2.3.1. A GRACE trend value of $\mathrm{7.6}\pm \mathrm{0.6}$ Gt yr^{−1} is obtained (Table 2). Similar values of −8.6 Gt yr^{−1} from (Scanlon et al., 2018) and −9.1 Gt yr^{−1} from (Kvas et al., 2019) have been reported. The differences can be attributed to processing strategies and different data releases used in these studies. The performance of the processing scheme is assessed by comparison to the time series obtained from the Level 3 release 6 CSR MASCON (CSRM) solution (Save et al., 2016). Although derived from the same level 1 data as SH solutions, MASCONS are constrained to reduce the leakage effect arising from the postprocessing step. Figure 4 also shows the GRACE SH and CSRM timeseries correlation for the Indus Basin. The high value of Pearson's correlation coefficient (R^{2}=0.95) indicates that the SH solution and processing strategy used are nearly as good as the MASCON solution for the Indus Basin. An R^{2} value of 0.91 is obtained when both time series are detrended, which confirms that the high values are not a spurious result.
Figure 5 shows the mass anomaly time series from the standard and integrated version of WGHM along with the GRACE SHbased time series. A significantly more negative trend is seen in the integrated version ($\mathrm{20.5}\pm \mathrm{0.4}$ Gt yr^{−1}) than in the standard version ($\mathrm{10.3}\pm \mathrm{0.4}$ Gt yr^{−1}), indicative of significant glacier melt in the Indus Basin, which has a dominating effect on the overall TWS trend of the Indus Basin. The negative trend from the standard version can be attributed to the humaninduced changes due to severe groundwater depletion for meeting the irrigation demands in the Indus plains (Rodell et al., 2009) and increasing losses from ET due to climateinduced changes from increasing mean annual temperatures (Shrestha et al., 2019).
The periodogram analysis of the standard and integrated version (Fig. 6b, c) shows clear peaks at annual and semiannual frequencies. The peak frequencies obtained from the model and GRACE (Fig. 6a) are almost identical, confirming the ability of GRACE to observe small seasonal signals in the Indus Basin, which is mainly dominated by trends. However, from a subsequent timeseries fit (Table 2), the annual amplitude (11.8±3 Gt) obtained from GRACE is much smaller than the semiannual amplitude (24.8±2.9 Gt). It is found to be the artifact of filtering, as explained later in Sect. 3.2. Figure 7 shows the spatial distribution of different components of mass changes from WGHM standard and integrated versions. Close to the northern basin boundary, the annual signal content is higher in the integrated version than in the standard version, while the semiannual component is almost the same. This is due to the annual component of glacier mass change contribution absent in the standard version. Respective amplitudes from the model time series fit (Table 2) also reflect this. The semiannual seasonality in TWS changes results from bimodal precipitation distribution in the Indus Basin. The interannual signal content (area of the frequency spectrum in the interannual bandwidth) is almost identical in both model versions and GRACE. This shows that adding glacier mass changes to WGHM does not contribute to the interannual variations. These interannual variations are probably the result of longterm groundwater depletion (Pradhan, 2014) and interannual variability resulting from winter/spring precipitation over the upper Indus Basin due to El Niño–Southern Oscillation (ENSO) (Krakauer, 2019).
3.2 Effect of filtering
Spatially, the effect of filtering is visualized in Fig. 8, which shows mean TWS anomalies for each calendar month from the unfiltered and filtered WGHM Integrated version. Local features and smallscale changes in TWS are smoothed out and reduced in amplitude. Largescale changes in TWS (right panel, Fig. 8) follow the pattern of bimodal precipitation distribution in the Indus Basin due to western disturbances and ISM (Hasson et al., 2016; Hussain et al., 2020). The northwest to southeast increase in storage (reduction of blue color intensity) during the winter months (December to February) and premonsoon months (March to May) can be attributed to western disturbances. It is worth mentioning that there is generally a 1month lag between precipitation and storage in this region. The southeast to the northwest increase of storage (decrease in blue and increase of yellow intensity) during the summer from July to October can be attributed to ISM. October and November show the retreat of monsoon and hence decreasing storage. Smallscale TWS changes (left panel, Fig. 8) are driven by heavy irrigation along the river, southeast Indus plains, and snow and glacier melt in the upper Indus Basin.
Temporally, the filtering dampens the trend, i.e., makes it less negative, in both the model versions (Fig. 9a, b). The dampening is by ∼ 13 % for the standard version and stronger, by ∼ 20 %, for the integrated version. The filtering brings the modelbased trends closer to the GRACEbased trends, partly explaining the trend differences between filtered GRACEbased and unfiltered modelbased time series.
The effect of filtering on seasonal signals can be analyzed from Fig. 6d, e. The peaks are obtained at the same frequencies as the unfiltered time series, which shows that filtering does not induce any additional periodicity of mass changes in the Indus Basin. From the relative peak heights of annual and semiannual frequencies, it is also clear that filtering affects annual and semiannual mass changes differently. The time series fit (Table 2) shows that filtering significantly reduces the annual amplitudes in both versions, while the semiannual amplitudes are reduced only slightly. Hence, both versions show annual mass variations to leak out, explaining why the annual peak from GRACE (Fig. 6a) is suppressed. This is a significant effect of filtering and must be restored by scaling.
3.3 Scaling factors
3.3.1 Basinscale factors
Basin scaling factors of k_{s}=1.14 from the standard version and k_{g}=1.22 from the integrated version are obtained. These values are greater than one, indicating that filtering causes the overall mass changes to leak out; hence, a factor greater than 1 is required to restore the signal. However, the values are small, indicating that the Indus Basin has a small leakage amount overall. Basin scale factors from two studies have reported similar values (Table 3). The addition of glacier mass changes in the model leads only to a minor effect on the basinscale factor. Glacier mass changes are located at the edge of the basin and suffer from more leakage out. This explains the slightly larger scaling factor from the integrated version. The similarity of our results to results from other studies that used different models indicates that basinscale factors are not very sensitive to the model used for the Indus Basin.
3.3.2 Gridded scaling factors
Figure 10 shows the gridded scaling factors in 1^{∘} resolution maps from the standard (left) and the integrated version (right). The maps include 108 scaling factors, one for each pixel inside the Indus Basin. From the standard version, the scaling factors ranged from −0.4 to 10.1. From the integrated version, the scaling factors ranged from −0.5 to 8.5. Significant differences between the two occur in the upper Indus Basin, where the glaciers are located. Table 4 lists the standard interpretation of gridded scaling factor values used in most studies (Long et al., 2015). Figure 11 shows the histogram of scaling factors from both model versions.
From the standard version, scaling factors for 13 grid cells are negative. These grid cells were excluded for scaling, considering outofphase behavior with respect to their neighboring grid cells. Most scaling factors are less than 1 (67 grid cells), out of which most are less than 0.5 (55 grid cells), indicating that a large part of the basin suffers from moderate to significant leakagein. These small scaling factors occur mainly in the lower reaches of the basin (Indus plains) and upper Indus Basin, where nonglacier mass changes are small in magnitude. Few grid cells stand out with larger scaling factors (>3), which depict large local mass changes in the pixel. These occur in the southeast Indus Basin region, which has a larger magnitude of TWS changes due to ISM and human intervention in the form of irrigation and groundwater depletion. These more considerable changes leak out to nearby dry regions of the upper Indus Basin (as represented in the standard version), requiring greater than 1 scaling factor.
From the integrated version, scaling factors at 11 grid cells are found negative and excluded. The number of grid cells with significant leakage in (k<0.5) is larger (62) than for the standard version. The grid cells with large scaling factors (k>3) are more distributed than the standard version. This is because the glacier mass changes cause large local mass variations in the upper Indus Basin, requiring larger scaling factors to account for leakage out to nearby dry regions (northeast of the Indus Basin with arid Tibetan plateau). The spatial variability of scaling factors obtained is higher (CV=1.69) from the integrated version than from the standard version (CV=1.48) due to a more heterogeneous representation of mass changes in the integrated version.
Figure 12 shows the histograms of grid scaling factors only for the nonglaciated region of the Indus Basin (i.e., the nonglaciated pixels in Fig. 1). The increase in the frequency of small scaling factors (<0.5) and decrease in large ones (>3) indicates that the distribution of values from the integrated version is shifted towards zero as compared to the standard version. This shows that the addition of a localized water storage compartment in the model (glacier in this case) not only affects the scaling factors in that region but in the entire basin.
3.3.3 Frequencydependent scaling factors
The frequencydependent scaling factors for the Indus Basin are also shown in Table 3. For the trend component, the scaling factors from both versions (1.14 from standard and 1.23 from integrated) are almost identical to the frequencyindependent basinscale factors. This reinforces that the Indus Basin is dominated by longterm trends rather than seasonal signals, which causes the basinscale factors to be driven by filtering the trend component. The scaling factors for the annual component from both versions (1.64 from standard and 1.59 from integrated) are significantly higher than 1 to account for leakage out of the annual mass changes described earlier. This also shows that the basinscale factor alone would underscale the annual component. The semiannual scaling factors are nearly identical and close to 1 in the standard and integrated version, showing that filtering has a negligible effect on this component. These deviations from a single basinscale factor reinforce the need for frequencydependent scaling factors in basins where mass changes occur at different frequencies and the filtering has a significantly different effect on these individual frequencies.
3.4 Scaled GRACE mass estimates
Table 5 shows the scaled GRACE time series parameters from each scaling scheme along with the corresponding values from an independent DDC approach. Trends become more negative (compared to −7.6 Gt yr^{−1} from unscaled GRACE) from all three scaling schemes. Similar scaled trends from the basin and frequencydependent scaling indicate the dominance of the trend component in the Indus Basin. Grid scaling leads to the most negative trends probably since grid cells with significant local mass changes contributing to the overall trend are scaled more with larger scaling factors than basinaveraged factors. All three scaling schemes restore the mass changes of annual frequency, lost due to filtering, as evident from increased annual amplitude (compared to 11.8 Gt from unscaled). Grid scaling seems to overestimate the annual amplitude when compared to frequencydependent scaling. This overestimation is larger when using the integrated model. Across the schemes, scaling using the integrated model version leads to more negative trends and larger seasonal amplitudes than the standard version, except in the case of semiannual amplitudes from grid scaling, which is reduced (compared to 24.8 Gt from unscaled). We studied these behaviors further in the simulation experiments and discuss their validity in detail in Sect. 3.5.
Different mass change components from gridscaled GRACE from both model versions along with CSR MASCONS are shown in Fig. 13. We can see that even though MASCONS restore some signal, they do it in such a way that signal is added mostly where signal already was in the SH solutions (notice the corresponding pixels in panel 1 and 2 from left). Rescaling, on the other hand, can be seen to offer much more spatial contrast, driven by the spatial distribution of scaling factors. Therefore, rescaling can allow downscaling of GRACE resolution along with leakage correction if an ideal, perfect model were to be used. For example, MASCONS cannot separate the glacier loss trends in the upper Indus pixels from GW depletion trends in pixels in the southern part of Indus, while grid scaling from the integrated model does. However, the tendency of gridded scaling factors to be driven by the dominant mass change component in the pixels (trends or seasonal) may lead to incorrect spatial patterns for one or more of these components. This is discussed further in detail in Sect. 3.6.
The top row of trend maps in Fig. 13 also provides some interesting observations for a phenomenon known as the Karakoram Anomaly through the lens of GRACE. The “Karakoram Anomaly” is termed as the stability or anomalous growth of glaciers in the central Karakoram, in contrast to the retreat of glaciers in other nearby mountainous ranges of the Himalayas and other mountainous ranges of the world (Dimri, 2021). Various mass balances over this region have shown to be balanced or slightly positive (Farinotti et al., 2020). There is, however, significant uncertainty over the reasons for this anomaly which seems to be caused by the absence of reliable in situ observations in this region. However, GRACE observations of such small (nearly stable) trends will be contaminated with the nearby large negative trends due to leakage. This can be seen in Fig. 13 (top panel), where the red bounding box demarcates the extent of Karakoram (Dimri, 2021). The scaled trends for these pixels lying within the Indus Basin (top rightmost) are extremely small but negative (ranging from −0.01 to −0.1 Gt yr^{−1}). The uncertainty associated with these values is, however, nearly double. The two major limitations for the uncertainty of these values come from the inability of GRACE to resolve this extremely small signal from this phenomenon and the inability of the underlying integrated WGHM model to simulate this anomalous behavior at the current 1^{∘} resolution, seen by large negative trends in the corresponding pixels (Fig. 7). However, scaling does seem to make a distinction between pixels of Karakoram with less negative trends and the adjacent pixels of Ladakh region with more negative trends, which cannot be seen in unscaled GRACE fields. We feel this to be a promising result indicating the presence of anomalous behavior in the Karakoram, but any further analysis is out of the scope of this study.
Figure 14 depicts the basinaveraged GRACE time series from grid and frequencydependent scaling using integrated WGHM compared to the DDC time series (Sect. 2.3.3).
Figure 14 shows that overall, frequencydependent scaling using WGHM agrees exceptionally well with the independent DDC method for the Indus Basin. The agreement is also depicted in the time series components in Table 5. The small differences are well within the uncertainty limits and can be attributed to the methods being entirely different. Larger differences are, however, seen between grid scaling and DDC, which highlight the limitation of grid scaling to over or underscale certain pixels. Grid scaling can be seen to overestimate the trends compared to DDC since the trendcontributing pixels get scaled with larger scaling factors. Thus, this effect is seen to be more pronounced in grid scaling from the integrated version than from the standard version. The seasonal amplitudes have large differences, but their nature cannot be generalized between annual and semiannual frequencies. This highlights the differences in physical processes governing the largescale and finescale behavior of these components. These results provide confidence in the new frequencydependent scaling using WGHM to restore the damaged signal at the catchment scale correctly.
A comparison of the gridscaled GRACE SH time series to the MASCON time series (CSRM) is shown in Fig. 15 to judge how realistic the values are. The R^{2} values are similar to those obtained with unscaled GRACE SH, which provides confidence about the scaling factors obtained. We do not use this high value to judge the performance of scaling factors but only to highlight their numerical authenticity. A similarly high R^{2} for the detrended time series shows that the correlation is not spurious.
Similarly, the result of frequencydependent scaling is compared with the MASCON time series in Fig. 16. The obtained R^{2} values are higher than those of grid scaling. This can be attributed primarily to the fact that the noise has not been scaled as in the grid scaling, making it more comparable to the MASCON time series. Even the real signals that may be present in the residual along with noise are not scaled. These real signals are mostly interannual, which are not much affected by filtering (Sect. 3.2), and hence not scaling them leads to a better correlation with the MASCON time series.
3.5 Simulation experiments
The results of the simulations are shown in the form of time series components of the rescaled WGHM fields compared to the original true values in Table 6 for the standard version and Table 7 for the integrated version. The rescaled spatial fields from grid scaling compared to true fields are shown in Fig. 17 for the standard version and Fig. 18 for the integrated version.
The simulation experiments establish the frequencydependent scaling as the best performing scheme in terms of recovering the true basinaveraged signal for both the model versions. The similarity of recovered trends from basin and frequencydependent scaling supports the inference that the basin scaling factors are driven by the dominant mass change component, which is the trend in the case of the Indus Basin. The spatial pattern of the trend component is recovered extremely well by grid scaling. However, at the basin scale, recovered trends from grid scaling are more negative compared to the true trends since the pixels with large trends (Figs. 17 and 18; top left) are scaled with larger scaling factors (Fig. 10). The recovered semiannual amplitude from grid scaling using integrated version is reduced, while increased using the standard version (compared to the corresponding filtered–corrupted semiannual component). This supports the behavior seen in scaled GRACE with the following explanation: most of the semiannual signal contribution comes from nonglaciated pixels in the trunk Indus and few pixels in the upper Indus (Figs. 17 and 18, bottom left). The scaling factors from the integrated version for the trunk Indus pixels are smaller than from the standard version, contributing to a decrease in semiannual amplitude, while at the same time, the scaling factors are larger for the few upper Indus pixels, contributing to an increase (Figs. 17 and 18; bottom right). However, the greater number of trunk Indus pixels outweighs the increase in semiannual signal of a few pixels of the upper Indus, leading to an overall decrease at the basin scale.
It can be seen (Tables 6 and 7) that grid scaling is unable to recover the true annual signal using both model versions, falling short by 29 % in the case of standard and 8 % in the case of the integrated version. Similar to the case made for semiannual amplitude, this is evident from Figs. 17 and 18, where the true spatial pattern of annual component is not recovered by scaling, contributed by greater loss from pixels in trunk Indus and smaller gain from pixels in upper Indus. This contradicts the behavior seen with scaled GRACE estimates where grid scaling overestimates the annual amplitude compared to frequencydependent scaling (equivalent to true amplitude). The above reasoning fails to explain the contradiction. However, comparing the spatial pattern of unscaled GRACE annual amplitude in Fig. 13 (second panel from left) with the pattern of annual amplitude in filtered–corrupted models in Figs. 17 and 18 (middle panel), it can be observed that the pixels from GRACE in the upper Indus Basin already hold much stronger annual signal compared to the corresponding pixels from the filtered–corrupted models. Therefore, larger scaling of these pixels is able to compensate for the loss from trunk Indus pixels, leading to an overall increase at basin scale (again, compared to frequencydependent scaling).
3.6 Residual leakage and scaled noise estimates
Table 8 shows the initial leakage error and the noise level in the unscaled GRACE basin averaged time series containing the random component of leakage. The initial leakage error determined from the standard version is less than that determined from the integrated version due to the smaller magnitude of TWS anomalies. The unscaled noise provides a baseline measure against scaled noise levels to evaluate the effect of scaling on random leakage components.
The scaled noise levels and the residual leakage errors from all three scaling schemes are shown in Table 9. These scaled noises implicitly contain the scaled random component of leakage, as explained in Sect. 2.3.3. Compared to unscaled GRACE (Table 8), the basin and grid scaling cause the noise to be scaled along with the signal, whereas frequencydependent scaling does not affect noise. Grid scaling leads to lower noise than basin scaling due to the suppression of random errors in most of the Indus Basin regions where smaller scaling factors are present. Even smaller gridded scaling factors from the integrated version in those regions may explain the lower scaled noise. Basin scaling from the integrated version leads to higher noise than basin scaling from the standard version due to the higher magnitude of the scaling factor.
The residual leakages (Table 9) are least in frequencydependent scaling, indicating its better performance. Residual leakages from grid scaling indicate their inability to reproduce the spatial pattern of seasonal signals in the basin. We must caution that although the residual leakage from grid scaling using the integrated version is lower, it is only because of the small magnitude of scaling factors. The underlying spatial patterns are no better recovered compared to grid scaling from the standard model.
The study aimed to derive and evaluate scaling factors for the Indus River basin from WGHM 2.2d using its standard and integrated versions, to account for the leakage effects in mass estimates derived from GRACE spherical harmonic solutions. Scaling factors were derived based on three schemes of different spatiotemporal characteristics: basin scaling factors that are spatially and temporally constant; gridded scaling factors that are spatially variable while temporally constant; and frequencydependent scaling factors that are spatially constant and temporally variable. The results of the scaling approach were compared with an independent DDC approach at the basin scale and with CSR MASCONS at both grid and basin scales. Inferences were made in an intercomparison framework evaluated by detailed simulation experiments designed using WGHM spatial fields as a proxy for true TWSA and GRACE errors derived using full error covariance information of the SH coefficients.
The new frequencydependent scaling outperforms others in terms of recovering the true basin average signal (including all its components), as evident by the simulation results and residual leakage levels. Its excellent agreement with an independent DDCbased estimate shows that it can be a viable alternative as a leakage correction approach for the Indus Basin. The frequencydependent scaling, as proposed here, keeps the noise levels unscaled, which can provide robust estimates for practical applications that require accurate knowledge of the seasonal cycle of TWSA, such as water availability studies by decisionmaking bodies to ensure a safe supply of water to a region every year. If the users interested in other basins can determine additional dominant frequencies with sufficient confidence, then the frequencydependent scaling can easily be applied to such frequencies.
Grid scaling allowed the recovery of the spatial pattern of trends but failed to capture the patterns associated with small seasonal amplitudes in the trunk Indus. It was found that with the integrated version, gridscaled GRACE SH fields, unlike CSR MASCONS, could resolve negative trends from glacier mass loss in upper Indus from groundwater loss in the Indus plains. The addition of the glacier mass component modifying the scaling factors of the nonglaciated grid cells along with the glaciated grid cells is an interesting finding of the study. However, it was also found that the characteristics of basin averages of gridscaled GRACE, such as trends and seasonal amplitudes, can be misleading in judging the performance of grid scaling due to over and underscaled pixels compensating, and leading to seemingly realistic basinaveraged values. Thus, we advise our readers to thoroughly assess the behavior of grid scaling for their basin of interest before directly using it for downstream applications.
The basin scaling scheme appears less sensitive to the model's mass distribution, which may justify the use of even worse models in their derivation. For the Indus Basin, basin scaling factors seem to be driven by the filtering effect on the trend and may be used for applications dealing with longterm trends in the basin. However, a better signal–noise separation must be achieved to minimize the scaling of noise. Frequencydependent scaling shows that using a single basin scaling factor for basins like the Indus Basin with mass changes occurring at different frequencies will lead to inappropriate scaling of one or more of these components.
It is obvious to realize the possibility of a fourth gridded and frequencydependent scaling scheme that would provide three scaling factors per grid cell. However, our initial experiments show that the time series decomposition at gridscale leads to vastly varying annual and semiannual components, and causes the scaling factors to be unrealistic. Hence, we leave that development for future studies. The study is limited in providing an external validation with independent TWS estimates. One approach is a comparison with TWS obtained from a water balance, but uncertainties in the derived TWS are much larger than the changes due to scaling. Hence, we stick to an intercomparison framework along with simulation experiments that can be realistically reproduced and applied to different applications and studies that utilize GRACE data in the Indus Basin, and gain insights for choosing an appropriate scaling scheme per requirements. It may be stressed here that although the study has been done only for GRACE data, it can naturally be extended to the ongoing GRACEFO observations (with the availability of the latest model outputs) without any loss of generality.
The MATLAB scripts written to conduct this study can be obtained upon request to the corresponding author.
All GRACE Level 2 release 6 spherical harmonic datasets used in this study are available at http://grace.jpl.nasa.gov (https://doi.org/10.5067/GRGSM20J06, NASA JPL, 2018; https://doi.org/10.5880/GFZ.GRACE_06_GSM, Dahle et al., 2018; https://doi.org/10.5067/GRGSM20C06, UTCSR, 2018), supported by the NASA MEaSUREs Program. CSR release 6 MASCONS were downloaded from http://www2.csr.utexas.edu/grace (Save et al., 2016). The ITSG2018 normal equations for deriving GRACE errors are available at https://www.tugraz.at/institute/ifg/downloads/gravityfieldmodels/itsggrace2018 (Kvas et al., 2019). The output of different WGHM model versions (Cáceres et al., 2020) used in this study can be obtained upon request to the corresponding author. Maps of scaling factors developed in this study can be obtained upon request to the corresponding author. Data required to reproduce all the time series in the article have been provided in the Supplement. Sufficient information has been provided in the article to reproduce the results of this study using the datasets mentioned here.
Data to reproduce all the time series, tabular data, and periodogram plots in the article have been provided in a separate supplementary file. The supplement related to this article is available online at: https://doi.org/10.5194/hess2645152022supplement.
VT and AG designed the study, VT conducted the analysis and wrote the article, MH supervised the GRACE data processing and article editing, and RR supervised the hydrological inferences and article editing. All the authors provided critical feedback, comments, and suggestions for the development of the article.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank Petra Döll and Denise Cáceres (Goethe University Frankfurt) for providing the WGHM v2.2d spatial grids of TWS anomalies and Benjamin D. Gutknecht (Technische Universität Dresden) for helping with their analysis. We also thank the reviewers, A. P. Dimri and Henryk Dobslaw, for their comments that helped us improve this article.
This research has been supported by the German Academic Exchange Service New Delhi (grant no. Combined Study and Practice Stays for Engineers from Developing Countries (KOSPIE) with Indian IITs 2019/20).
This paper was edited by Narendra Das and reviewed by A. P. Dimri and Henryk Dobslaw.
A, G., Wahr, J. and Zhong, S.: Computations of the viscoelastic response of a 3D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada, Geophys. J. Int., 192, 557–572, https://doi.org/10.1093/gji/ggs030, 2013.
Blazquez, A., Meyssignac, B., Lemoine, J., Berthier, E., Ribes, A., and Cazenave, A.: Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets, 215, 415–430, https://doi.org/10.1093/gji/ggy293, 2018.
Cáceres, D., Marzeion, B., Malles, J. H., Gutknecht, B. D., Müller Schmied, H., and Döll, P.: Assessing global water mass transfers from continents to oceans over the period 1948–2016, Hydrol. Earth Syst. Sci., 24, 4831–4851, https://doi.org/10.5194/hess2448312020, 2020.
Caron, L., Ivins, E. R., Larour, E., Adhikari, S., Nilsson, J., and Blewitt, G.: GIA Model Statistics for GRACE Hydrology, Cryosphere, and Ocean Science, Geophys. Res. Lett., 45, 2203–2212, https://doi.org/10.1002/2017GL076644, 2018.
Cheema, M. J. M. and Bastiaanssen, W. G. M.: Land use and land cover classification in the irrigated Indus Basin using growth phenology information from satellite data to support water management analysis, Agr. Water Manage., 97, 1541–1552, https://doi.org/10.1016/j.agwat.2010.05.009, 2010.
Cheng, M. and Ries, J.: The unexpected signal in GRACE estimates of C_{20}, J. Geod., 91, 897–914, https://doi.org/10.1007/s0019001609955, 2017.
Dahle, C., Flechtner, F., Murböck, M., Michalak, G., Neumayer, H., Abrykosov, O., Reinhold, A., and König, R.: GRACE Geopotential GSM Coefficients GFZ RL06 (6.0), GRACE [data set], https://doi.org/10.5880/GFZ.GRACE_06_GSM, 2018.
Dimri, A. P.: Decoding the Karakoram Anomaly, Sci. Total Environ., 788, 147864, https://doi.org/10.1016/j.scitotenv.2021.147864, 2021.
Dimri, A. P., Kumar, D., Chopra, S., and Choudhary, A.: Indus River Basin: Future climate and water budget, Int. J. Climatol., 39, 395–406, https://doi.org/10.1002/joc.5816, 2019.
Ditmar, P.: How to quantify the accuracy of mass anomaly timeseries based on GRACE data in the absence of knowledge about true signal?, J. Geod., 96, 54, https://doi.org/10.1007/s0019002201640x, 2022.
Döll, P., Kaspar, F., and Lehner, B.: A global hydrological model for deriving water availability indicators: model tuning and validation, J. Hydrol., 270, 105–134, https://doi.org/10.1016/S00221694(02)002834, 2003.
Farinotti, D., Immerzeel, W. W., de Kok, R. J., Quincey, D. J., and Dehecq, A.: Manifestations and mechanisms of the Karakoram glacier Anomaly, Nat. Geosci., 13, 8–16, https://doi.org/10.1038/s4156101905135, 2020.
Flechtner, F., Neumayer, K.H., Dahle, C., Dobslaw, H., Fagiolini, E., Raimondo, J.C., and Güntner, A.: What Can be Expected from the GRACEFO Laser Ranging Interferometer for Earth Science Applications?, Surv. Geophys., 37, 453–470, https://doi.org/10.1007/s107120159338y, 2016.
Frenken, K.: Irrigation in Southern and Eastern Asia in Figures: Aquastat survey, 2011, Food and Agriculture of the United Nations, FAO, Rome, 487 p., ISBN 9789251072820, 2012.
Garg, S., Motagh, M., Indu, J., and Karanam, V.: Tracking hidden crisis in India's capital from space: implications of unsustainable groundwater use, Sci. Rep., 12, 651, https://doi.org/10.1038/s41598021041939, 2022.
GRACE Tellus: Monthly Mass Grids – Land  Get Data – GRACE Tellus, https://grace.jpl.nasa.gov/data/getdata/monthlymassgridsland, last access: 17 October 2019.
Groh, A., Horwath, M., Horvath, A., Meister, R., Sørensen, L. S., Barletta, V. R., Forsberg, R., Wouters, B., Ditmar, P., Ran, J., Klees, R., Su, X., Shang, K., Guo, J., Shum, C. K., Schrama, E., and Shepherd, A.: Evaluating GRACE Mass Change Time Series for the Antarctic and Greenland Ice Sheet – Methods and Results, Geosciences, 9, 415, https://doi.org/10.3390/geosciences9100415, 2019.
Harris, I., Jones, P. D., Osborn, T. J., and Lister, D. H.: Updated highresolution grids of monthly climatic observations – the CRU TS3.10 Dataset: Updated HighResolution Grids Of Monthly Climatic Observations, Int. J. Climatol., 34, 623–642, https://doi.org/10.1002/joc.3711, 2014.
Hsu, C.W. and Velicogna, I.: Detection of sea level fingerprints derived from GRACE gravity data: Detecting Sea Level Fingerprints, Geophys. Res. Lett., 44, 8953–8961, https://doi.org/10.1002/2017GL074070, 2017.
Huang, Z., Jiao, J. J., Luo, X., Pan, Y., and Zhang, C.: Sensitivity Analysis of Leakage Correction of GRACE Data in Southwest China Using APriori Model Simulations: InterComparison of Spherical Harmonics, Mass Concentration and In Situ Observations, Sensors, 19, 3149, https://doi.org/10.3390/s19143149, 2019.
Hussain, D., Kao, H.C., Khan, A. A., Lan, W.H., Imani, M., Lee, C.M., and Kuo, C.Y.: Spatial and Temporal Variations of Terrestrial Water Storage in Upper Indus Basin Using GRACE and Altimetry Data, IEEE Access., 8, 65327–65339, https://doi.org/10.1109/ACCESS.2020.2984794, 2020.
ICIMOD: Major River Basins in the Hindu Kush Himalaya (HKH) Region, ICIMOD, 22 April 2021, https://doi.org/10.26066/RDS.2732, 2021.
Jiang, D., Wang, J., Huang, Y., Zhou, K., Ding, X., and Fu, J.: The Review of GRACE Data Applications in Terrestrial Hydrology Monitoring, Adv. Meteorol., 2014, 1–9, https://doi.org/10.1155/2014/725131, 2014.
Klees, R., Revtova, E. A., Gunter, B. C., Ditmar, P., Oudman, E., Winsemius, H. C., and Savenije, H. H. G.: The design of an optimal filter for monthly GRACE gravity models, Geophys. J. Int., 175, 417–432, https://doi.org/10.1111/j.1365246X.2008.03922.x, 2008.
Krakauer, N. Y.: Yearahead predictability of South Asian Summer Monsoon precipitation, Environ. Res. Lett., 14, 044006, https://doi.org/10.1088/17489326/ab006a, 2019.
Kvas, A., Behzadpour, S., Ellmer, M., Klinger, B., Strasser, S., Zehentner, N., and MayerGürr, T.: ITSGGrace2018: Overview and Evaluation of a New GRACEOnly Gravity Field Time Series, J. Geophys. Res.Solid Earth, 124, 9332–9344, https://doi.org/10.1029/2019JB017415, data available at https://www.tugraz.at/institute/ifg/downloads/gravityfieldmodels/itsggrace2018 (last access: 1 April 2022), 2019.
Landerer, F. W. and Swenson, S. C.: Accuracy of scaled GRACE terrestrial water storage estimates: ACCURACY OF GRACETWS, Water Resour. Res., 48, W04531, https://doi.org/10.1029/2011WR011453, 2012.
Lehner, B., Liermann, C. R., Revenga, C., Vörösmarty, C., Fekete, B., Crouzet, P., Döll, P., Endejan, M., Frenken, K., Magome, J., Nilsson, C., Robertson, J. C., Rödel, R., Sindorf, N., and Wisser, D.: Highresolution mapping of the world's reservoirs and dams for sustainable riverflow management, Front. Ecol. Environ., 9, 494–502, https://doi.org/10.1890/100125, 2011.
Long, D., Longuevergne, L., and Scanlon, B. R.: Global analysis of approaches for deriving total water storage changes from GRACE satellites, Water Resour. Res., 51, 2574–2594, https://doi.org/10.1002/2014WR016853, 2015.
Marzeion, B., Jarosch, A. H., and Hofer, M.: Past and future sealevel change from the surface mass balance of glaciers, The Cryosphere, 6, 1295–1322, https://doi.org/10.5194/tc612952012, 2012.
Müller Schmied, H., Adam, L., Eisner, S., Fink, G., Flörke, M., Kim, H., Oki, T., Portmann, F. T., Reinecke, R., Riedel, C., Song, Q., Zhang, J., and Döll, P.: Variations of global and continental water balance components as impacted by climate forcing uncertainty and human water use, Hydrol. Earth Syst. Sci., 20, 2877–2898, https://doi.org/10.5194/hess2028772016, 2016.
NASA JPL (Jet Propulsion Laboratory): GRACE Static Field Geopotential Coefficients JPL Release 6.0, GRACE [data set], https://doi.org/10.5067/GRGSM20J06, 2018.
Peltier, W. R., Argus, D. F., and Drummond, R.: Comment on “An Assessment of the ICE6G_C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al.: The ICE6G_C (VM5a) GIA model, J. Geophys. Res.Solid Earth, 123, 2019–2028, https://doi.org/10.1002/2016JB013844, 2018.
Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J.O., Hock, R., Kaser, G., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radić, V., Rastner, P., Raup, B. H., Rich, J., Sharp, M. J., and The Randolph Consortium: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–552, https://doi.org/10.3189/2014JoG13J176, 2014.
Pradhan, G. P.: Understanding interannual groundwater variability in North India using GRACE, Master's thesis, University of Twente, https://purl.utwente.nl/essays/84053 (last access: 1 July 2020), 2014.
Rodell, M., Velicogna, I., and Famiglietti, J. S.: Satellitebased estimates of groundwater depletion in India, Nature, 460, 999–1002, https://doi.org/10.1038/nature08238, 2009.
Sasgen, I., MartínEspañol, A., Horvath, A., Klemann, V., Petrie, E. J., Wouters, B., Horwath, M., Pail, R., Bamber, J. L., Clarke, P. J., Konrad, H., Wilson, T., and Drinkwater, M. R.: Altimetry, gravimetry, GPS and viscoelastic modeling data for the joint inversion for glacial isostatic adjustment in Antarctica (ESA STSE Project REGINA), Earth Syst. Sci. Data, 10, 493–523, https://doi.org/10.5194/essd104932018, 2018.
Save, H., Bettadpur, S., and Tapley, B. D.: Highresolution CSR GRACE RL05 mascons, J. Geophys. Res.Solid Earth, 121, 7547–7569, https://doi.org/10.1002/2016JB013007, data available at http://www2.csr.utexas.edu/grace (last access: 1 July 2021), 2016.
Scanlon, B. R., Zhang, Z., Save, H., Sun, A. Y., Müller Schmied, H., van Beek, L. P. H., Wiese, D. N., Wada, Y., Long, D., Reedy, R. C., Longuevergne, L., Döll, P., and Bierkens, M. F. P.: Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data, Proc. Natl. Acad. Sci. USA, 115, E1080–E1089, https://doi.org/10.1073/pnas.1704665115, 2018.
Scargle, J. D.: Studies in astronomical time series analysis. II – Statistical aspects of spectral analysis of unevenly spaced data, ApJ, 263, 835, https://doi.org/10.1086/160554, 1982.
Schneider, U., Becker, A., Finger, P., MeyerChristoffer, A., Rudolf, B., and Ziese, M.: GPCC full data reanalysis version 7.0 at 0.50^{∘}: monthly landsurface precipitation from raingauges built on GTSbased and historic data, https://doi.org/10.5676/dwd_gpcc/fd_m_v7_050, 2015.
Shrestha, A. B., Wagle, N., and Rajbhandari, R.: A Review on the Projected Changes in Climate Over the Indus Basin, in: Indus River Basin, Water Sec. Sustain., 2019, 145–158, https://doi.org/10.1016/B9780128127827.000072, 2019.
Sun, Y., Riva, R., and Ditmar, P.: Optimizing estimates of annual variations and trends in geocenter motion and J_{2} from a combination of GRACE data and geophysical models, J. Geophys. Res.Solid Earth, 121, 8352–8370, https://doi.org/10.1002/2016JB013073, 2016.
Swenson, S. and Wahr, J.: Postprocessing removal of correlated errors in GRACE data, Geophys. Res. Lett., 33, L08402, https://doi.org/10.1029/2005GL025285, 2006.
Swenson, S., Chambers, D., and Wahr, J.: Estimating geocenter variations from a combination of GRACE and ocean model output: Estimating Geocenter Variations, J. Geophys. Res., 113, B08410, https://doi.org/10.1029/2007JB005338, 2008.
Tapley, B. D., Watkins, M. M., Flechtner, F., Reigber, C., Bettadpur, S., Rodell, M., Sasgen, I., Famiglietti, J. S., Landerer, F. W., Chambers, D. P., Reager, J. T., Gardner, A. S., Save, H., Ivins, E. R., Swenson, S. C., Boening, C., Dahle, C., Wiese, D. N., Dobslaw, H., Tamisiea, M. E., and Velicogna, I.: Contributions of GRACE to understanding climate change, Nat. Clim. Chang., 9, 358–369, https://doi.org/10.1038/s4155801904562, 2019.
ul Hasson, S.: Future Water Availability from HindukushKarakoramHimalaya upper Indus Basin under Conflicting Climate Change Scenarios, Climate, 4, 40, https://doi.org/10.3390/cli4030040, 2016.
UTCSR (University Of Texas Center For Space Research): Grace Static Field Geopotential Coefficients Csr Release 6.0, GRACE [data set], https://doi.org/10.5067/GRGSM20C06, 2018.
Velicogna, I. and Wahr, J.: Timevariable gravity observations of ice sheet mass balance: Precision and limitations of the GRACE satellite data: Understanding Grace Ice Mass Estimates, Geophys. Res. Lett., 40, 3055–3063, https://doi.org/10.1002/grl.50527, 2013.
Vishwakarma, B. D., Horwath, M., Devaraju, B., Groh, A., and Sneeuw, N.: A DataDriven Approach for Repairing the Hydrological Catchment Signal Damage Due to Filtering of GRACE Products: Repairing Signal Damage Due To Filtering, Water Resour. Res., 53, 9824–9844, https://doi.org/10.1002/2017WR021150, 2017.
Wahr, J., Molenaar, M., and Bryan, F.: Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE, J. Geophys. Res., 103, 30205–30229, https://doi.org/10.1029/98JB02844, 1998.
Willen, M. O., M. Horwath, A. Groh, V. Helm, B. Uebbing, and J. Kusche. : Feasibility of a global inversion for spatially resolved glacial isostatic adjustment and ice sheet mass changes proven in simulation experiments, J. Geod., https://doi.org/10.1007/s00190022016518, accepted, 11 August 2022.
Zhang, L., Dobslaw, H., and Thomas, M.: Globally gridded terrestrial water storage variations from GRACE satellite gravimetry for hydrometeorological applications, Geophys. J. Int., 206, 368–378, https://doi.org/10.1093/gji/ggw153, 2016.