What is the spatial resolution of GRACE satellite products for hydrology? Remote Sensing

: The mass change information from the Gravity Recovery And Climate Experiment ( GRACE ) satellite mission is available in terms of noisy spherical harmonic coefﬁcients truncated at a maximum degree (band-limited). Therefore, ﬁltering is an inevitable step in post-processing of GRACE ﬁelds to extract meaningful information about mass redistribution in the Earth-system. It is well known from previous studies that a number can be allotted to the spatial resolution of a band-limited spherical harmonic spectrum and also to a ﬁltered ﬁeld. Furthermore, it is now a common practice to correct the ﬁltered GRACE data for signal damage due to ﬁltering (or convolution in the spatial domain). These correction methods resemble deconvolution, and, therefore, the spatial resolution of the corrected GRACE data have to be reconsidered. Therefore, the effective spatial resolution at which we can obtain mass changes from GRACE products is an area of debate. In this contribution, we assess the spatial resolution both theoretically and practically. We conﬁrm that, theoretically, the smallest resolvable catchment is directly related to the band-limit of the spherical harmonic spectrum of the GRACE data. However, due to the approximate nature of the correction schemes and the noise present in GRACE data, practically, the complete band-limited signal cannot be retrieved. In this context, we perform a closed-loop simulation comparing four popular correction schemes over 255 catchments to demarcate the minimum size of the catchment whose signal can be efﬁciently recovered by the correction schemes. We show that the amount of closure error is inversely related to the size of the catchment area. We use this trade-off between the error and the catchment size for deﬁning the potential spatial resolution of the GRACE product obtained from a correction method. The magnitude of the error and hence the spatial resolution are both dependent on the correction scheme. Currently, a catchment of the size ≈ 63,000 km 2 can be resolved at an error level of 2 cm in terms of equivalent water height.


Introduction
The Gravity Recovery And Climate Experiment (GRACE) satellite mission has provided valuable information towards understanding the continental to regional scale hydrology [1,2]. The time-variable gravity field of the Earth from the GRACE satellite mission can be obtained from various organizations, at various levels of complexity and time resolution. Monthly, weekly and even daily fields are provided either in terms of band-limited spherical harmonic coefficients or global grids of mass change from Center for Space Research (CSR), GeoForschungZentrum (GFZ), Jet Propulsion Laboratory (JPL), and several other processing centers. The most commonly used GRACE products are the monthly scales. Rowlands et al. [31] proposed a limit of ≈150,000 km 2 . On the other hand, Lorenz et al. [7] demonstrated that many catchments smaller than these limits were well observed by GRACE provided they have a strong seasonal cycle in terms of water storage changes. Tourian et al. [32] showed that GRACE was able to capture the water loss in Urmia basin, which has an area of ≈52,000 km 2 . In a recent study, Khaki et al. [33] proposed a two-step Kernel Fourier Integration (KeFIn) filter to reduce errors in high-frequency mass changes and also to decrease spatial leakage errors. While authors show that the filter reduces mass estimation errors in many small and medium size river basins, they do not evaluate the spatial resolution of GRACE data.
In this contribution, we discuss the ideal spatial resolution of the GRACE product and how it changes with the post-processing strategy we choose. We use popular methods and tools developed in the past decade by fellow researchers to help us set a limit to the catchment size that can be observed with an accuracy better than a certain limit. We carry out the investigation in a closed-loop simulation environment that emulates GRACE satellite products. In order to be comprehensive in our analysis, we analyze the error behaviour with respect to the catchment size for four popular repairing strategies over 255 catchments. In general, we find that the error increases as the catchment size decreases, but we observe that the amount of error varies from one correction method to the other. We also find that one can obtain better GRACE time-series over smaller catchments with the data-driven correction scheme by Vishwakarma et al. [25].

Spherical Harmonic Coefficients and Their Corresponding Spatial Resolution
In this section, we will revisit the idea of spatial resolution of a field given in the form of spherical harmonic coefficients by recapitulating the state-of-the-art in this domain. We will start with a continuous field f (θ, λ) and its representation in terms of spherical harmonics f (θ, λ) = ∞ ∑ l=0 l ∑ m=0P lm (cos θ) (C lm cos mλ +S lm sin mλ) , 2π λ=0 f (θ, λ)P lm (cos θ) cos mλ sin mλ sin θ dθ dλ, where θ, λ are the co-latitude and longitude of a point in the field f (·) andP lm (cos θ) are the fully normalized associated Legendre functions of degree l and order m normalized [34]. The computation of f (θ, λ) from the spherical harmonic coefficientsC lm ,S lm is called spherical harmonic synthesis (1) and the computation ofC lm ,S lm from the field is called spherical harmonic analysis (2).

