1 Introduction

The Gravity Recovery and Climate Experiment (GRACE) satellite mission is managed by the US National Aeronautics and Space Administration (NASA) and the German Aerospace Centre (DLR) (e.g., Tapley et al. 2004). Its goal is to map the Earth’s time-varying gravity field with a temporal resolution of 1 month or better and a spatial resolution of about 400 km. The mission was launched in March 2002 and is expected to provide data until at least 2010. Temporal gravity variations at these spatial and temporal scales are mainly caused by mass redistribution in the atmosphere and oceans, tides, postglacial rebound, and terrestrial water cycling. Through the data pre-processing, the contributions of tides, atmosphere, and oceans are largely removed using models of the underlying geophysical processes. Therefore, monthly gravity models mainly reflect changes in terrestrial water storage, as well as the change in snow/ice mass of the polar ice sheets and mountain glaciers, with respect to a long-term mean.

Monthly gravity field solutions are computed at the University of Texas at Austin Center for Space Research (CSR), the GeoForschungsZentrum Postsdam (GFZ) and the Jet Propulsion Laboratory (JPL) (see http://podaac.jpl.nasa.gov/grace or http://isdc.gfz-potsdam.de/grace), the Centre National d’Etudes Spatiales (CNES) (Lemoine et al. 2007), and the Delft Institute of Earth Observation and Space Systems (DEOS) at Delft University of Technology (Liu 2008), among others. These solutions are expressed in the form of spherical harmonic coefficients up to some maximum degree, representing monthly average values, although temporal sampling and averaging intervals are not completely uniform. The analysis centres follow different data pre-processing, processing, and post-processing strategies, which cause differences in the monthly sets of spherical harmonic coefficients.

Some analysis centres such as DEOS (e.g., Klees et al. 2007a), the Institute of Theoretical Geodesy (ITG) at the University of Bonn (e.g., Eicker 2008), the Goddard Space Flight Center (GSFC) (e.g., Rowlands et al. 2005; Luthcke et al. 2006), and others, compute solutions for specific areas of interest using spherical radial basis functions (SRBFs) or single layer densities (mascons) as an alternative to the spherical harmonic representation of surface mass change. These so-called regional solutions use overflight data exclusively over the region of interest. They are attractive for two reasons: (1) they require only a limited set of GRACE data to be processed and a relatively small number of basis functions, and (2) they are expected to exploit the resolution of the GRACE observations better than the spherical harmonic basis functions. Again, each analysis centre uses its own data retrieval procedure.

The main objectives of the paper are (1) to quantify the differences between the latest release of monthly spherical harmonic models of terrestrial water mass change, (2) to quantify the differences between monthly spherical harmonic models and monthly regional (SRBF and mascon) models, and (3) to investigate whether regional models provide more accurate estimates of water storage variations and better spatial resolutions than global models.

The paper is organized as follows: in Section 2, we briefly describe the GRACE gravity models, the main target areas, and the hydrological data used for GRACE model validation. The main results are presented in Section 3, which consist of a series of comparisons between regional and global GRACE models for different target areas, as well as a validation of these models. The latter consists of two parts: (1) we compare the GRACE models with the output of a regional hydrological model for the Zambezi River basin developed at Delft University of Technology; (2) we use the driest part of the Sahara desert as a validation area, with the assumption that the monthly mean mass change over this area is close to zero. In Section 4 we discuss the importance of a proper regularization or post-processing filtering of the regional GRACE solutions. The main results of the study are summarized and some conclusions are drawn in Section 5.

2 GRACE Gravity Models, Target Areas, and Hydrological Models

2.1 GRACE Gravity Models

We compare monthly global GRACE models provided by the CSR (GRACE Level-2 products, version CSR-RL04), GFZ (GRACE Level-2 products, version GFZ-RL04), JPL (GRACE Level-2 products, version JPL-RL04), CNES, and DEOS. The cut-off degree for most of the models is 120 with the exception of the CNES models, which is 50. The models are analyzed for the period February 2003 to February 2006. To allow for a fair comparison, we apply the same filters to all models in order to remove correlated noise in the spherical harmonic coefficients. We use the destriping technique of Swenson and Wahr (2006) in combination with a 400 km halfwidth Gaussian smoothing. This filter is referred to as DS400. The decision to apply the DS400 also to the CNES models has been made after a comparison with the unsmoothed CNES models. This comparison reveals that the unsmoothed CNES models still suffer from significant noise artefacts and perform worse compared with the DS400-filtered CNES models. The models are complemented by a global DEOS solution, which is filtered with an anisotropic non-symmetric (ANS) filter. The ANS filter has been designed to minimize the global mean of the mean-square error (Klees et al. 2008). The filter exploits full signal and noise variance–covariance information in an iterative least-squares approach. The corresponding model is referred to as DEOS-ANS, whereas the DEOS model that uses the DS400 filter is referred to as DEOS.

Two types of regional solutions are included in the study, which differ in terms of the parameterization used to model water mass change over the target area. The DEOS-SRBF solutions use Poisson wavelets (Holschneider et al. 2003; Holschneider and Iglewska-Nowak 2007) on an order-35 Reuter grid (Freeden et al. 1998), for a total of 1,542 basis functions. The order of the Reuter grid is an indication of the spatial resolution: order 35 corresponds to a cutoff spherical harmonic degree of 35. The optimal bandwidth of the Poisson wavelets was estimated in an iterative least-squares approach using Generalized Cross Validation (Golub et al. 1979). The approach is very similar to the one of Klees et al. (2007a) with the only difference being that no local refinement needs to be applied. The latter is justified by the relatively homogeneous coverage of the target area with satellite data. The solutions were regularized using the inverse of the ANS signal covariance matrix as a regularization matrix. This is equivalent to using the ANS post-processing filter, as shown in (Klees et al. 2008). The DEOS-MASCON solutions represent the time-varying gravity field as the potential of a single layer. The single layer density is modelled as a piecewise constant function over blocks of 4 × 4 deg. This is the same parameterization as that currently being used in the mascon solutions of GSFC (Rowland et al. 2005). The main difference between the DEOS-MASCON solutions and the GSFC solutions is that the latter use (1) a short arc analysis technique to relate GRACE range-rate data to the unknown mascons, and (2) an exponential signal covariance function in space and time to stabilize the solution and to enhance the spatio-temporal resolution. The DEOS-MASCON solutions are regularized in the same way as the DEOS-SRBF and the DEOS-ANS solutions.

The global and regional DEOS models use linear combinations of range measurements as observations in the functional model, which can be interpreted as average range accelerations. They are linearly related to the spherical harmonic coefficients, which express time variations in the gravity field. Assuming the noise in range observations is white, the noise in linear range combinations is colored, indicating that there are noise correlations in successive linear range combinations. These noise correlations were properly taken into account in the least-squares estimation procedure.

In all models to be studied the coefficients of degree 1 have been forced to zero and the coefficients of degree 2 were estimated together with the other spherical harmonic coefficients. The only exception are the CNES models, which use LAGEOS 1/2 satellite laser ranging data to fix the degree 2 coefficients.

2.2 Target Areas

The focus will be on three target areas: Zambezi (1.4 million km2), La Plata (3.0 million km2), and Sahara desert (3.0 million km2). Moreover, some results will also be presented for the Amazon (6.9 million km2), Mississippi (4.8 million km2), and Ob (3.0 million km2) drainage systems.

The Zambezi River basin is the fourth-largest basin of Africa, after the Congo/Zaire, Nile and Niger basins (see Fig. 1). The basin has a strong annual amplitude of about 0.20 m in terms of equivalent water height.

Fig. 1
figure 1

Target areas, which are used in the analysis: the Zambezi River basin and the part of the Sahara desert (left) and the La Plata River basin (right)

The La Plata River basin is the fifth largest river basin in the world (see Fig. 1) covering an area of about 3.0 million km2. The area is of interest because we expect that the GRACE estimates of the La Plata monthly mean water storage variations suffer from leakage of the Amazon River basin, the Tocantins catchment basin, and the San Francisco catchment basin, all of them being located to the north and north-east of the La Plata River basin. Leakage is caused by the application of a spatial filter (DS400 or ANS). Because the filter function is actually non-zero outside the target area (though its amplitude decreases quickly with increasing distance from the target area), mass variations outside the target leak into the target area, leading to a wrong estimate of the mean mass variation over the target area (see e.g., Klees et al. 2007b).

The part of the Sahara desert being used in this study comprises only the driest part of the desert. In particular, the Nile River basin in the eastern part and the Atlas mountain range in the western part of the Sahara desert have been excluded (see Fig. 1). The area extends over about 3.5 million km2 and will be referred to as Sahara desert for reasons of simplicity. When accepting that the monthly mean mass variations over this area are almost zero, then the GRACE estimates represent residual noise artefacts and/or leakage from surrounding areas, in particular from the Mediterranean Sea, the Intertropical Convergence Zone, and the Nile River basin. This would make the area well suited to validate the global and regional GRACE models.

2.3 LEW Regional Hydrological Model

For the Zambezi River basin, the Lumped Elementary Watershed (LEW) regional hydrological model is used as reference for the GRACE models. The LEW approach has been presented in a previous study by Winsemius et al. (2006) (see also Klees et al. 2007a, b). The LEW approach incorporates the redistribution of surface runoff in downstream located model units, called LEW’s, that represent e.g., a wetland, lake or man-made reservoir. The model has been forced by data from the Climate Research Unit (CRU) (New et al. 2002) and monthly runoff data from several data sources, including the Global Runoff Data Centre and the Zambian Department of Water Affairs. The LEW water storage estimates have been generated using rainfall estimates from the Famine Early Warning System (FEWS) (Herman et al. 1997). These estimates were lumped over the LEWs to provide a time series from February 2001 to February 2006. The first 2 years of simulation (February 2001 to February 2003) have been taken as a warm-up period to stabilize the state variables of the LEW model structures.

To assess the uncertainties in storage, related to uncertainties in the hydrological model, a Monte–Carlo framework was applied. In each simulation, uncertainty was assumed in the model’s parameters, for which storage and rainfall input, \(\hat{P,}\) are the most sensitive. The simulations were carried out in two steps. First, the parameters were perturbed 10,000 times. Each perturbed parameter set was used to run the model. Since the model is calibrated on discharge records, only the model results that gave a satisfactory resemblance with the calibrated model in terms of discharge were selected as being behavioral. With respect to model calibration, this yields a realistic simulation of parameter uncertainty, because prior knowledge about the discharge is now used to constrain it. In the second step, input uncertainty was considered. True rainfall was assumed to be:

$$ P = \hat P + \varepsilon_P,\quad \varepsilon_P = \varepsilon\,P, $$
(1)

where \(\varepsilon \sim {\mathbf{N}}(0,0.05).\) All the behavioral models were run with randomly perturbed rainfall, following the error model of Eq. 1.

The LEW model has a box-type conceptual model structure, in which rainfall is routed through a number of boxes where processes such as evaporation, runoff partitioning and stream discharge are represented by relations with the storage in the boxes. The parameters of the model and their estimation procedure is elaborately described by Winsemius et al. (2006), while the spatial discretization of model units is given in Klees et al. (2007b). The parameters that have the most influence on storage and its variability are typically the ones that determine the storage variability of the unsaturated zone. Thus, the following parameters were perturbed:

  1. 1.

    D (mm/month), which determines the amount of the monthly rainfall that is intercepted by vegetation and top-soil. D was allowed to vary around the calibrated value by 50% (uniformly distributed) with a fixed minimum of 30 mm/month;

  2. 2.

    Su,max (mm), determining the maximum soil moisture capacity. Su,max was allowed to vary with 25% with a fixed minimum of 0 (a negative soil moisture storage does not exist);

  3. 3.

    B (−), describing the spatial variability of the soil moisture capacity. B was allowed to vary with 50% with a fixed minimum of 0.

In this manner, 4.4% of all 10,000 Monte–Carlo simulations were retained and used as a proxy for the model parameter uncertainty. These models were then used to provide simulations with perturbed rainfall. The resulting time-variable storage estimations, lumped over the whole basin, are used in the following sections to provide confidence bounds for comparisons with the different GRACE models.

3 Comparison and Validation

Figures 2, 3, and 4 show geographic plots of estimated water mass change over, respectively, Africa, the Sahara desert, and South America for April 2004. April is typically a month of extreme water storage variations compared to the yearly mean with maxima and minima for river basins on the northern hemisphere and southern hemisphere, respectively. The plots indicate some significant differences among the global and regional GRACE solutions in terms of amplitude (see e.g., the Orinoco River basin and the Amazon River basin) and spatial pattern (see e.g., the Sahara desert, the African intertropical convergence zone, and the Congo/Zambezi River basins).

Fig. 2
figure 2

Water mass change (m) over Africa for April 2004 as inferred from global and regional GRACE models. For comparison, the output of the regional hydrological model LEW for Southern Africa is shown, as well

Fig. 3
figure 3

Water mass change (m) over the Sahara desert for April 2004 as inferred from global and regional GRACE models

Fig. 4
figure 4

Water mass change (m) over South America for April 2004 as inferred from global and regional GRACE models

3.1 Sahara

In order to get more insight into these differences, we first analyzed the GRACE solutions for the Sahara desert. Figure 5 depicts a time series of monthly mean water mass change over the Sahara desert spanning the period from February 2003 to February 2006. With the exception of the CNES models, all GRACE solutions show comparable RMS signals between 0.008 m (DEOS-SRBF, DEOS-MASCON) and 0.014 m (DEOS) (see Fig. 6). The time series of CNES models differs significantly in terms of RMS signal (0.023 m) and temporal pattern. Striking is the large maximum amplitude of −0.06 m attained in March 2003, which is not visible in all other GRACE solutions.

Fig. 5
figure 5

Time series of monthly mean mass change (m) in terms of equivalent water heights for the Sahara desert. Top global GRACE solutions; bottom regional GRACE solutions

Fig. 6
figure 6

Sahara desert: RMS of monthly mean water mass change (m) over the period February 2003 to February 2006. The regional models DEOS-SRBF and DEOS-MASCON have the smallest RMS values, whereas the CNES models suffer from strong leakage effects

A further analysis of the 3-year time series reveals that the CNES time series is the only one with a statistically significant yearly cycle. The amplitude is 0.02 ± 0.006 m. This yearly cycle is likely the effect of leakage from surrounding areas, e.g., the intertropical convergence zone, the Mediterranean Sea, or the Nile River system. The maximum of the yearly cycle is attained around December/January each year, which indicates a dominant influence from the intertropical convergence zone. The yearly amplitudes of the other global and regional models do not differ statistically from zero. A representative example is the DEOS time series with an amplitude of 0.007 ± 0.004 m.

3.2 Zambezi

GRACE solutions of monthly mean water mass change over the Zambezi river basin seem to fit very well (cf. Fig. 7). There are no phase differences between the solutions. However, a closer look reveals that there are significant differences in amplitude during periods of extreme water mass change. For particular months, the differences between the GRACE solutions can take up values of 0.07 m in terms of equivalent water heights. This is about 45% of the total signal of 0.15 m water mass change observed over the period from February 2003 to February 2006.

Fig. 7
figure 7

Time series of monthly mean mass change (m) in terms of equivalent water heights for the Zambezi River basin from February 2003 to February 2006 as inferred from different global GRACE solutions. Top global GRACE solutions; bottom regional GRACE solutions

The 3-year time series of monthly mean water mass change is dominated by the yearly cycle. The results of fitting a yearly sinusoidal through the time series of the Zambezi River basin are shown in Table 1. The amplitudes fit fairly well, ranging from 0.088 m (DEOS) to 0.118 m (DEOS-SRBF) with standard deviations between 0.006 and 0.010 m. There are no statistically significant phase shifts. The a posteriori standard deviations range from 0.029 m (JPL) to 0.044 m (DEOS-MASCON). They give an upper bound of the error in monthly mean water mass change estimates for the Zambezi River basin, since the residuals still contain unmodelled (in particular semi-annual and shorter period) signal, which are not accounted for by a sinusoidal model.

Table 1 Zambezi River basin: fit of a yearly sinusoidal through a time series spanning the period from February 2003 to February 2006

An indicator of GRACE model quality is provided by a comparison with the output of the LEW regional hydrological model. Figure 8 shows the differences between the LEW model and the GRACE models for monthly mean water mass change over the period February 2003 to February 2006. The standard deviations range from 0.025 m (CNES, CSR) to 0.036 m (DEOS-MASCON) (cf. Fig. 10). The RMS differences, which are also shown in Fig. 10, are significantly larger. This is due to the fact that the mass variations of the LEW model and the GRACE models refer to a different datum, which introduces a constant shift between the LEW and GRACE model outputs. The standard deviations are slightly smaller than the a posteriori standard deviations of a yearly sinusoidal fit (cf. Table 1), because the latter contain also semi-annual or shorter period variability.

Fig. 8
figure 8

Monthly mean water mass change (m) over the Zambezi River basin between February 2003 and February 2006: difference between the output of the LEW regional hydrological model and global (top) and regional (bottom) GRACE models. The violet curve shows the LEW model. All curves are centred to zero mean value

Note that the standard deviations and the RMS of the differences refer to GRACE models which have not been corrected yet for a bias. A more realistic comparison with LEW is expected if the GRACE estimates of monthly mean water storage variations are corrected for this bias. Therefore, the bias has been estimated for all models using the method of Klees et al. (2007b). In the bias estimation procedure, the DEOS-ANS models were used to represent the signal. Figure 9 shows the estimated bias. Note that the bias introduced by the DS400 filter is significantly larger than the bias of the ANS filter; extreme values are about 0.03 m for the DS400 filter and about 0.01 m for the ANS filter. After bias correction, the fit with the LEW model improves significantly for the JPL and DEOS solutions; for the other models, no significant changes compared with the solutions without bias correction are observed (see Fig. 10).

Fig. 9
figure 9

Time series of the bias in monthly mean mass change over the Zambezi River basin from February 2003 to February 2006. ANS filter bias (blue); DS400 filter bias (brown). Units are in (m) of equivalent water height

Fig. 10
figure 10

Statistics of differences between GRACE and LEW estimates of monthly mean water mass change over the Zambezi River basin from February 2003 to February 2006: RMS without bias correction (blue), standard deviations without bias correction (orange), and standard deviations after bias correction (yellow). Units are in (m) equivalent water height

The annual amplitude of the LEW model is 0.129 ± 0.007 m; about 11% of the signal variance cannot be explained by the yearly cycle. The RMS of the un-modelled part is 0.031 m, which is about 20% of the peak amplitude. GRACE underestimates the annual signal amplitude over the period February 2003 to February 2006; the differences with respect to the LEW model range from −0.041 m (DEOS) to −0.011 m (DEOS-SRBF).

In order to decide whether the GRACE models differ significantly from LEW, we conducted a Monte–Carlo simulation to quantify the uncertainty of the LEW model. Figure 11 shows the range between the 5% quantile and the 95% quantile. A look at the corresponding GRACE estimates reveals that the latter fall outside this band during periods of maximum mass surplus (spring), whereas they fall inside the band during periods of maximum mass deficit (fall). This misfit could partially be attributed to the effect of leakage; however, GRACE amplitudes are significantly larger than the LEW model in spring 2005, and correcting GRACE for the effect of leakage would increase the GRACE amplitudes even further (see Klees et al. 2007b). We have no explanation for this exceptional discrepancy between the GRACE and LEW models for spring 2005.

Fig. 11
figure 11

Monthly mean water mass change over the Zambezi River basin. The yellow band indicates the range between the 0.05 quantile and the 0.95 quantile of the LEW model from Monte–Carlo simulations. Units are in (mm) of equivalent water height

Another indicator of goodness of fit between GRACE models and LEW model is the spatial RMS difference. We subdivided the Zambezi River basin into 4 × 4 deg blocks, computed the standard deviation over all blocks, and averaged over all months between February 2003 and February 2006. The blocks are the same as those used in the DEOS-MASCON solution. Moreover, a mean spatial correlation coefficient over the 3-year times series was computed analogously as a relative measure of fit. In this way, a potential bias and scaling factor in LEW is removed. The results are shown in Fig. 12. The fit of the GRACE models with the LEW model is rather poor; the standard deviations range from 0.061 (DEOS) to 0.088 (DEOS-MASCON), i.e., they are significantly larger than the standard deviations of the mean over the river basin. This poor fit is also supported by rather low values for the mean spatial correlation coefficient, which are around 0.60 for all models except for the DEOS-MASCON models, which have an even lower correlation coefficient of only 0.42. This implies that none of the GRACE models are able to resolve spatial patterns at or below 4 × 4 deg blocks, which corresponds to a 400 km half wavelength.

Fig. 12
figure 12

Zambezi River basin: RMS difference (top) and spatial correlation coefficient (bottom) between GRACE models and LEW model in units of (m) equivalent water heights computed over 4 × 4 deg blocks and averaged over the period February 2003 to February 2006. Every block has been centred to zero mean value

3.3 Intertropical Convergence Zone

A striking difference among the GRACE solutions is visible in the April 2004 solution for the intertropical convergence zone (see Fig. 2). The DEOS-ANS global solution and the DEOS-SRBF regional solution show significantly more signal in that area than any other model. A qualitative indication of whether this signal is physical or not can be obtained when rainfall data are accumulated over the period from mid-March 2004 to mid-April 2004, and are compared with the difference between the GRACE solutions for April 2004 and March 2004. Figure 13 depicts such a comparison. Accumulated rainfall data from the Tropical Rainfall Measuring Mission (TRMM) (Kummerow et al. 1998) show a strong mass change signal, which is clearly visible in the global DEOS-ANS solution and in the regional DEOS-SRBF solution. This effect is less pronounced in the other solutions, probably due to smearing effects caused by the application of the DNS400 filter. A similar situation holds for the DEOS-MASCON solution, although another filter has been used.

Fig. 13
figure 13

South-African intertropical convergence zone: accumulated TRMM rainfall data from mid-March 2004 to mid-April 2004 (upper left) and the difference between April 2004 and March 2004 CSR (upper right), DEOS-ANS (bottom left), and DEOS-SRBF (bottom right) solutions. The TRMM data indicate the presence of a significant mass change signal (accumulated rainfall). This signal is not visible in the CSR solution, but clearly visible in the DEOS-ANS and DEOS-SRBF solutions

3.4 La Plata

The plots for South America (cf. Fig. 4) are dominated by the strong mass redistribution signal in the Orinoco and Amazon river basins, which is the reason why the mass change patterns provided by the various GRACE solutions seem to be very similar. However, a closer look reveals that the amplitudes of the DEOS solution for Orinoco and Amazon are significantly smaller than the amplitudes provided by the other models. Also visible are the more pronounced differences among the models for the La Plata River basin compared with other river basins like the Zambezi, Amazon, and Ob River basins. When changing the color bar to emphasize mass variations outside the Amazon and Orinoco River basins (not shown here), we observe that all geographic plots show bumpy patterns and some East–West smearing effects, which may represent residual noise and, as far as the global models CSR, GFZ, JPL, DEOS, and CNES are concerned, the effect of the DS400 filter.

The maximum water mass change in metres of equivalent water height for the La Plata River basin over the period from February 2003 to February 2006 is about 0.04 m. There are a couple of periods of extreme water mass change with significant differences among the global models, e.g., March 2003, April 2004, and September 2005 (cf. Fig. 14). In March 2003, the estimated amplitudes range from −0.01 m (DEOS) to 0.07 m (GFZ); in April 2004 from −0.03 m (DEOS) to +0.04 cm (GFZ); in September from −0.04 (CNES) and +0.05 (DEOS).

Fig. 14
figure 14

Time series of monthly mean mass change (m) in terms of equivalent water heights for the La Plata River basin. Top global GRACE solutions; bottom regional GRACE solutions

Moreover, Fig. 14 reveals periods of significant phase shifts, e.g., during Spring 2004 and Winter 2005/2006. Correspondingly, the phases of the yearly sinusoidals that were fit through the 3-year time series do not match well (see Table 2). They range from 0.4 ± 0.3 months (GFZ) to 3.8 ± 0.3 months (DEOS) and 4.4 ± 0.2 months (DEOS-ANS). The phases of the JPL, DEOS, and DEOS-ANS models differ statistically significant from the phases of the CNES, CSR, and GFZ models. The phase shifts cannot be attributed to the post-processing filtering, because (1) the same filter (DS400) has been applied to all global models (except for DEOS-ANS), and (2) the phase shifts between DEOS (DS400 filter) and DEOS-ANS (ANS filter) are statistically not significant. Therefore, we attribute the phase shifts to the different GRACE (pre-) processing strategies. We observe also that the overall fit of the time series to a yearly sinusoidal is weak, in particular when we compare the results with the ones obtained for the Zambezi (Table 1), Amazon (Table 3), and Ob River basins (Table 5). Only between 27% (DEOS-ANS) and 62% (CSR) of the overall signal variance can be explained by a yearly sinusoidal. For the Zambezi River basin, values around 80% were obtained for all global models (cf. Table 1).

Table 2 La Plata River basin: fit of a yearly sinusoidal through the time series spanning the period from February 2003 to February 2006
Table 3 Amazon River basin: fit of a yearly sinusoidal through the time series from February 2003 to February 2006

We explain these larger differences among GRACE models for the La Plata river basin with complex hydrological cycle, which makes the results sensitive to the differences in data processing and post-processing. Berbery and Barros (2002) have shown that there are different precipitation regimes within the La Plata river basin. Each of these different regimes has a well-defined annual cycle but with different phases. Correspondingly, the mean annual cycle over the basin is rather small. The limited spatial resolution of the filtered GRACE solutions does not allow to distinguish between the different precipitation regimes. The significant spatio-temporal variability inside the La Plata river basin is supported by the very low value of the signal variance that can be explained by the yearly sinusoidal for the DEOS-ANS solutions (27%), because the spatial resolution of DEOS-ANS is usually better than that of the other global models due to a better performance of the ANS filter compared with the DS400 filter (Klees et al. 2008). Another reason for the larger differences is leakage from outside the river basin, in particular from the southernmost extension of the monsoon system, which is close to the northern boundary of the river basin. The larger differences could also be attributed to residual noise artefacts, due to the relatively small amplitude of monthly mean water mass change over the La Plata River basin (Klees et al. 2008).

3.5 Amazon, Mississippi, and Ob

We also analyzed the GRACE solutions for the Amazon River basin, the Mississippi River basin, and the Ob River basin. The Figs. 15, 16, and 17 show the monthly mean mass variations over the period February 2003 to February 2006. The Tables 3, 4, and 5 contain the statistics of the fit of a yearly sinusoidal. All solutions provide very similar amplitudes and no phase shifts relative to each other. The maximum difference among the various GRACE solutions of the yearly amplitude is 0.051 m (30%) for Amazon, only 0.006 m (13%) for Mississippi, and 0.014 m (25%) for Ob. The formal standard deviations for the yearly amplitude do not exceed 0.009 m (Amazon), 0.006 m (Mississippi), and 0.005 m (Ob), and cannot explain the larger differences between the yearly amplitudes. Extreme values of the RMS of the posterior residuals are 0.02 m (DEOS-MASCON) and 0.04 m (CSR) for Amazon, 0.02 m (DEOS) and 0.03 m (CNES) for Mississippi, and 0.01 m (DEOS) and 0.03 m (CSR) for Ob. For Amazon, more than 90% of the total signal can be explained by the yearly cycle; these values are significantly smaller for Mississippi and Ob (between 75 and 85%). The CNES time series for the Mississippi River basin do not fit well with the annual signal. Although the estimated amplitude is comparable with the amplitude of other GRACE solutions, only 37% of the total signal variance is explained by the annual signal, whereas typical values for other GRACE models are above 70%.

Fig. 15
figure 15

Time series of monthly mean mass change in (m) of equivalent water height for the Amazon River basin. Top global GRACE solutions; bottom regional GRACE solutions

Fig. 16
figure 16

Time series of monthly mean mass change in (m) of equivalent water height for the Mississippi River basin. Top global GRACE solutions; bottom regional GRACE solutions

Fig. 17
figure 17

Time series of monthly mean mass change in (m) of equivalent water height for the Ob River basin. Top global GRACE solutions; bottom regional GRACE solutions

Table 4 Mississippi River basin: fit of a yearly sinusoidal through the time series from February 2003 to February 2006
Table 5 Ob River basin: fit of a yearly sinusoidal through the time series from February 2003 to February 2006

4 The Effect of Filtering

A direct comparison of the DEOS time series with the DEOS-ANS time series shows the effect of filtering. Figures 2 and 4 indicate that the use of an ANS filter (DEOS-ANS) provides a smoother solution and, at the same time, improved spatial resolution and higher amplitudes than the DS400 filter. This is clearly visible in the South African intertropical convergence zone, on Madagascar island, in the Orinoco River basin, and in the La Plata River basin. For instance, the DEOS-ANS solution is the only global solution which resolves the mass change signal on Madagascar island for April 2004, whereas this signal is hidden in all global solutions that use the DS400 filter due to some pronounced smearing effects. This improved spatial resolution of the DEOS-ANS solution is due to the optimal filter, which exploits information about noise and signal variances and covariances. Leakage is more pronounced in the DS400 solutions and the application of a DS400 filter yields over-smoothed solutions. The superior performance of the DEOS-ANS models is also supported by a smaller RMS signal for the Sahara region (cf. Fig. 6). The ANS filter is extensively discussed and analyzed in Klees et al. (2008).

Regional solutions are attractive, because they require a parameterization of only the target area and the surroundings, only use overflight data, and allow for a more flexible adaptation of the parameterization to the signal over the target area. Global solutions have to be filtered in order to remove the high-frequency correlated noise. Filtering is mostly done as a post-processing operation (e.g., DS400); the ANS-filtered solution used in the DEOS-ANS, DEOS-SRBF, and DEOS-MASCON time series can also be interpreted as a regularized least-squares solution of GRACE K-band ranging data. One of the objectives of this study was to determine whether regularization (or post-processing filtering) is also needed when computing regional solutions. This question is justified, because one could argue that regional solutions do not require (additional) regularization, as long as the parameterization is suitably chosen. Note that selecting the number and halfwidth of SRBFs can always be interpreted as a regularization by parameterization. When mascons are used, the choice of the size of the block, represented by a mascon, is nothing else but a regularization by parameterization. On the other hand, when choosing an insufficient number of SRBFs (or too large of blocks for each mascon), the resolution of the solution may be suboptimal, leading to a loss of information. This is equivalent to over-smoothing. Alternatively, one may select a larger number of SRBFs with a smaller halfwidth and adjust for the instability of the normal equations by adding a regularization matrix. If the latter is chosen equal to the inverse of the signal variance–covariance matrix (as has been done for the DEOS-SRBF and DEOS-MASCON solutions presented in this study), the solution will be the same as when applying the optimal ANS filter as a post-processing operation.

Figure 18 shows three different SRBF solutions for April 2004: the previously described DEOS-SRBF solution (which uses 1,542 SRBFs of halfwidth 700 km), this same solution but without regularization, and a SRBF solution with only 502 SRBFs of halfwidth 1,050 km. The optimal halfwidths of 700 and 1,050 km, respectively, have been determined using the approach of Klees et al. (2007a, b). Figure 18 shows that using only the case with 502 SRBFs has a significant smoothing effect: most of the signal in the intertropical convergence zone is spread-out over a larger area, and the mass change signal on Madagascar island is not visible anymore. Using 1,542 SRBFs without additional regularization enhances the signal in the intertropical convergence zone and on Madagascar island, similar to the DEOS-SRBF solution. At the same time, however, the highly correlated noise is not sufficiently suppressed as the stripy patterns indicate. The results clearly show that in order to obtain optimal gravity models, regional solutions, like the global ones, require a careful regularization or post-processing filtering.

Fig. 18
figure 18

Three regional SRBF solutions for April 2004. Left SRBF solution with 1,542 basis functions (bandwidth 700 km) and optimal ANS regularization; middle: as before but without regularization; right SRBF solution with 502 basis functions (bandwidth 1,050 km) and no regularization

5 Summary and Conclusions

We analyzed a range of global and regional GRACE solutions to quantify the differences between them and to answer the question whether regional GRACE solutions are to be preferred to global ones. The main result of the study is that current regional modelling techniques do not perform better than global ones for the river basins investigated in this study, as long as a proper spatial filter is applied. This could be different for high-latitude target areas, because more data are available there than for equatorial regions. This data coverage may result in a better signal-to-noise ratio for particular areas, which can be better exploited using a regional parameterization. This is a subject currently being investigated.

More important than the question of a regional or global parameterization is the application of a proper regularization or a suitable post-processing filter. For instance, Klees et al. (2008) have shown that the ANS filter performs significantly better than the DS400 filter. The latter performs well at the mid-latitudes, but has some weaknesses at the equatorial regions and the higher latitudes (cf. Swenson and Wahr 2006). The differences between the DEOS solutions (which were filtered with the DS400 filter) and the DEOS-ANS solutions are significant. Moreover, the regional DEOS solutions perform worse if other filters or regularization matrices are used. This includes the exponential signal covariance function, which is being used in the regional mascon solutions of the Goddard Space Flight Centre (GSFC) (e.g., Rowlands et al. 2005) as several tests with this regularization matrix (not presented here) have shown.

When we take the differences between GRACE models as a lower bound for the accuracy, we believe that 0.02 m equivalent water height is currently a reasonable estimate of accuracy of GRACE estimates of monthly mean water storage variations over river basins above 1 million km2. The amplitude of the annual signal can be estimated with an accuracy of several millimetres equivalent water height. We could not find a relation between the accuracy and the size of the river basin, even though river basins which varied significantly in size, i.e., between 1 million km2 (Zambezi) and 6 million km2 (Amazon), were examined. In general, the annual signal accounts for 70–80% of the total signal variance, with only the La Plata River basin showing exceptions in both directions. For the La Plata River basin, between 40 and 70% of the signal cannot be explained by the annual cycle, which may be seen as an indicator of strong out-of-phase mass variations outside the river basin (leakage) or significant variability at shorter spatio-temporal scales inside the river basin.

A comparison of the various GRACE solutions with the output of the LEW regional hydrological model for the Zambezi River basin does not allow us to conclude which model is to be preferred for this particular target area. The monthly mean mass change estimates derived from the filtered CSR solutions are closer to the unfiltered LEW model than the other ones; on the other hand, the GFZ models show the highest spatial correlation coefficient, though the correlation coefficients of the other GRACE models are not much smaller. Moreover, the RMS differences over blocks of 4 × 4 deg are significantly larger than the RMS differences of the mean over the Zambezi River basin; this is the effect of residual noise, which is left after filtering. Moreover, during periods of maximum mass surplus, the GRACE models are outside the 5–95% quantile band of the LEW model, which can be attributed to leakage effects.

The Sahara test also has its shortcomings because the signal is close to zero and no information about annual or shorter period signals is provided. Nevertheless, the analysis of the Sahara data has revealed significant differences between the GRACE solutions. In particular, the CNES solutions suffer from strong leakage effects, which lead to a statistically significant annual amplitude of 0.02 m equivalent water height, whereas the other GRACE solutions show random signals with a RMS value around or even below 0.01 m.