Articles | Volume 26, issue 17
Research article
12 Sep 2022
Research article |  | 12 Sep 2022

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, and 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.

1 Introduction

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 km2). The scaling factor approach has been the most widely used for leakage correction. Its simplicity in application to the gridded TWS products (10×10) 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.

Figure 1The Indus Basin. The blue squares represent the glaciated areas as cells of 1 size. Basin boundaries by the International Centre for Integrated Mountain Development (ICIMOD). Main background map © Google Maps and inset map © Carto.

2 Materials and methods

2.1 Study area

The Indus River basin, shown in Fig. 1 (basin boundaries by ICIMOD, 2021), covers an area of 1.14 million km2. The basin spans over four nations, India, Afghanistan, Pakistan, and China, and supports over 215 million people with an approximate water availability of 1329 m3 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).

Figure 2The land use land cover map (100 m resolution) of the Indus Basin for 2016, obtained from the Copernicus Global Land Cover Service. Background map © CARTO.

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 C20 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.

Table 1Summary of WGHM 2.2d versions used in this study. The integrated WGHM version is standard WGHM with an additional glacier module from GGM. An average of four variants in each version were used in the study.

Download Print Version | Download XLSX

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 km3 (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),

(1) V r , θ , λ , t = GM R e n = 0 N R e r n + 1 m = 0 n P n m cos θ C n m t cos m λ + S n m t sin m λ ,

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, Pnm are the associated Legendre's functions, Cnmt and Snmt are the spherical harmonic coefficients at month t, and Re (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),

(2) Δ σ ( θ λ t ) = R e ρ avg 3 n = 0 N 2 n + 1 1 + k n m = 0 n P n m ( cos θ ) Δ C n m t cos m λ + Δ S n m t sin m λ ,

where ρavg is the average density of Earth (5515 kg m−3) and kn 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),

(3) Δ σ ( θ λ t ) = R e ρ avg 3 n = 0 N 2 n + 1 1 + k n m = 0 n P n m ( cos θ ) W n Δ C n m t cos m λ + Δ S n m t sin m λ ,

where Wn 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),

(4) RF ( θ λ ) = 1 0 ( θ , λ ) R ( θ , λ ) Ω - R ,

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),

(5) Δ M ( θ λ t ) = R e ρ avg 3 n = 0 N 2 n + 1 1 + k n m = 0 n P n m cos θ W n R n m c Δ C n m t cos m λ + R n m s Δ S n m t sin m λ ,

where Rc and Rs 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),

(6) Δ M region t = Δ σ θ , λ , t RF θ , λ d Ω ,

where dΩ=Re2sinθdθdλ 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),

(7) L = t = 1 T Δ S t u - k Δ S t f 2 ,

over the entire period of study. In Eq. (7), ΔStu is the unfiltered and ΔStf is the filtered model mass for the tth month. A scaling factor for each grid cell in a 10×10 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),

(8) Δ S = Δ S long term + Δ S sesonal + Δ S residual .

For this purpose, a time series model is fit to the data using least squares as in Eq. (9),

(9) y i = β 0 + β 1 t + k = 1 n β k c cos 2 π T k t + β k s sin 2 π T k t + ν t ,

where β1 is the linear trend, βkc and βks are the amplitudes of the cosine and sine terms of the periodic component with a period Tk, and t represents the month with respect to the mean of the observation period. However, the periods (Tk) 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 (10-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),

(10) L = j = 1 N Δ S j u - k j Δ S j f ,  where  j =  trend   annual   semi-annual  ,

where kj is the scaling factor for jth component.

Figure 3The schematic diagram for deriving and applying the three scaling schemes.


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.