Sampling and Half-Wavelength of the Field
Practically, a continuous field is always discretely sampled, where sampling controls the information content that can be ascertained about the field. Spherical harmonic coefficients computed from a discretely sampled dataset are band-limited, i.e., the spectrum is limited to a finite degree L. Thus, Equation (1) becomes In the case of the GRACE time-variable gravity field, data are conventionally disseminated in terms of spherical harmonic coefficients that have been analysed from the GRACE range-rate observables. The maximum degree and order of the spherical harmonic expansion is dictated by the modified Colombo-Nyquist sampling rule for satellite gravimetry [35]. On the other hand, if the spherical harmonic coefficients are computed from a regularly gridded dataset, then the sampling in the latitude direction must be N θ ≥ L + 1 and that in the longitude direction must be N λ ≥ 2L + 1 [36,37]. It is this sampling that gives us the first idea of resolution, the half-wavelength of a field: It should be clear from Equation (4) that the half-wavelength describes the sampling distance and not the resolution itself. This is all the more obvious along the longitude as it clearly follows the Shannon-Nyquist rule in signal processing.
The half-wavelength of a field has been used as the de facto value for describing the resolution of a field. Since it is derived from sampling on an equiangular grid [36], the spatial resolution of a band-limited field is also interpreted in terms of the equiangular grid, especially in terms of its pixel size. A complication in using the pixel area/size as the spatial resolution value is that, due to the convergence of the spherical coordinates at the poles, the pixel size along the latitude circles becomes smaller as we go closer to the poles, and finally vanishes at the poles. Due to this reason, the pixel size cannot be treated as the resolution homogeneously throughout the sphere. Furthermore, if the half-wavelength value is treated as the resolution, it will be a highly optimistic value (cf. Table 1) because it is just the minimum spatial sampling required on an equiangular grid to estimate the spherical harmonic coefficients up to a certain harmonic degree: It is worthwhile to note here that the sampling theorem on the sphere is not unique [38].

Ideal Spatial Resolution
The problem of spatial resolution can be dealt with by treating the band-limit of the spherical harmonic spectrum as a spectral filter, and then by ascertaining the ideal spatial resolution of the filter [17]. Ideal spatial resolution of a filter is defined as the minimum distance required between two signals of equal magnitude, such that they are still seen as two different signals after filtering. The spatial resolution is closely tied to the contrast of the field, i.e., the prominence of the signals after filtering (cf. Figure 1), and, therefore, resolution is not just a number, but it is devised as a function of contrast (modulation) called the modulation transfer function (MTF) [17]. The MTF is a function of modulation transfer and separation between two signals of interest (cf. Figure 1): The MTF helps us identify the minimum spherical distance (ideal resolution) between two Dirac pulses above which they are resolved as two distinct signals. Furthermore, from the slope of the MTF, the level of contrast that will be retained after filtering can also be ascertained-the steeper the slope, the more the contrast. Essentially, the ideal spatial resolution tells us that below this value any two signals of interest are indistinguishable. The MATLAB scripts (Mathworks, Natick, MA, USA) to compute MTF can be downloaded from http://gracebundle.tuxfamily.org.
The band-limited field in Equation (3) can be rewritten as a box-car filtered field as follows (see Appendix A for details): The degree-dependent spectrum B l when propagated to the spatial domain becomes the well-known Shannon kernel [39]. The MTF of the Shannon kernel for the typical GRACE bandwidths provides their ideal spatial resolution (cf. Figure 2), and, obviously, those numbers are higher than the half-wavelengths (sampling distance) that are usually used to designate their resolution (cf. Table 1). Since the Shannon kernel is a homogeneous and isotropic kernel (see Appendix A for definitions of homogeneous and isotropic), its ideal spatial resolution is also homogeneous and isotropic throughout the sphere [14]. Thus, with the ideal spatial resolution values, we are able to define the nominal size of features that can be clearly distinguished without ambiguity, and that they are valid throughout the sphere.  Figure 1. Illustration of the ideal spatial resolution computation for filter windows. Two Dirac pulses (gray lines in (a-e)) are set up at points P and Q and the field is filtered (black lines in (a-e)). For a given filter, they are resolved as two different signals at a certain separation (c)-the ideal spatial resolution. However, the signals are more prominent after filtering (contrast), if the distance between the signals is larger (d,e). The contrast is quantified via the modulation transfer function (MTF), which is a function of the modulation transfer (cf. (6)) and signal separation (f). Due to their high noise levels, the GRACE monthly solutions are filtered prior to their use. Filtering suppresses noise, but, at the same time, has a damaging effect on the signal-it attenuates the signal, introduces leakage and changes the resolution. A Gaussian filter of half-width radius 400 km (3.6 • ) has an ideal resolution of 680 km (6.1 • ), which is far higher than the innate resolution of nearly all the typical band-limited GRACE fields (cf. Table 1). Due to signal attenuation, filtering also reduces signal contrast, which can be seen from the unfavourably gentle slope of the Gaussian MTF curve (cf. Figure 2). Table 1. Comparison of the half-wavelength and ideal spatial resolution values for typical GRACE bandwidths. Also shown are the typical area associated with them. The area of an equi-angular grid cell centered on the equator (A g ) is associated with the half-wavelength, and the spherical area (A s ) is associated with the ideal spatial resolution. N/A stands for not applicable.

