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 hydro-climatic 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/Follow-On (GRACE/GRACE-FO) monthly gravity field models, distributed as spherical harmonic (SH) coefficients, are the only means of observing this state variable. The well-known 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 frequency-dependent scaling for leakage correction of GRACE TWSC in a unique, basin-specific assessment for the Indus Basin. We harness the characteristics of significant heterogeneity in the Indus Basin due to climate and human-induced 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 GRACE-like errors using full monthly error covariance matrices. We find that frequency-dependent 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. Frequency-dependent scaling can provide a robust estimate of the seasonal cycle of TWSC for practical applications such as regional-scale 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) -
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). Regional-scale 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 point-scale 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, global-scale 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 continental-scale 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 inter-satellite 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 inter-satellite 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 hydro-geodetic communities, and have been traditionally carried
out using three hydrological or land surface model-dependent 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 model-independent 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 timescale-dependent 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 timescale-dependent approach is far less studied. The timescale-dependent (or frequency-dependent) 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 sea-level 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 GRACE-like errors to validate the inferences made in a region-specific inter-comparison 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 climate-induced and human-induced 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 hydro-climatic 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 basin-specific assessment, develop a novel frequency-dependent 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
data-scarce 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
semi-arid to temperate sub-humid 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
hydro-climatic 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 evapo-transpiration (ET) is also high in these regions. Hence, human intervention is the dominant short-scale 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 ICE-6G_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 (TN-13) on the Tellus website, were used. The replacement
*C*_{20} coefficients from (Cheng and Ries, 2017), distributed as TN-11 on
the Tellus website, were used. To derive GRACE-like errors for the design of
simulation experiments, the monthly normal equations provided by TU Graz for
the ITSG-2018 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; ERA-Interim
reanalysis (WFDEI) applied with WATCH Forcing Data (WFD) methodology and
bias-corrected 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
near-surface 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 one-layer soil water storage compartment characterized by a land cover and soil-specific 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 sub-compartments: 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 groundwater-depleted 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. Country-specific 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 co-latitude 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 ICE-6G 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 degree-dependent 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 mid-latitude 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 degree-dependent factors of Gaussian filter in
spectral-domain (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 spectral-domain (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
*t*th 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 frequency-dependent scaling factors, the total mass changes from unfiltered and filtered model versions are decomposed into long-term 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 semi-annual 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 semi-annual 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 *j*th 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 data-driven correction approach

The DDC approach or the method of deviation (Vishwakarma et al., 2017) uses twice-filtered 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 basin-averaged 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 GRACE-like 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, ITSG-2018) 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 GRACE-like 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 ITSG-based 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 un-modeled inter-annual
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 18-month 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 (CSR-M) 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 post-processing step. Figure 4 also shows
the GRACE SH and CSR-M time-series 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 de-trended, 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 SH-based 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 human-induced
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 climate-induced 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 semi-annual 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 time-series fit (Table 2), the annual amplitude (11.8±3 Gt) obtained from GRACE is much smaller than the semi-annual 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 semi-annual 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 semi-annual seasonality in TWS changes results from bimodal precipitation distribution in the Indus Basin. The inter-annual signal content (area of the frequency spectrum in the inter-annual bandwidth) is almost identical in both model versions and GRACE. This shows that adding glacier mass changes to WGHM does not contribute to the inter-annual variations. These inter-annual variations are probably the result of long-term groundwater depletion (Pradhan, 2014) and inter-annual 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 small-scale changes in TWS are smoothed out and reduced in amplitude. Large-scale 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 pre-monsoon months (March to May) can be attributed to western disturbances. It is worth mentioning that there is generally a 1-month 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. Small-scale 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 model-based trends closer to the GRACE-based trends, partly explaining the trend differences between filtered GRACE-based and unfiltered model-based 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 semi-annual frequencies, it is also clear that filtering affects annual and semi-annual mass changes differently. The time series fit (Table 2) shows that filtering significantly reduces the annual amplitudes in both versions, while the semi-annual 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 Basin-scale 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 basin-scale 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 basin-scale 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 out-of-phase
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 leakage-in. These small scaling factors occur mainly
in the lower reaches of the basin (Indus plains) and upper Indus Basin,
where non-glacier 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 (north-east 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
non-glaciated region of the Indus Basin (i.e., the non-glaciated 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 Frequency-dependent scaling factors

The frequency-dependent 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 frequency-independent basin-scale factors. This reinforces that the Indus Basin is dominated by long-term trends rather than seasonal signals, which causes the basin-scale 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 basin-scale factor alone would under-scale the annual component. The semi-annual 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 basin-scale factor reinforce the need for frequency-dependent 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 frequency-dependent 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 basin-averaged 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 frequency-dependent
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 semi-annual 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 grid-scaled 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 basin-averaged GRACE time series from grid and frequency-dependent scaling using integrated WGHM compared to the DDC time series (Sect. 2.3.3).

Figure 14 shows that overall, frequency-dependent 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 under-scale certain pixels. Grid scaling can be seen to overestimate the trends compared to DDC since the trend-contributing 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 semi-annual frequencies. This highlights the differences in physical processes governing the large-scale and fine-scale behavior of these components. These results provide confidence in the new frequency-dependent scaling using WGHM to restore the damaged signal at the catchment scale correctly.

A comparison of the grid-scaled GRACE SH time series to the MASCON time
series (CSR-M) 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
de-trended time series shows that the correlation is not spurious.

Similarly, the result of frequency-dependent 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
inter-annual, 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 frequency-dependent scaling as the best performing scheme in terms of recovering the true basin-averaged signal for both the model versions. The similarity of recovered trends from basin and frequency-dependent 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 semi-annual amplitude from grid scaling using integrated version is reduced, while increased using the standard version (compared to the corresponding filtered–corrupted semi-annual component). This supports the behavior seen in scaled GRACE with the following explanation: most of the semi-annual signal contribution comes from non-glaciated 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 semi-annual 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 semi-annual 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 semi-annual 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 frequency-dependent 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 frequency-dependent 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 frequency-dependent 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 frequency-dependent 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 spatio-temporal characteristics: basin scaling factors that are spatially and temporally constant; gridded scaling factors that are spatially variable while temporally constant; and frequency-dependent 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 inter-comparison 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 frequency-dependent 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 DDC-based estimate shows that it can be a viable alternative as a leakage correction approach for the Indus Basin. The frequency-dependent 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 decision-making 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 frequency-dependent 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, grid-scaled 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 non-glaciated 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 grid-scaled GRACE, such as trends and seasonal amplitudes, can be misleading in judging the performance of grid scaling due to over- and under-scaled pixels compensating, and leading to seemingly realistic basin-averaged 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 long-term trends in the basin. However, a better signal–noise separation must be achieved to minimize the scaling of noise. Frequency-dependent 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 frequency-dependent scaling scheme that would provide three scaling factors per grid cell. However, our initial experiments show that the time series decomposition at grid-scale leads to vastly varying annual and semi-annual 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 inter-comparison 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 GRACE-FO 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/GRGSM-20J06, NASA JPL, 2018; https://doi.org/10.5880/GFZ.GRACE_06_GSM, Dahle et al., 2018; https://doi.org/10.5067/GRGSM-20C06, 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 ITSG-2018 normal equations for deriving GRACE errors are available at https://www.tugraz.at/institute/ifg/downloads/gravity-field-models/itsg-grace2018 (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/hess-26-4515-2022-supplement.

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 3-D 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/hess-24-4831-2020, 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/s00190-016-0995-5, 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 time-series based on GRACE data in the absence of knowledge about true signal?, J. Geod., 96, 54, https://doi.org/10.1007/s00190-022-01640-x, 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/S0022-1694(02)00283-4, 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/s41561-019-0513-5, 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 GRACE-FO Laser Ranging Interferometer for Earth Science Applications?, Surv. Geophys., 37, 453–470, https://doi.org/10.1007/s10712-015-9338-y, 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/s41598-021-04193-9, 2022.

GRACE Tellus: Monthly Mass Grids – Land | Get Data – GRACE Tellus, https://grace.jpl.nasa.gov/data/get-data/monthly-mass-grids-land, 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 high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset: Updated High-Resolution 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 A-Priori Model Simulations: Inter-Comparison 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.1365-246X.2008.03922.x, 2008.

Krakauer, N. Y.: Year-ahead predictability of South Asian Summer Monsoon precipitation, Environ. Res. Lett., 14, 044006, https://doi.org/10.1088/1748-9326/ab006a, 2019.

Kvas, A., Behzadpour, S., Ellmer, M., Klinger, B., Strasser, S., Zehentner, N., and Mayer-Gürr, T.: ITSG-Grace2018: Overview and Evaluation of a New GRACE-Only 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/gravity-field-models/itsg-grace2018 (last access: 1 April 2022), 2019.

Landerer, F. W. and Swenson, S. C.: Accuracy of scaled GRACE terrestrial water storage estimates: ACCURACY OF GRACE-TWS, 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.: High-resolution mapping of the world's reservoirs and dams for sustainable river-flow 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 sea-level change from the surface mass balance of glaciers, The Cryosphere, 6, 1295–1322, https://doi.org/10.5194/tc-6-1295-2012, 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/hess-20-2877-2016, 2016.

NASA JPL (Jet Propulsion Laboratory): GRACE Static Field Geopotential Coefficients JPL Release 6.0, GRACE [data set], https://doi.org/10.5067/GRGSM-20J06, 2018.

Peltier, W. R., Argus, D. F., and Drummond, R.: Comment on “An Assessment of the ICE-6G_C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al.: The ICE-6G_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.: Satellite-based estimates of groundwater depletion in India, Nature, 460, 999–1002, https://doi.org/10.1038/nature08238, 2009.

Sasgen, I., Martín-Españ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/essd-10-493-2018, 2018.

Save, H., Bettadpur, S., and Tapley, B. D.: High-resolution 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., Meyer-Christoffer, A., Rudolf, B., and Ziese, M.: GPCC full data reanalysis version 7.0 at 0.50^{∘}: monthly land-surface precipitation from rain-gauges built on GTS-based 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/B978-0-12-812782-7.00007-2, 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.: Post-processing 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/s41558-019-0456-2, 2019.

ul Hasson, S.: Future Water Availability from Hindukush-Karakoram-Himalaya 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/GRGSM-20C06, 2018.

Velicogna, I. and Wahr, J.: Time-variable 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 Data-Driven 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/s00190-022-01651-8, 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.