(11)Initial Leakage Error=stdΔStu-ΔStfc×RMSΔGtRMSΔStf(12)Residual Leakage Error=stdΔStu-kΔStfc×RMSΔGtRMSΔStf,

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, ΔGt represents GRACE mass time series, and ΔStfc 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 C20. 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 C20 (<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.

Figure 4Time series of mass anomalies (mean removed: 2002–2016) in the Indus Basin derived from GRACE spherical harmonic solution (mean of JPL, CSR, and GFZ), corrected with ICE6G, and filtered with Swenson destriping and 300 km Gaussian filter.


Table 2Summary of parameters from time series fitting of basin averages.

Download Print Version | Download XLSX

3 Results and discussion

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 -7.6±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 (R2=0.95) indicates that the SH solution and processing strategy used are nearly as good as the MASCON solution for the Indus Basin. An R2 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 5Basin-averaged time series from GRACE SH, standard WGHM, and integrated WGHM for the Indus Basin.


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 (-20.5±0.4 Gt yr−1) than in the standard version (-10.3±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).

Figure 6Lomb–Scargle periodograms for basin-averaged time series from (a) GRACE SH, (b) unfiltered standard WGHM, (c) unfiltered integrated WGHM, (d) filtered standard WGHM, and (e) filtered integrated WGHM. The peaks marked were selected as dominant for time series fit.


Figure 7Spatial patterns of trend, annual, and semi-annual mass changes from WGHM standard (left) and integrated (right) versions.

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).

Figure 8Mean mass anomalies (in Gt) for each month from April 2002–December 2016 obtained from WGHM Integrated: unfiltered (top) and filtered (bottom). Refer to Fig. 2 for land cover information.

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.

Figure 9Unfiltered and filtered basin-averaged time series from (a) standard WGHM and (b) integrated WGHM version. Note the different y axis scales.


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.

Table 3Basin-scale factors and frequency-dependent scale factors for the Indus Basin. Basin-scale factors from two different studies are given for reference.

Download Print Version | Download XLSX

3.3 Scaling factors

3.3.1 Basin-scale factors

Basin scaling factors of ks=1.14 from the standard version and kg=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.

Figure 10Gridded scaling factor map from the standard (a) and integrated (b) version of WGHM. The grid size is 1 equiangular.

Figure 11Frequency distributions of grid scaling factors from (a) standard and (b) integrated WGHM versions for the Indus Basin. The mean, standard deviation, and coefficient of variation (CV) are inset. Negative scaling factors are excluded.


Table 4Interpretation of gridded scaling factors.

Download Print Version | Download XLSX

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 12Frequency distributions of grid scaling factors from (a) standard and (b) integrated WGHM versions for the non-glaciated region of the Indus Basin. Notice the increased frequency of scaling factors <0.5 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.

Table 5Parameters along with their formal uncertainties from time series fit to basin averages from scaled GRACE SH estimates under different scaling schemes. Also shown are corresponding parameters of DDC-corrected time series.

Download Print Version | Download XLSX

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.

Figure 13Spatial distributions of the components mass anomalies from (a) CSR MASCONS, (b) unscaled SH solutions, (c) grid-scaled SH solutions using standard, and (d) grid-scaled SH solutions using integrated WGHM (left to right panels) in the Indus Basin. The red bounding box indicates the extent of the Karakoram region.

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 14Basin-averaged GRACE time series for Indus Basin with corrected leakage using DDC, grid scaling, and frequency-dependent scaling using integrated WGHM.


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.

Figure 15Comparison of basin-averaged time series from grid-scaled GRACE SH using standard, and integrated WGHM with basin-averaged time series from CSR-M.


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 R2 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 R2 for the de-trended time series shows that the correlation is not spurious.

Figure 16Comparison of basin-averaged time series from frequency-dependent scaling of GRACE SH using standard and integrated WGHM with basin-averaged time series from CSR-M.


Similarly, the result of frequency-dependent scaling is compared with the MASCON time series in Fig. 16. The obtained R2 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.

Table 6Time series components of basin averages from the simulation using standard WGHM.

Download Print Version | Download XLSX

Table 7Time series components of basin averages from the simulation using integrated WGHM.

Download Print Version | Download XLSX

Figure 17True (a) and recovered (c) spatial patterns of different mass change components from grid scaling of filtered–corrupted standard WGHM fields (b) in the simulation experiment.