L
Half-Wavelength Ideal Spatial Resolution

Catchment Averages, Post-Filtering Corrections and Resolution
Thus far, the spatial resolution has been treated as a spherical distance, but it has to be translated into an areal quantity. As resolution of a band-limited field is isotropic, the area of the smallest resolvable feature will be a spherical cap with the diameter of the ideal spatial resolution. The spherical cap areas for typical GRACE bandwidths are given in Table 1.
The gridded GRACE products can be synthesized at a grid spacing smaller than the ideal spatial resolution. Therefore, it is logical to study the hydrological signal not at point scale, but rather at catchment scale [5,7,21]. It is common practice to compute a catchment averagef c from the filtered fieldf (θ, λ), which is written as Since area aggregation is an operation performed over the spatial coordinates, the spherical harmonic spectrum is not tampered with. Therefore, it can be said with confidence that the minimum resolvable area does not change with the area aggregation operation.
The catchment average is computed for every epoch to obtain a time-series demonstrating the hydrological behaviour of the catchment. Since filtering damages the signal, the catchment scale products from filtered GRACE fields are corrected [5,8,[19][20][21]24,25]. Filtering can be written as a convolution integral [11,14,21,40], and, loosely speaking, the repair scheme imitates a deconvolution process. Therefore, if the repair method does a good job, the resolution must improve. However, it is not a perfect process and that means we will not reach the ideal spatial resolution of the bandwidth, but it will definitely be better than that of the filtered field.

Some Exceptions
As mentioned earlier, satellite gravimetry does not image the field like optical/microwave remote sensing satellites, but provides it via geophysical inversion. Therefore, the values of spatial resolution are intricately tied to the inversion scheme. In the context of band-limited spherical harmonic representation of the GRACE observations, some key aspects have to be kept in mind while using or interpreting the values in Table 1 for catchment aggregated values. Small catchments 1.
that have a strong seasonal variation in their water storage [7], or 2.
that are similar in magnitude and temporal phase as their neighboring catchments, i.e., without spatial contrast [23] (it should not be confused with the signal contrast/modulation mentioned earlier), or 3.
that are isolated or dominant in terms of their signal strength, for example huge reservoir volume changes [41] can still be seen despite the resolution limitations of the given bandwidth. For example, a small but strong signal source will be detected by the GRACE satellites, but, due to the band-limited nature of the spherical harmonic representation, the signal is spread over the area equivalent to the corresponding ideal resolution. The area aggregated value over such a catchment will contain a fair part of the signal from the catchment, and, depending on the above-mentioned conditions, it might as well be mixed with other sources in the vicinity. Therefore, such signal sources-catchments-can be studied with GRACE, but a priori information and signal segregation methods are needed to extract their full signal strength. This might be the reason that, for some very small catchments, Lorenz et al. [7] observe good correlation between the GRACE data and the mass change rate from the water balance equation, but with strong bias and significant differences in the amplitude.

