The Gravity Recovery and Climate Experiment (GRACE) satellite mission measures the combined gravity signal of several overlapping processes. A common approach to separate the hydrological signal in previous ice-covered regions is to apply numerical models to simulate the glacial isostatic adjustment (GIA) signals related to the vanished ice load and then remove them from the observed GRACE data. However, the results of this method are strongly affected by the uncertainties of the ice and viscosity models of GIA. To avoid this, Wang et al. (Nat Geosci 6(1):38–42, 2013. https://doi.org/10.1038/NGEO1652; Geodesy Geodyn 6(4):267–273, 2015) followed the theory of Wahr et al. (Geophys Res Lett 22(8):977–980, 1995) and isolated water storage changes from GRACE in North America and Scandinavia with the help of Global Positioning System (GPS) data. Lambert et al. (Postglacial rebound and total water storage variations in the Nelson River drainage basin: a gravity GPS Study,
Geological Survey of Canada Open File, 7317, 2013a, Geophys Res Lett 40(23):6118–6122, https://doi.org/10.1002/2013GL057973, 2013b) did a similar study for the Nelson River basin in North America but applying GPS and absolute gravity measurements. However, the results of the two studies in the Nelson River basin differ largely, especially for the magnitude of the hydrology signal which differs about 35%. Through detailed comparison and analysis of the input data, data post-processing techniques, methods and results of these two works, we find that the different GRACE data post-processing techniques may lead to this difference. Also the GRACE input has a larger effect on the hydrology signal amplitude than the GPS input in the Nelson River basin due to the relatively small uplift signal in this region. Meanwhile, the influence of the value of \(\alpha \), which represents the ratio between GIA-induced uplift rate and GIA-induced gravity-rate-of-change (before the correction for surface uplift), is more obvious in areas with high vertical uplift, but is smaller in the Nelson River basin. From Gaussian filtering of simulated data,
we found that the magnitude of the peak gravity signal value can decrease significantly after Gaussian filtering with large average radius filter, but the effect in the Nelson River basin is rather small. More work is needed to understand the effect of amplitude restoration in the post-processing of GRACE g-dot signal. However, it is encouraging to find that both the methodologies of Wang et al. (2013, 2015) and Lambert et al. (2013a, b) can produce very similar results if their inputs are the same. This means that their methodologies can be applied to study the hydrology in other areas that are also affected by GIA provided that the effects of post-processing of their inputs are under control.

The Gravity Recovery and Climate Experiment (GRACE) satellites, launched in March 2002, aim to measure the Earth’s gravity change caused by the redistribution of mass within the Earth systems, including the oceans, solid Earth, atmosphere and hydrosphere, with a spatial resolution of approximately 300 km every 30 days. Because of the ranging technique used between the two twin satellites and their relatively low orbit, GRACE time-variable gravity data are not only able to reveal more spatial details than previous satellite gravity missions, but also allow gravity changes (e.g. due to global changes) to be monitored over time (Wahr et al. 1998).

The GRACE signal observed over North America and northern Europe has two main contributions: one is from the mass redistribution due to the glacial isostatic adjustment (GIA) process, and the other is from today’s hydrology or water storage change. Our interest is in identifying the latter (Lambert et al. 2013a, b; Wang et al. 2013, 2015). A common approach is to employ numerical models to simulate the GIA contribution and then remove them from the observed GRACE data. However, the results of this approach are strongly affected by the uncertainties of the ice and mantle viscosity models of GIA (Steffen et al. 2010; Wang et al. 2013). To avoid this, Wang et al. (2013, 2015) modified the method of Wahr et al. (1995) to investigate the water storage change in North America and northern Europe. Lambert et al. (2013a, b) also did a similar study for the Nelson River basin in North America. However, their results on the Nelson River basin are different, especially for the magnitude of the hydrology signal (see Fig. 4a, b below) which differs by about 35%. This disagreement is problematic because it could mean that there is an error in the hydrology signal separation procedure or that the separated hydrology signal errors have been underestimated. Thus, the purpose of this paper is to find out the causes of their difference. Hereafter, we shall refer to Wang et al. (2013, 2015) as W1315 and Lambert et al. (2013a, b) as L2013. In this paper, the focus is on the hydrology change in the Nelson River basin, but if successful, the methodology can be applied to hydrological studies in other areas that are also affected by GIA.

In the next section, we give a detailed comparison between different aspects of these two works following the sequence in Table 1. This is followed by an analysis regarding the possible causes that may lead to the differences in the estimated hydrology signal; then, a conclusion will be given at the end.

The key differences between W1315 and L2013 are summarised in Table 1. In this section, the comparison shall follow the order listed in Table 1.

Data

W1315 used GRACE and GPS data for their method, while L2013 used absolute gravity data in addition to GRACE and GPS data. They both used the GRACE Release 4 monthly solutions from 2002 to 2011, but the time series in L2013 is slightly longer by 10 months.

For both W1315 and L2013, the GRACE data were computed to spherical harmonic degree and order 60, equivalent to a spatial resolution of about 335 km.

Post-processing

W1315 applied a de-striping filter and an isotropic Gaussian filter with 340 km averaging radius to reduce the unwanted noise. The linear trend is obtained from the time series of the local data which is assumed to consist of periodic signals from the annual, 2.5 year and S2-tide variations in addition to the trend. Then, the trend coefficients were transformed into gravity perturbation on a regular \(1^{\circ }\times 1^{\circ }\) grid. Figure 1 shows the original gravity rates of GRACE and GPS-derived vertical velocities from W1315.

On the other hand, L2013 implemented a two-step method to preserve the signal amplitude (Huang et al. 2012): first, the GRACE data are de-striped with the additional criterion to minimise the loss of signal; second, a non-isotropic Gaussian filter is implemented only for the residual coefficients after identifying the prevailing coefficients. Lastly, they incorporated a gravity trend grid with a resolution of \(0.5^{\circ }\times 0.5^{\circ }\).

Table 1

Comparison of input data and post-processing methods of the two studies analysing water storage change in the Nelson River basin

GRACE g-dot data

Figure 2 compares the GRACE g-dot between W1315 and L2013 for the Nelson River basin and surrounding areas. An inspection of Fig. 2 shows that the patterns are similar: with low values on the southwest and increasing values towards the northeast. However, the maximum value at the northeast corner is about 1.4 \(\upmu {\hbox {Gal}}/{\hbox {yr}}\) for Wang et al. (2013) and 1.8 \(\upmu {\hbox {Gal}}/{\hbox {yr}}\) for Lambert et al.(2015). Overall, the g-dot value in Fig. 2b is approximately 0.4 \(\upmu {\hbox {Gal}}/{\hbox {yr}}\) larger than that of Fig. 2a. This is mainly due to the different GRACE post-processing techniques, where L2013 used a two-step method to try to minimise the loss of signal. Another difference may be due to the fact that L2013 have included the degree 1 coefficient in their GRACE data, while W1315 did not. However, as shown in “Appendix 2,” the long-term trend contribution of degree 1 is less than 0.008 \(\upmu {\hbox {Gal}}/{\hbox {yr}}\) in the prairies, so this effect is too small to be noticed in Fig. 2.

Uplift rate data from GPS

As for the GPS data, W1315 derived the vertical velocities from continuous GPS data from 1993 to 2006 in North America. The reference frame used is IGb00 for GPS data in North America (Sella et al. 2007). Although the time periods of the GPS data and GRACE data are not the same, the constant trend in the GPS data is able to resolve those difficulties. For the areas beyond the GPS networks, W1315 used the GRACE data there to acquire the corresponding vertical displacement rates to avoid edge effects. Besides, they implemented the same idea to eliminate the influence of ice melting in Greenland and Alaska. On the other hand, L2013 combined the 27 continuous GPS data from 1996 to 2010 and 50 Canadian Base Network (CBN) data over time period 1995–2011, using a 2-D adaptive Gaussian interpolation function to obtain a grid of \(2^{\circ }\times 2^{\circ }\) vertical velocity map. Moreover, they calculated the GPS-derived g-dot corrections for a larger area to avoid the possible edge effects. The GPS-derived maps for vertical velocities from W1315 and L2013 are shown in Fig. 3. It can be seen that both maps are quite similar in the northeast, but the negative contours and the contour shape in the west and south are different.

Absolute gravity data

In addition to GRACE and GPS data, L2013 also used absolute gravity surveys at 8 sites in mid-North America during 1987–2011. After removing the effects of local hydrology from the absolute gravity data, they are combined with collocated GPS data to infer the \(1/\alpha \) value in the empirical relationship between the GIA-induced uplift rate and g-dot which varies with geographic location and viscoelastic properties. Since the absolute gravimeters make measurements on the deformed surface of the Earth, L2013 corrected for the free air effects before using it to infer the hydrology signal.

Theory and methodology

Theory

Both approaches of W1315 and L2013 are based on the empirical relationship derived by Wahr et al. (1995):

Suppose we operate both a GPS receiver and a gravimeter near the centre of rebound and if the observed GPS vertical velocity signals are only caused by GIA, then from this empirical relationship, the GIA-induced g-dot can be obtained from the GPS observations (\(\frac{{\mathrm {d}} u}{{\mathrm {d}} t}\)). Therefore, the gravity signal due to water storage change (or hydrology) can be obtained by removing the GIA-induced g-dot signal from the observed GRACE data:

The assumption that the GPS vertical velocity signals are only caused by GIA has been verified in the centres of GIA in Canada and northern Europe (Wang et al. 2013; van der Wal et al. 2011).

Methodology

The methodology of W1315 is based on a spectral inversion method that uses the equation: \(\dot{\sigma }_{lm}^{\mathrm{TW}}= \frac{1}{\beta _{l}}\left( \dot{u}_{lm}^{\mathrm{GPS}}+{\alpha _{l}} \Delta \dot{g}_{lm}^{\mathrm{GRACE}}\right) \) where \(\dot{\sigma }_{lm}^{\mathrm{TW}}\) are the rate-of-change of the coefficients of surface mass density of terrestrial water, \(\dot{u}_{lm}^{\mathrm{GPS}}\) are the coefficients of GPS uplift rate, \(\Delta \dot{g}_{lm}^{\mathrm{GRACE}}\) are the coefficients of g-dot as observed by GRACE, \(\beta _{l}\) depends on the density structure and the yield properties of the Earth, and \(1/\alpha _{l}\) is a scaling constant determined theoretically to be \(+\,0.150\,\upmu {\hbox {Gal}}/{\hbox {mm}}\) after correcting for the uplift effect (Wang et al. 2015; Wahr et al. 1995). The equivalent water thickness (EWT) of the water storage change is obtained by dividing the surface mass density by the density of water. The result of such inversion is shown in Fig. 4a.

Fig. 4

Hydrology trend rates in equivalent water thickness (EWT) superposed on surface topography. a From W1315 with \(1/\alpha = -\,0.15\,\upmu {\hbox {Gal/mm}}\). b From L2013. c From W1315 with \(1/\alpha = -\,0.165\,\upmu {\hbox {Gal/mm}}\). d From W1315 with \(1/\alpha = -\,0.177\,\upmu {\hbox {Gal/mm}}\). Note the colour scale refers to the surface topography and the black contours are for the EWT trend rates

On the other hand, the methodology of L2013 is different. From gravity rate-of-change as determined by absolute gravity measurements at 8 absolute gravity stations that are collocated with continuous GPS stations in the area, the value of \(1/\alpha \) is found to be \(-\,0.165\pm 0.012\)\(\upmu {\hbox {Gal}}/{\hbox {mm}}\). Thus, using this value of \(1/\alpha \) and taking into account the free-air effect (\(-0.3086\,\upmu {\hbox {Gal}}/{\hbox {mm}}\)), L2013 converted GPS land uplift rate to g-dot due to GIA. The g-dot due to changes in terrestrial water storage can thus be obtained by subtracting the g-dot due to GIA from the “observed” g-dot signal from GRACE. The residual g-dot is then converted to EWT by a scaling constant \(24.3\,{\hbox {mm}}/\upmu {\hbox {Gal}}\) (Lambert et al. 2013b). The result of this method is shown in Fig. 4b.

Comparison of Fig. 4a, b shows that both have a peak near the centre of the map, in the upper Assiniboine River watershed east of Saskatoon. However, the peak amplitudes differ significantly, with a magnitude of \(34\,{\hbox {mm}}/{\hbox {yr}}\) in Fig. 4b but only \(24\,{\hbox {mm}}/{\hbox {yr}}\) in Fig. 4a. Furthermore, the high value centre in Fig. 4a extends further south than that in Fig. 4b.

In this section, we investigate the reason(s) why the results of L2013 differ from W1315 (see Fig. 4a, b). One possibility is due to the value of \(\alpha\) used. Another possibility is due to their difference in the g-dot and uplift rates used which may be due to the difference in their post-processing used. It is also possible that the uncertainties of the results have previously been underestimated. In the following subsections, we will study each of these possibilities individually.

Value of \(\alpha \)

As we noted earlier, the value of \(1/\alpha \) used by W1315 is \(-\,0.150\, \upmu {\hbox {Gal}}/{\hbox {mm}}\) which is slightly larger than the upper limit of range used in L2013 (i.e. \(-\,0.165\pm 0.012\,\upmu {\hbox {Gal}}/{\hbox {mm}}\)). So, can the difference in the inferred amplitude of the hydrology signal be due to the difference in the value of \(1/\alpha \) used?

To answer this, we used the methodology of W1315, but with \(1/\alpha \) equal the middle and lower limit of the range in L2013 (i.e. \(-\,0.165\) and \(-\,0.177\,\upmu {\hbox {Gal}}/{\hbox {mm}}\), respectively). The results are shown in Fig. 4c, d, which show that the value of \(1/\alpha \) has little effect on the signal amplitude in the Nelson River basin. This is due to the fact that the GPS uplift rate in the Nelson River basin is quite small (see Fig. 3). However, near the centre of GIA in Hudson Bay, where the uplift rate is much higher (the red area in Fig. 5d), the effect of the value of \(1/\alpha \) is more significant (see Fig. 5). For example, in the west and east side of Hudson Bay, the amplitude of the blue troughs is \(-\,10\) and \(-\,16\,{\hbox {mm}}/{\hbox {yr}}\), respectively in Fig. 5a, but \(-\,14\) and \(-\,20\,{\hbox {mm}}/{\hbox {yr}}\), respectively in Fig. 5c. Thus, the magnitude in the west and east side of Hudson Bay increases by 40 and 20%, respectively.

g-dot and uplift rates

If the value of \(1/\alpha \) is not the main cause of the difference in the inferred amplitude of the hydrology signal in the Nelson River basin, what else can be the cause? As we have seen earlier, due to the post-processing of the GRACE data, the g-dot values in W1315 are \(0.4\,\upmu {\hbox {Gal}}/{\hbox {yr}}\) smaller than that for L2013 (see Fig. 2a, b). This difference is approximately \(9.72\,{\hbox {mm}}/{\hbox {yr}}\) in equivalent water thickness (EWT) and can possibly explain the difference in the maximum value in Fig. 4.

To check this, we investigate the effects of the post-processed GRACE and GPS inputs on the resulting predicted EWT rates. This time, we follow the approach of L2013, but use different combinations of post-processed GRACE and GPS input data from W1315 or L2013. For visual comparison purpose, the EWT predictions along an east-west profile along \(52^{\circ }{\hbox {N}}\) are shown in Fig. 6. Here, the red solid curve is computed with the GRACE and GPS data both from L2013, and the blue solid curve uses GRACE and GPS data from W1315. On the other hand, the green dashed curve uses the data from Figs. 2a and 3b, while the orange dot-dash curve used the data from Figs. 2b and 3a.

Fig. 6

EWT profile computed with the method of Lambert et al. (2013a, b) using different GRACE & GPS input. See text for details

Inspection of Fig. 6 shows that: (1) the peak value of the red solid curve agrees with the peak value in Fig. 4b. This is expected because they all use the method of L2013, and their inputs are almost the same. (2) The values of the blue solid curve (including the peak) agree reasonably well with that along the profile at 52N in Fig. 4a. Note that the blue solid curve and the profile in Fig. 4a are computed with different approaches, but the inputs are almost the same. This means that both methodologies give almost the same EWT prediction when the input files are the same. (3) If the GRACE input of L2013 is used with the GPS input file of W1315, then the orange dot-dash curve shows that the peak is slightly displaced but the peak amplitude is slightly smaller than the red curve. On the other hand, when the GRACE input of W1315 is used with the GPS input file of L2013, then the green dashed curve is closer to the blue solid curve. This means that the GRACE input has a larger effect on the EWT amplitude than the GPS input. The reason is that at the location of the EWT rate (east of Saskatoon), the difference in the GPS data is relatively small, so its effect is smaller.

As for the difference in the pattern south of the peak (see Fig. 4a), the cause is probably due to the difference in the density and location of the GPS network used (see Table 1) which results in more negative contours in the south in Fig. 3b.

In short, the main cause of the difference in the peak amplitude of the inverted EWT rate between W1315 (Fig. 4a) and L2013 (Fig. 4b) is mostly due to the post-processing of the input GRACE data, rather than the methodology to determine the hydrology signal or the value of \(1/\alpha \).

Post-processing step

So, the next question is: which steps in the processing of GRACE data are responsible for the difference? From Table 1, we see that the main difference between the processing of GRACE data of W1315 and L2013 is that the latter tried to restore the signal lost in de-striping filtering, but the former did not. Furthermore, an isotropic Gaussian filter is applied to the GRACE data in W1315, while a non-isotropic Gaussian filter is applied in L2013. This could cause more amplitude reduction in L2013 than in W1315.

Thus, the next question is: how does the isotropic Gaussian filter in W1315 or the non-isotropic Gaussian filter of L2013 affect the amplitude of the GRACE g-dot data in the Nelson River basin? The answer of the first part may be found in “Appendix 1”, where it is shown that the magnitude of the peak g-dot values can decrease significantly after Gaussian filtering if the filter radius is large. However, these peak g-dot values lie outside the Nelson River basin. Figure 7 in “Appendix 1” also shows that for the synthetic GIA-induced g-dot data in the Nelson River basin, their amplitudes are not significantly affected, although the contours become smoother when higher frequency signal is reduced. Figure 8 in “Appendix 1” confirms that this is also true for the GRACE g-dot data that contains both GIA and hydrology contributions. Therefore, the isotropic Gaussian filter in W1315 should not be the main cause of the difference in amplitude of the GRACE g-dot data or the EWT rate around the Nelson River basin.

What about the non-isotropic Gaussian filter? A comparison of Fig. 4a, b shows that the pattern of the EWT results for W1315 resembles the topography map of the Nelson River basin more than that of L2013. Since basin topography should affect the pattern of the hydrology signal, the distortion of the hydrology signal may be caused by the non-isotropic Gaussian filter of L2013. This hypothesis should be verified in future studies.

Unfortunately, the synthetic and real GRACE g-dot data in “Appendix 1” cannot be used to study the effect of the de-striping filter on amplitude reduction nor the faithfulness of L2013’s statistical test in restoring the signal lost in de-striping because the whole procedure of L2013 requires the GRACE time series as input data but not the spatial distribution of g-dot data (private communication of J. Huang).

So far, we have traced the main source of difference between the amplitude of the hydrology signals inferred by W1315 and L2013 to be due to the difference in their post-processing of GRACE data. However, a thorough understanding of the effects of post-processing to the amplitude of the GRACE data is beyond the scope of this paper and should be addressed in future research.

Uncertainties

Other important questions that arise in view of that are now: how large are the uncertainties of the L2013 and W1315 results? Are their uncertainties large enough that the difference in the inverted hydrology signal can be explained? How large are the uncertainties of the inferred total water storage? “Appendix 2” shows that the uncertainty of the hydrology component in EWT is about \(2.3\,{\hbox {mm}}/{\hbox {yr}}\) for L2013 and \(2.5\,{\hbox {mm}}/{\hbox {yr}}\) for W1315. Since the peak EWT value is about \(34\,{\hbox {mm}}/{\hbox {yr}}\) in Fig. 4b and \(24\,{\hbox {mm}}/{\hbox {yr}}\) in Fig. 4a, that means the differences between the peak EWT value in Fig. 4a, b are larger than their uncertainties and thus enough to be resolved. Since we only consider the measurement uncertainties here, but do not consider the uncertainties from the separation approach and irregular GPS distribution, maybe we underestimate the uncertainties. The uncertainties details are illustrated in “Appendix 2”.

Overall, the publications W1315 and L2013 both recognise the decrease in water storage in the Canadian prairies from 2002 to 2004, as well as the recovery from a drought that started in 2004 and the floods that begin in 2010 due to heavy summer rainfall. Their comparisons with other hydrological data also verify that their monthly variations of the separated GRACE signal are related to changes in hydrology (Lambert et al. 2013b; Wang et al. 2013).

Through further analysis, we can see that the methodologies of W1315 and L2013 give similar EWT predictions in this part of Western Canada if the GRACE and GPS input files are the same. The difference between the amplitude of the EWT peaks between W1315 and L2013 is mainly due to the diverse post-processing techniques applied to the GRACE data: L2013 tried to restore the signal lost during de-stripe filtering while W1315 did not. In addition, Gaussian filtering on simulated GRACE data shows that the values in the Nelson River basin east of Saskatoon may not be significantly affected.

The difference in GPS land uplift rates also affects the amplitude of the predicted EWT, but its effect is relatively minor in the Nelson River basin compared to the GRACE g-dot data in this area with small uplift rates.

The good comparison between the methodologies of W1315 and L2013 means that their methodologies can be applied to the study of hydrology of other areas affected by GIA—although care must be exercised in getting the input data processed properly. However,
for the study of hydrology in Fennoscandia, the situation may be slightly different, since the centre of the separated hydrology signal in Fennoscandia is located at a region with high land uplift rate (Wang et al. 2013); thus, errors in the GPS measured land uplift rates and the empirical scaling factor \(\alpha \) will have larger effects on the separated hydrology signal results. Other than that the methodologies are expected to work fine.

Authors’ contributions

PW proposed and led the project; TL and
PW conducted the analysis and wrote the paper. HW and LJ processed and provided the GRACE, GPS and separated hydrology data. HW and HS contributed to the discussion and the preparation of the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We would like to thank two anonymous reviewers for their useful comments and suggestions. We are grateful to Dr. Jianliang Huang for his communication and for providing the data files for Figs. 2b, 3b and 4b. This study was conducted in support of the Groundwater Recharge in the Prairies (GRIP) project funded by Alberta Innovates. Patrick Wu is supported by GRF Grants 17315316 and 17305314 from the Hong Kong Research Grants Council. Hansheng Wang is supported by National Key R&D Program of China (2017YFA0603103), Key program of frontier science CAS (QYZDJ-SSW-DQC042 and QYZDB-SSW-DQC027) and National Natural Science Foundation of China (Grant Nos. 41431070, 41590854). We thank Liming Jiang for plotting the Graphical Abstract map.

Competing interests

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Not applicable.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Effect of Gaussian filter on signal magnitude

This appendix demonstrates the influence of isotropic Gaussian filtering on the magnitude of the GRACE signal. For the purpose of the first illustration in Fig. 7, the g-dot data used as input to the filter are synthetic and are computed from a spherical self-gravitating GIA model with ICE-6G_C (VM5a) (Argus et al. 2014; Peltier et al. 2015). The isotropic Gaussian filters we apply to the g-dot synthetic data in North America are 340 km averaging radius, 220 km averaging radius and 167 km averaging radius. The results are shown in Fig. 7b, c and d, respectively. The reason why we choose isotropic Gaussian filter rather than non-isotropic Gaussian filters is that the latter has the higher risk of removing real signals (Swenson and Wahr 2008).

From Fig. 7, we see that the magnitude of the peak decreases significantly after filtering: the peak value retains only 77.5, 88.4 and 92.5% of its original value with 340, 220 and 167 km averaging radius Gaussian filter, respectively. However, for the area around the Nelson River basin, the magnitude of the signal is not strongly affected, but the contours become smoother and high frequency noise is eliminated. Filtering may also eliminate some detailed information, like the small positive value centred at the north of Hudson Bay.

In the second illustration in Fig. 8, the g-dot input is the g-dot of real GRACE data and contains both GIA and hydrology contributions. The same Gaussian filters for Fig. 7 are also applied here. An inspection of Fig. 8 shows that the conclusions of Fig. 7 are also true for real GRACE g-dot data.

Therefore, the selection of a filter is important and that should depend on the situation of the study region and the aim of the investigation. Like in Fennoscandia, it has been revealed that Gaussian filter of radius 350–450 km yield superior results when GRACE data are used to determine the Earth viscosity structure (Steffen et al. 2010).

Appendix 2: Difference in uncertainty estimates for the two approaches

In this appendix, we describe how the uncertainties based on different data input are estimated. We follow the method described in L2013. The relation between the residual g-dot, where the hydrology component is estimated, and the g-dot from the GRACE data and GPS data are given by (2) which is rewritten here as:

Therefore, the uncertainty of the residual g-dot comes from two sources: (1) the uncertainty of the GRACE g-dot, and (2) the uncertainty of g-dot from GPS or GIA. According to the propagation law of the measurement errors, the residual g-dot uncertainty is given by:

Here, F represents the empirical conversion factor when g-dot from GPS vertical velocity (v) is derived at a fixed height. Thus, the error in GIA g-dot has two sources: (1) the error in F, (2) the error in GPS vertical velocity and can be expressed as:

Now, L2013 used \(-\,0.165\pm 0.012\,\upmu {\hbox {Gal}}/{\hbox {mm}}\) as \(1/\alpha \) value for the linear relationship between GPS vertical velocity and g-dot on the surface. So, after correcting for the effect of vertical surface motion (\(-0.3086 \upmu {\hbox {Gal}}/{\hbox {mm}}\)), the \(\frac{1}{{\alpha }'}\) becomes \(+\,0.1436\pm 0.012\)\(\upmu {\hbox {Gal}}/{\hbox {mm}}\) for L2013. However, the \(\frac{1}{{\alpha }'}\) is approximately \(0.15\, \mu{\rm{Gal}}\verb|/|{\rm{mm}}\) for W1315. As shown in Fig. 3, the GPS map for the two approaches is very similar, so we employ the same average GPS velocity value and GPS error for both two papers, \(v^{2}=23.16\,{\hbox {mm}}^2/{\hbox {yr}}^2\) and \(\sigma _{v}=0.1\,{\hbox {mm}}/{\hbox {yr}}\) (Lambert et al. 2013b). The errors and corresponding parameter values are listed in Table 2.

According to Table 2, the uncertainty of the residual g-dot differs around 7.2% in these two papers. If we convert the mean g-dot to equivalent water thickness through the factor \(24.3\,{\hbox {mm}}/\upmu {\hbox {Gal}}\) and integrate over 9.5 years, then we get the mean equivalent water thickness errors of 21.7 and 23.4 mm. In addition, W1315 did not include degree 1 in their GRACE data, so their separation method only applies to degree 2 to 60 (Wang et al. 2013; Wahr et al. 1995; Purcell et al. 2011). The effect of degree 1 on the EWT trend rates from 2002 to 2011 can be computed using the degree 1 coefficients of monthly gravity changes from the Tellus website (Swenson et al. 2008). The result is shown in Fig. 9, where the effect is less than \(0.2\,{\hbox {mm}}/{\hbox {yr}}\) or \(0.008\,\upmu {\hbox {Gal}}/{\hbox {yr}}\) in the prairies. This can be included in the uncertainty estimate for W1315’s separated EWT signal over 9.5 years, which is \(\sqrt{23.4^{2}+(0.2\times 9.5))^{2}}=23.5\,{\hbox {mm}}\).

Fig. 9

EWT trend rates of degree 1 GRACE gravity changes in North America from August 2002 to March 2011

Finally, numerically integration over the Nelson River basin reveals an error for the total water accumulation of 74.7 Gt for L2013 and 80.9 Gt for W1315 over 9.5 years from June 2002 to October 2011.

Note that we only consider the measurement uncertainties here, but the uncertainties from the separation approach and irregular GPS distribution which have been illustrated in Wang et al. (2013) have not been taken into account.