Figure 18True (a) and recovered (c) spatial patterns of different mass change components from grid scaling of filtered–corrupted integrated WGHM fields (b) in the simulation experiment.

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).

Table 8Noise and initial leakage error in unscaled GRACE basin-averaged time series of Indus Basin.

Download Print Version | Download XLSX

Table 9Scaled noise and residual leakage in scaled GRACE basin-averaged time series of Indus Basin.

Download Print Version | Download XLSX

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.

4 Summary and conclusion

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.

Code availability

The MATLAB scripts written to conduct this study can be obtained upon request to the corresponding author.

Data availability

All GRACE Level 2 release 6 spherical harmonic datasets used in this study are available at (, NASA JPL, 2018;, Dahle et al., 2018;, UTCSR, 2018), supported by the NASA MEaSUREs Program. CSR release 6 MASCONS were downloaded from (Save et al., 2016). The ITSG-2018 normal equations for deriving GRACE errors are available at (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:

Author contributions

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.

Competing interests

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.

Financial support

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).

Review statement

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,, 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,, 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,, 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,, 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,, 2010. 

Cheng, M. and Ries, J.: The unexpected signal in GRACE estimates of C20, J. Geod., 91, 897–914,, 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],, 2018. 

Dimri, A. P.: Decoding the Karakoram Anomaly, Sci. Total Environ., 788, 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,, 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,, 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,, 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,, 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,, 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,, 2022. 

GRACE Tellus: Monthly Mass Grids – Land | Get Data – GRACE Tellus,, 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,, 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,, 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,, 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,, 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,, 2020. 

ICIMOD: Major River Basins in the Hindu Kush Himalaya (HKH) Region, ICIMOD, 22 April 2021,, 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,, 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,, 2008. 

Krakauer, N. Y.: Year-ahead predictability of South Asian Summer Monsoon precipitation, Environ. Res. Lett., 14, 044006,, 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,, data available at (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,, 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,, 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,, 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,, 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,, 2016. 

NASA JPL (Jet Propulsion Laboratory): GRACE Static Field Geopotential Coefficients JPL Release 6.0, GRACE [data set],, 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,, 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,, 2014. 

Pradhan, G. P.: Understanding interannual groundwater variability in North India using GRACE, Master's thesis, University of Twente, (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,, 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,, 2018. 

Save, H., Bettadpur, S., and Tapley, B. D.: High-resolution CSR GRACE RL05 mascons, J. Geophys. Res.-Solid Earth, 121, 7547–7569,, data available at (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,, 2018. 

Scargle, J. D.: Studies in astronomical time series analysis. II – Statistical aspects of spectral analysis of unevenly spaced data, ApJ, 263, 835,, 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,, 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,, 2019. 

Sun, Y., Riva, R., and Ditmar, P.: Optimizing estimates of annual variations and trends in geocenter motion and J2 from a combination of GRACE data and geophysical models, J. Geophys. Res.-Solid Earth, 121, 8352–8370,, 2016. 

Swenson, S. and Wahr, J.: Post-processing removal of correlated errors in GRACE data, Geophys. Res. Lett., 33, L08402,, 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,, 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,, 2019. 

ul Hasson, S.: Future Water Availability from Hindukush-Karakoram-Himalaya upper Indus Basin under Conflicting Climate Change Scenarios, Climate, 4, 40,, 2016. 

UTCSR (University Of Texas Center For Space Research): Grace Static Field Geopotential Coefficients Csr Release 6.0, GRACE [data set],, 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,, 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,, 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,, 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.,, 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,, 2016. 

Short summary
GRACE/GRACE-FO provided global observations of water storage change since 2002. Scaling is a common approach to compensate for the spatial filtering inherent to the results. However, for complex hydrological basins, the compatibility of scaling with the characteristics of regional hydrology has been rarely assessed. We assess traditional scaling approaches and a new scaling approach for the Indus Basin. Our results will help users with regional focus understand implications of scaling choices.