Data and Method
In order to use the relation between the accuracy and catchment size for identifying the spatial resolution of GRACE products in a comprehensive manner, we analyze more than 250 catchments shown in Figure 3. The largest catchment is the Amazon in South America with an area of ≈4,672,876 km 2 and the smallest is Cavally catchment in West Africa with an area of ≈30,744 km 2 . The study is carried out in a closed-loop simulation environment used by Vishwakarma et al. [21] and Vishwakarma et al. [25]. It consists of total water storage anomaly from the Global Land Data Assimilation System (GLDAS) Noah Land Surface Model as the background truth [42]. These global hydrological fields are contaminated with noise extracted from GFZ GRACE fields [3]. In order to extract noise, we first filter the GFZ GRACE fields with a destriping and a 400 km Gaussian filter [10], then we subtract the filtered fields from the corresponding noisy GRACE fields to derive a realistic GRACE-type noise. The GRACE-type noise is then added to the GLDAS model fields for 72 months to obtain GRACE-type noisy fields with known truth. Further details about the simulated field can be found in Vishwakarma et al. [21]. We can process the GRACE-type fields to obtain time-series for a catchment and it can be validated against the time-series from GLDAS fields. In this setup, we can assess the magnitude of error accurately because we have control over each component. We filter the simulated GRACE-type noisy fields f (θ, λ) with a Gaussian 400 km filter and then compute catchment averagesf c for 255 catchments to obtain the respective time-series. Let us denote the corrected time-series byf c . The true catchment averagef c from the band-limited field is not known, but the regional average from the filtered fieldf c is known. One may approach true catchment average, which is approximately equal to the catchment average from the band-limited field, from catchment average of the filtered fieldf c with the help of correction methods. Out of many available methods, we choose four methods: the Multiplicative method by Longuevergne et al. [19], the scaling method by Landerer et al. [24], the additive method by Klees et al. [20], and the data-driven method by Vishwakarma et al. [25]. The mathematical relation that helps you correct for the signal damage due to filtering is given below: [25].
The multiplicative method approaches corrected time-seriesf c by first removing leakage l m c from catchment aggregate of filtered GRACE fieldf c , and then amplifying it by a scale factor s. Leakage is obtained from a hydrological model and the scale factor s is obtained from catchment mask and filtered catchment mask: where R(θ, λ) is the catchment mask: 1 inside and 0 outside,R(θ, λ) is the filtered catchment mask, Ω is the domain of the surface of the Earth, and dΩ is the infinitesimal surface element sin θdθdλ. The additive method approaches corrected time-seriesf c by subtracting leakage l m c from the catchment aggregate of filtered GRACE fieldf c and adding bias b m c , where both leakage and bias are obtained from a hydrological model. The scaling method approaches corrected time-series by multiplying GRACE products by a scale factor k that is obtained by estimating a multiplicative factor, between a hydrological model and its filtered version, via least squares estimation.
The data-driven method corrects time-series by removing leakage l c and deviation integral δF c from the catchment average of filtered GRACE fieldsf c . To compute the leakage and the deviation integral, we do not rely on hydrological models, but use the GRACE fields only. Here, one may argue that the GRACE fields are noisy, thus the leakage and the deviation integral computed from noisy fields will not be accurate. This is the reason they are computed from filtered GRACE fields. Since we know that the filtered GRACE fields are not equal to the truth, the leakage and the deviation integral from these fields are also not accurate. However, Vishwakarma et al. [25] demonstrated that if we multiply the leakage and the deviation integral from filtered GRACE fields, by a scalar ratio between the leakage and the deviation integral from once filtered GRACE fields and twice filtered GRACE fields, we can approach near-truth leakage and deviation integral. Although this approximation works for most of the catchments, it has limitations over arid regions and the estimated leakage and deviation integral for these regions is less accurate. Nevertheless, the accuracy of data-driven methods is still either at par or better than what we can achieve from other methods [25]. Please refer to Vishwakarma et al. [25] and Vishwakarma et al. [21] for more information. Since the multiplicative, the additive, and the scaling approach use a hydrological model to estimate corrected time-series, we call them model-dependent approaches, while the data-driven approach uses GRACE fields only and do not depend on models, we can also call it model-independent.
Global hydrological models have huge uncertainties that vary in space and in time, which is responsible for their differences with GRACE fields [43]. Using these models for correcting GRACE products brings their uncertainty to the final GRACE product. In order to imitate the impact of using a model for correcting GRACE, we prefer not to use the GLDAS model that is also the background truth; instead, we use the WaterGAP hydrological model (WGHM model [44]) to compute model-dependent correction terms employed by model-dependent methods. The data-driven method computes its correction terms, leakage l c and the deviation integral δF c , from once filtered and twice filtered GRACE-type noisy fields. The MATLAB scripts for the data-driven method can be downloaded from https://www.researchgate.net/publication/324804360_Matlab_scripts_Data_ driven or from Institute of Geodesy, Stuttgart web page http://www.gis.uni-stuttgart.de/research/ projects/DataDrivenCorrection/.
The corrected time-series from the above-mentioned four methods is compared to the truth obtained from GLDAS model fields for 255 catchments. The following statistical measures are computed to characterize the behavior of difference (error) between the corrected time-series and the true time-series with respect to the catchment size: • Root Mean Square of Error (RMSE): • Cyclostationary Nash-Sutcliffe Efficiency (NSE seas ) [45,46]: where f c represents the true value obtained from GLDAS fields,f c is the mean annual behaviour (mean monthly values) and m is the number of epochs. RMSE can attain any positive value, a RMSE close to zero represents excellent agreement between f c andf c . NSE seas can attain any value between −∞ and 1. A positive NSE seas value indicates that the difference between the true signal ( f c ) and the corrected time-series (f c ) is well below the non-seasonal variations in the time-series. Typically, hydrological signals have a clear seasonal signal, i.e., change in the amplitude and/or sign of the signal from the summer months to the winter months. However, every summer/winter is not the same and there is a random variation in the amplitude of the signal from year-to-year, which is the non-seasonal variation (the natural variability of the signal). When the error in corrected time-series is smaller than this natural variability, then we can conclude that the signal is very well restored. If we were to use the NSE as proposed by Nash and Sutcliffe [47], without accounting for the cyclostationary seasonal signal, the differences will be compared to the amplitudes of the annual signal. This will not provide a proper indicator for the efficacy of the repairing schemes.

Results and Discussion
In Figures 4 and 5, we have plotted the RMSE and the NSE seas for four corrected time-series and time-series from filtered fields over 255 catchments. The catchments are sorted by their area. Since the scatter of these statistical measures is wide, we fit a smooth line, obtained by using locally weighted scatter-plot smoothing (LOESS), to identify the general pattern. LOESS is a non-parametric method that uses iterative locally weighted regression to fit a second degree polynomial through points in a scatter plot [48]. Since the number of large catchments is smaller than the number of small catchments, the data points on the x-axis, i.e., catchment area, suffer from unequal intervals. Therefore, we have defined bins to counter this irregular data distribution. The first bin is for catchments from 30,000 km 2 to 100,000 km 2 with sample points every 1000 km 2 , the next bin is for catchments from 100,000 km 2 to 1,000,000 km 2 with sample point every 10,000 km 2 , and the last bin is for catchment from 1,000,000 km 2 to 4,000,000 km 2 with sample points every 100,000 km 2 . A window size of 51 is used for smoothing.
We have chosen these parameters to obtain a smooth fit that would represent the general behaviour of the scatter.
We can observe that the error, in general, increases as the area of the catchment decreases, irrespective of the correction method used. However, the rate at which the error increases varies from one method to the other. The multiplicative approach, due to a larger scale factor for smaller catchment, experiences a steep decay in performance. The additive approach, the scaling approach and the time-series from filtered fields are competitive. However, the data-driven method is able to provide better mass change estimates at all catchment scales.
Although, in general, the disagreement between the corrected and the true time-series increases as the catchment size decreases, several small catchments exhibit less error in comparison to a few relatively large catchments. The magnitude of error corresponding to a catchment is pertinent to this simulation setup, and it will change a little bit if we change the background models and simulation setup. Nevertheless, the general pattern that error increases as catchment size decreases will hold.
Ideally, with no noise and no approximations, we expect to recover the full band-limited signal after applying the data-driven method (see Figure 3, Vishwakarma et al. [25]). However, in reality, we lose some accuracy and the error is not zero even for the largest catchment Amazon. The error in the corrected time-series, for catchments smaller than the resolution of the filtered field and larger than the ideal resolution of band-limited GRACE data, should be close to the error in time-series from filtered fields for larger catchments. We can use the trade-off between accuracy and the catchment size to discuss the potential GRACE resolution: those catchments with an error less than a defined threshold can be categorized as resolvable.
Let us say we can accept the GRACE products for all the catchments with an RMSE better than a certain value, say 1 cm or 2 cm. We leave it to the user to define the maximum tolerable error for their application and then decide whether the catchment or region of interest is suitable for analysis with GRACE products corrected with a certain repair scheme. For example, if we choose 1 cm as the error limit, then filtered products can be used to monitor catchments larger than ≈750,000 km 2 , a multiplicative approach can be used for catchments larger than ≈1,563,000 km 2 , an additive approach for catchments larger than ≈1,103,000 km 2 , a scaling approach for catchments larger than ≈563,000 km 2 and the data-driven method for catchments larger than ≈250,000 km 2 . If the tolerable error is 2 cm, then the limit is ≈152,000 km 2 for filtered and additive, ≈810,000 km 2 for multiplicative and ≈63,000 km 2 for both scaling and the data-driven method. In Table 2, we have summarized the resolvable catchment corresponding to an error level. Table 2. The approximate resolvable catchment size (in 1000 km 2 ) for a repair method performing better than a given RMSE. The values given here are obtained from the fit, and one should be more careful while carrying out studies for catchments close to the limit as the scatter is wide. Our suggestion is to avoid any interpretation after the fit tends to become noisy. When the fit for a method becomes noisy, we represent it by writing N/A.

RMSE (cm) Filtered Multiplicative Additive
Scaling Data-Driven Since NSE seas accounts for both the difference in the magnitude and the correlation between the corrected and the true time-series, it is a powerful indicator and popularly used in time-series analysis. One may choose an acceptable NSE seas value for their application to determine the corresponding spatial resolution. For example, if we choose 0.8 as an acceptable limit, then filtered products can be used to monitor catchments larger than ≈1,100,000 km 2 , a multiplicative approach can be used for catchments larger than ≈1,750,500 km 2 , an additive approach for catchments larger than ≈1,100,500 km 2 , a scaling approach for catchments larger than ≈600,000 km 2 and the data-driven method for catchments larger than ≈200,000 km 2 . In Table 3, we have provided the approximate spatial resolution corresponding to a NSE seas value. It is to be noted that the LOESS fit for the multiplicative approach and for the additive approach never attain NSE seas value of 0.9. Furthermore, the spatial resolution improves with increasing NSE seas values, while the corresponding RMSE decreases. Table 3. The approximate resolvable catchment size (in 1000 km 2 ) for a repair method performing better than a given NSE seas . The values given here are obtained from the fit, and one should be more careful while carrying out studies for catchments close to the limit as the scatter is wide. Our suggestion is to avoid any interpretation after the fit tends to become noisy. When the fit for a method becomes noisy, we represent it by writing N/A.   Tables 2 and 3 are only meaningful if the LOESS fit is smooth; as soon as the fit seems noisy, we should be careful in interpreting the resolvable catchment size. Therefore, in this analysis, we have ignored all the catchments for which the fit tends to be noisy-for example, the fit for the data-driven method becomes noisy for catchments below ≈600,000 km 2 . We can find a lower resolvable catchment area for the data-driven method corresponding to an RMSE level of 2.5 cm, but this part of the fit is noisy and we avoid any interpretation. Please note that, within these suggested spatial resolution limits, where one can choose from many methods, one method may perform better than the others as per the fit. Such a general rule can be used but with the understanding that individual catchments may deviate from the defined error behaviour. This is shown by the disagreement between the scatter and the fit.

Conclusions
The GRACE satellite mission has helped us observe continental scale hydrology. Moreover, with improved data processing skills, we have been able to use GRACE products for monitoring catchment scale hydrology. Although the spatial resolution of GRACE is accepted to be coarse, it was necessary to put a number to the spatial resolution of the band-limited GRACE fields and of the filtered GRACE fields. It was established by previous contributions that ideal spatial resolution can be used to define their spatial resolutions.
Notwithstanding this, there have been many attempts to identify the potential spatial resolution of the GRACE products in terms of catchment size, but every attempt concluded with a different number. The difference in the research findings can be attributed to usage of different processing methods and selective regional analysis. For example, one can choose a filter out of many available ones and then a corrective method for repairing the signal damage due to the filtering. Thus, the end product is influenced by combined effects of different processes and it is hard to quantify whether the final product is able to resolve a catchment efficiently. With developments in post-processing algorithms for GRACE data, it was necessary to understand the contemporary spatial resolving power of GRACE products. This study provides a comprehensive analysis of GRACE products corrected with different repairing schemes to demarcate the catchment size that can be effectively observed. We have carried out such an analysis for 255 catchments of small to very large size in a closed-loop simulation environment, in order to answer the question of what is the minimum size of catchment that can be observed with the help of improved GRACE products.
We found that in general the error increases as the catchment size decreases, but the error level varies from one repair method to the other. GRACE time-series corrected with multiplicative approach shows the highest amount of error and the GRACE time-series from the data-driven approach has the lowest amount of error. We discussed the effective resolution of corrected GRACE products with respect to an acceptable error threshold. We found that the data-driven method and scaling method were able to approach the ideal GRACE resolution with an acceptable error RMSE level of 2 cm. The data-driven method has minimum divergence in terms of the spatial resolution with respect to the performance indicators (RMSE and NSE seas ).
Therefore, based on our investigation and its findings, we recommend the following: 1.
The spatial resolution of the band-limited GRACE spherical harmonics is not the half-wavelength at the equator, but the ideal spatial resolution. The spatial resolution of filtered GRACE data can also be described by the ideal spatial resolution.

2.
The users have to be wary that the spatial resolution of the corrected dataset is dependent on the adopted method. Furthermore, the spatial resolution is associated with the error tolerance required by the application and it has to be defined by the user.

3.
Given the fact that with enhanced processing techniques GRACE is able to see some catchments smaller than the spatial resolution and many catchments close to the band-limit resolution, it is worthwhile to provide spherical harmonics up to a maximum degree of 120 or higher.
The impending launch of the GRACE-Follow On mission, which is a near replica of the GRACE mission with only the Laser Ranging Instrument as the additional instrument, raises the question of the relevance of this study for GRACE-FO data. Recently, a simulation study by Flechtner et al. [49] indicated that the expected improvement from GRACE-FO is rather moderate, in which case we expect our quantitative results to hold even for GRACE-FO data. Having said that, the mathematical foundation and the understanding of the idea of spatial resolution developed in this study will always be relevant for data disseminated in terms of spherical harmonic coefficients.
P lm (cos θ) cos mλ B nkc lmcPnk (cos θ ) cos kλ + P lm (cos θ) cos mλ B nks lmcPnk (cos θ ) sin kλ + P lm (cos θ) sin mλ B nkc lmsPnk (cos θ ) cos kλ + P lm (cos θ) sin mλ B nks lmsPnk (cos θ ) sin kλ , where f (θ , λ ) is the noisy field; b(θ, λ, θ , λ ) is the filter kernel; C nk ,S nk and B nkc lmc , B nks lmc , B nkc lms , B nks lms are the corresponding spectra of the field and the filter;P lm (cos θ) are the normalized associated Legendre functions; andf (θ, λ) is the filtered field [8]. The filter function b(θ, λ, θ , λ ) is a two-point function that describes the relative weight between the calculation point (θ, λ) and the data point (θ , λ ). The calculation point is where we want the filtered f (θ, λ), which is computed by averaging all the data points by applying the weights b(·, ·) on the data points. The filter function can also be represented in polar co-ordinate system b(θ, λ, ψ, A), where the pole is the calculation point, and the data points are related to the calculation point by ψ the angular distance and A the azimuth.
The filter function described in Equation (A3) is the most general form of the filter, where the weights of the data points depend on the location of the calculation point, the angular distance as well as the azimuth. However, in the case of the GRACE data processing, the most commonly used Gausssian filter is a special type of filter that is classified as homogeneous isotropic filter. The weights of such special filters depend only on the spherical distance between the calculation and data points. They are independent of the location of the calculation point (homogeneous) and also of the azimuth (isotropic) (cf. Figure A1).
The spectrum of the homogeneous isotropic filters, due to their dependence only on the spherical distance, becomes dependent only on the spherical harmonic degree [14,50]: where P l (cos ψ) are the unnormalized Legendre polynomials of spherical harmonic degree l and B l are the coefficients of the homogeneous isotropic filter. Notice that there is only one index as opposed to four in Equation (A3). A field filtered with a homogeneous isotropic filter will have the following spectrum: where we see that the unfiltered spectrum C nk ,S nk is scaled by the filter spectrum B l for every spherical harmonic degree, as opposed to the general form of the filter (cf. (A2)). Equation (A5) is the same as Equation (7), and, therefore, all the band-limited fields have the properties of a field filtered with a homogeneous and isotropic filter. Since the fields filtered with a homogeneous isotropic filter have a homogeneous and isotropic spatial resolution, the spatial resolution of the band-limited fields is also homogeneous and isotropic.