A statistical filtering approach for Gravity Recovery and Climate Experiment (GRACE) gravity data

[ 1 ] We describe and analyze a statistical filtering approach for Gravity Recovery and Climate Experiment (GRACE) data that uses a parameterized model for the temporal evolution of the GRACE coefficients. After least squares adjustment, a statistical test is performed to assess the significance of the estimated parameters. If the test is passed, the parameters are used by the filter in the reconstruction of the field; otherwise, they are rejected. The test is performed, and the filter is formed, separately for annual components of the model and the trend. This new approach is distinct from Gaussian smoothing since it uses the data themselves to test for specific components of the time-varying gravity field. The statistical filter appears inherently to remove most of the ‘‘stripes’’ present in the GRACE fields, although destriping the fields prior to filtering seems to help the trend recovery. We demonstrate that the statistical filter produces reasonable maps for the annual components and trend. We furthermore assess the statistical filter for the annual components using ground-based GPS data in South America by assuming that the annual component of the gravity signal is associated only with groundwater storage. The undestriped, statistically filtered field has a c 2 value relative to the GPS data consistent with the best result from smoothing. In the space domain, the statistical filters are qualitatively similar to Gaussian smoothing. Unlike Gaussian smoothing, however, the statistical filter has significant sidelobes, including large negative sidelobes on the north-south axis, potentially revealing information on the errors, and the correlations among the errors, for the GRACE coefficients.


Introduction
[2] NASA's Gravity Recovery and Climate Experiment (GRACE) mission [e.g., Tapley et al., 2004] is providing a unique opportunity for studying the Earth's global gravity variations. At seasonal timescales, the most significant contributions to gravity changes are associated with climate-driven mass transport of water, in its various forms, on the surface of the Earth. At longer timescales, mass transport below the surface of the Earth due to glacial isostatic adjustment (GIA) is also important. A number of studies have used GRACE data to examine the long-term and seasonal variability in the various components of the Earth's water balance, including continental water storage, cryosphere systems, and oceans [e.g., Chambers et al., 2004;Wahr et al., 2004;Davis et al., 2004;Tamisiea et al., 2005;Velicogna and Wahr, 2006;Swenson and Milly, 2006;Yeh et al., 2006], as well as GIA [e.g., Tamisiea et al., 2007].
[3] These analyses are accomplished using spherical harmonic (Stokes) coefficient data sets, each representing a reduction of approximately one month of GRACE data by the GRACE analysis centers. For CSR Release 4 (RL04) GRACE fields, these Level-2 data sets [Bettadpur, 2007b] are complete to degree and order 60, although the higher degrees and orders within this range are expected to be noisier [e.g., Wahr et al., 1998Wahr et al., , 2006] and therefore require some kind of filtering. A number of approaches have been suggested. The most straightforward approach, other than simply termination at some degree, is Gaussian smoothing. Jekeli [1981] developed the spherical harmonic coefficients for this filter, and Wahr et al. [1998] described their application to GRACE data. This has become the standard approach because it is conceptually straightforward, simple to implement, and at the suitable averaging radius (typically greater than 500 km) produces normally distributed errors .
[4] Other methods have been proposed. Swenson and Wahr [2002] and Seo and Wilson [2005] proposed smoothing functions that are adapted to the specific region being studied and optimized for prior estimates of data uncertainties. Han et al. [2005] described another anisotropic averaging function, based on the Gaussian averaging function, that uses an averaging radius that depends on the spherical harmonic order. The specific order dependence they advocated yielded an averaging function with a Gaussian profile that is elongated in the east-west direction (1000 km smoothing radius) compared to the north-south direction (500 km).
[5] The east-west elongation of the Han et al. [2005] averaging function deals with the problem of ''striping'' in the GRACE data, errors that manifest themselves as generally north-south linear features. Swenson and Wahr [2006] discovered that these stripes were associated with correlations among the errors of the GRACE coefficients for a given order and degree parity (even/odd), and they presented an ad hoc method for reducing these correlations. Their destriping procedure decreases significantly the stripes in the GRACE data and the averaging radius required, but it does not replace Gaussian smoothing. Constraining the adjustments for the Stokes coefficients [Lemoine et al., 2007] also seems to reduce the error associated with stripes. (Throughout this paper, we use the term ''destriping'' to refer to the a procedure based on that described by Swenson and Wahr [2006], and ''destriped'' to refer to fields that have been processed using this procedure. We refer to fields processed without this procedure as ''undestriped''.) [6] Davis et al. [2004] made use of a technique that used a statistical test to determine whether the time series for each spherical harmonic degree and order possessed a significant annual component. Those that did not were not used in the reconstruction of the field. Gaussian smoothing was not performed. In this paper, we describe the approach in detail, and extend it to estimate secular variations. We use GPS data to assess the accuracy of this technique relative to Gaussian smoothing, investigate the impact of destriping, and calculate the spatial representation of our filter. We begin by describing the theoretical basis for the filter.

Model for Temporal Variations of Gravity
[7] Changes in the Stokes coefficients from a GRACE Level-2 data set can be used to calculate changes in the geoid DN(f, l, t) at a point on the surface of the Earth with latitude f and longitude l and at epoch t [e.g., Wahr et al., 1998;Bettadpur, 2007b]: where a is the radius of the Earth, P 'm (sin f) are the fully normalized associated Legendre functions [Bettadpur, 2007b], DC 'm (t) and DS 'm (t) are changes in the dimensionless Stokes coefficients at degree ', order m, and epoch t. Similar expressions have also been developed that give the surface mass density or load thickness (when the load is assumed to be water) [Wahr et al., 1998] and the radial deformation associated with the surface mass load [Davis et al., 2004].
[8] If the geoid varies in time with both seasonal and secular influences, then we might model it using the expression where DN°is a reference value, D _ N is a rate term, DN a,i and DN a,o are annual amplitudes, and Dt = t À t o with t o an (integer year) reference epoch. (Superscript i stands for ''inphase'' o for ''out of phase''.) The frequency f is one cycle per year. The dependence on position (f, l) in equation (2) has been omitted in DN, DN o , D _ N , DN a,i , and DN a,o for clarity. Since ''seasonal'' climate variability is not purely annual, equation (2) in reality will represent an approximation, and unmodeled variability during the time span of the data set would be interpreted as noise.
[9] Because of the linearity of equation (1), equation (2) implies that each Stokes coefficient can be written as and similarly for DS 'm (t). We can therefore estimate the geoid rate and annual amplitudes by fitting equation (3) to the time series of Stokes coefficients for each degree and order, and then using equation (1) to reconstruct the spatial field. To linearly filter the field, the coefficient is first multiplied by the Stokes coefficient for the spherical harmonic expansion of the filter (C 'm F , S 'm F ) before reconstructing via equation (1). The filtered geoid rate field, for example, is [10] For an isotropic filter the coefficients are independent of order and C 'm F = S 'm F . The salient feature of the isotropic Gaussian smoothing filter [e.g., Jekeli, 1981;Wahr et al., 1998], for example, is a monotonic decrease in the filter coefficients with increasing degree.

Statistical Filtering Approach
[11] The Stokes coefficients used to derive the parameters of the time-dependent model are noisy, and this noise propagates into the fields reconstructed from the estimated parameters. The motivation for smoothing of the reconstructed fields is that the errors in the GRACE coefficients and derived quantities (annual amplitudes and rates) at high degrees are larger than at lower degrees, or at least larger compared to the real signal. However, this assumption may not necessarily be correct. Figure 1 shows time series for the GRACE coefficients C 3,0 and C 4,0 . Both are what would be termed ''low degree and order'' terms. The (3,0) term shows a distinct annual pattern, whereas the (4,0) term does not. In any reconstruction of the global annual gravity variation, we would of course wish to include the annual amplitudes for the (3,0) term. On the other hand, for the (4,0) term we might ask whether it is better to include the estimated annual amplitudes, or to assume that the value is zero.
[12] The problem can be approached as a statistical hypothesis test between two models: one in which all the parameters are freely adjusted (the ''full'' model, subscript f below), and one in which some of the parameters are constrained to zero (the ''restricted'' model, subscript r) [e.g., Peixoto, 1993;Bennett et al., 1996]. The appropriate statistic is an F statistic which may be written as where the c f 2 and c r 2 are the postfit unit weight residual sums of squares and the n f and n r are the degrees of freedom for the solutions. The unit weight residuals are calculated by normalizing the residuals by the standard deviation of the observation error. Under the null hypothesis represented by the restricted model, the statistic and F is F distributed with m = n r À n f (numerator) and n = n f (denominator) degrees of freedom [Peixoto, 1993]. (The F distribution describes random variables that are the ratio of two random variables each being c 2 distributed.) A large difference between the c 2 values for the restricted and full models would indicate that the parameters in question ''absorb'' residual power. This situation would lead to a large F value compared to the expectation of unity under the null hypothesis.
[13] In applying this hypothesis test to one or more parameters, we form the restricted model by fixing (at a value of zero) these parameters from the full model in equation (3). If we are interested in the rates, then the restricted model has the rate term equal to zero, and n r À n f = 1. If we are interested in the annual amplitudes, then those two terms are set equal to zero in the restricted model and n r À n f = 2. For each degree and order, we can perform weighted least squares adjustments using each coefficient (C 'm or S 'm ) for the full model in equation (3) and for the appropriate restricted model. We can then form the F statistic, and accept or reject the relevant parameters if F exceeds some value based on the likelihood for this occurring under the null hypothesis.

Results
[14] We applied our statistical filtering process to 57 GRACE RL04 Level-2 monthly data sets between epochs 2002.3 and 2007.3, produced by the Center for Space Research (CSR) at the University of Texas. The GRACE data are first corrected for mass variations in the atmosphere and oceans using models as described in Bettadpur [2007aBettadpur [ , 2007b. We used the calibrated standard deviations in the least squares fit. (Calibrated errors are described in the CSR RL04 Release Notes available at ftp://podaac.jpl.nasa.gov/ pub/grace/doc/ReleaseNotes_csr_RL04.txt.) The maximum degree used in the solutions is 60. We performed an adjustment of the full and restricted models for the time series of each Stokes coefficient, formed the and F statistic, and rejected parameters if the statistic does not exceed the 99.99% confidence level under the null hypothesis. Thus our filter coefficients (equation (4)) take on values of one (F exceeds 99.99% confidence limit) or zero (F fails to exceed this limit).
[15] We begin by investigating whether the F values are simply F distributed, which would indicate that they reflect only random error, or whether they deviate from an F distribution, indicating that a rejection of the null hypothesis is warranted in some cases. Figure 2 shows the observed distributions of F (lumping the cosine C and sine S terms together) and the theoretical distributions, based on the null hypothesis. For low values of F, the observed distributions decrease rapidly as do the theoretical distributions. The observed distributions have much longer tails, however, indicating that the parameter would be accepted as being nonzero more frequently than predicted under the null hypothesis. The 99.99% confidence limits for F for this data set are 10.98 for the annual term F test and 17.74 for the rate F test, and 0.01% of the trials (<1 for maximum degree and order 60) are expected to exceed these values under the null hypothesis. In fact, 877 values (in-and outof-phase combined) exceed the 99.99% limit for the annual terms and 330 for the rate. (See the auxiliary mat erial.) 1 Several F values for both the rates and the annual terms exceed 100. For comparison, the F values for the annual terms for the components shown in Figure 1 are 74.8 (3,0) and 7.1 (4,0). For the rate terms, the F values are 18.2 (3,0) and 28.2 (4,0). Thus, both D _ C 3,0 and D _ C 4,0 would be included in the filtered rate reconstruction, whereas only DC 3,0 a,i and DC 3,0 a,o would be included in the filtered annual amplitude reconstruction. In sections 4.1-4.4, we examine these reconstructions.

Statistical Filter for the Annual Components
[16] We have calculated the statistically filtered annual components for surface water thickness [Wahr et al., 1998] after applying the destriping algorithm of Swenson and Wahr [2006] and without destriping, and compared it to the fields determined using Gaussian smoothing of 200 km smoothing radius with and without destriping ( Figure 3). (The Gaussian smoothing used in this study is an implementation of the Wahr et al. [1998] algorithm.) The longitudinal striping effect is evident in the Gaussian smoothed, undestriped fields (Figure 3a), although the largest-amplitude features can be identified. Destriping the field prior to Gaussian smoothing (Figure 3b) successfully removes most of the striping effects, except for a few (reduced amplitude) stripes near the equator, consistent with the findings of Swenson and Wahr [2006]. Although the stripes are still visible even with a smoothing radius of 400 km (not shown in Figure 3), the RL04 CSR GRACE fields exhibit a marked improvement in this regard relative to the RL01 fields.
[17] The statistically filtered field shows no evidence of longitudinal stripes, even without applying destriping ( Figure 3c). Destriping the fields prior to statistical filtering ( Figure 3d) has no obvious effect. The maximum difference between the fields in Figures 3c and 3d is 100 mm, which is $12% of the peak-to-peak maximum amplitude. (The color scale in Figure 3 saturates at ±300 mm to stress low-level features. A color scale that would avoid saturation for all plots would range from -378 mm to 444 mm.) [18] As can be seen from Figures 3b and 3c, the statistical filtering technique and Gaussian smoothing both smooth the fields. The source of the similarities and differences can be inferred from Figure 4. In Figure 4, we show the spherical harmonic coefficients for the statistical filter for the annual components and for the isotropic Gaussian filter of smoothing radius 400 km. Since the statistical filter has a dependence on order, for each degree ' we calculate and plot in Figure 4 the average value H ' of the filter coefficients C 'm F and S 'm F : where the superscript i is for in-phase, and o for out of phase. For the statistical filter, the C 'm F and S 'm F take on the values of either zero or one. If all the coefficients for a particular degree are accepted by the filter, then the value of H ' would be one; if half are accepted, then H ' has a value of 0.5, and so on.
[19] Aside from isotropy versus anisotropy, the statistical filter differs from the Gaussian smoothing filter in several ways. The Gaussian smoothing filter is greatest for the lowest degrees. Coefficients for degrees 0 and 1 are not included in the GRACE solutions, so the statistical filter has no defined values at these degrees. However, the statistical filter for the annual components has a value of 0.2 at degree 2 for this data set, indicating that all degree 2 components except two (S 2,2 a,i and S 2,2 a,o ) are rejected since they failed the hypothesis test for annual variations. The (2,0) component was rejected by the filter along with most of the other degree 2 terms. (Recall that these coefficients are estimated relative to atmosphere and ocean models. These models themselves have large annual amplitudes for several low-degree terms.) [20] The statistical filter does maximize between degrees 3 and 6, although its value varies between 0.59 and 0.86. These low values relative to the Gaussian smoothing filter (for nearly all degrees) imply that a number of components are being rejected by the statistical filter. The statistical filter value drops more rapidly with increasing degree than the Gaussian smoothing filter with this smoothing radius (400 km).
[21] In the spatial domain, when destriping is not performed (Figures 5a and 5b) the statistical filter for the annual components is similar to the anisotropic Gaussian smoothing filter of Han et al. [2005] (Figure 5c). Whereas the profiles are similar to the Gaussian profile, the east-west axes are elongated relative to the north-south axes. The reconstruction of the statistical filter is position-dependent, so in Figure 5 we evaluate it at two positions on the equator that illustrate the C 'm F terms (0°longitude) and S 'm F terms (W 90°longitude). Also shown for comparison is the anisotropic Gaussian smoothing filter of Han et al. [2005] evaluated on the equator with smoothing radii r = 400 km, r 1 = 800 km, and with m 1 = 15. When destriping is performed, combined representation of the destriping and the statistical filter (Figures 5d and 5e) has a smaller main lobe footprint than the anisotropic Gaussian filter, approaching the footprint of the isotropic Gaussian (Figure 5f). However, the statistical filter contains negative sidelobes, especially on the north-south axes, and considerable power in the NE-SW and NW-SE directions. This feature is shared with the Gaussian-smoothed destriping filter [Swenson and Wahr, 2006, Figure 4].
[22] Despite the radical differences between the statistical filters from the striped and destriped fields, the maps of the annual amplitudes that are produced are similar (Figures 3c  and 3d). Clearly, the destriping process is allowing the statistical filter to perform spatial averaging over a much smaller area, as it should. The destriping process also introduces other effects. From the negative sidelobes in the statistical filter, we infer that the destriping process may add negative latitudinal correlations. From the off-axis lobes, we infer that the destriping process does not remove some positive correlations that are not strictly longitudinal.

Comparison of Statistically Filtered Annual Gravity Field and GPS
[23] The expression for the radial deformation associated with a gravity change that is assumed to be due to a surface mass [Davis et al., 2004] is where k ' E and h ' E are the elastic Love numbers. The deformation inferred from the GRACE data can therefore be compared to determinations of crustal deformation made using (for example) the Global Positioning System (GPS) [Davis et al., 2004;King et al., 2006;van Dam et al., 2007]. Davis et al. [2004] used the large gravity variations in South America due to annual water storage in the Amazon Basin to compare the annual component of the radial deformation inferred from GRACE to that measured at GPS sites in the . Spherical harmonic coefficients, averaged over all orders for each degree (see text), for the statistical filter described in the paper for the annual components (solid line) and the isotropic Gaussian smoothing filter with smoothing radius 400 km (dashed line). The coefficients are normalized so that a filter that passes everything would have a uniform value of 1.
region. The comparison requires estimation of the geocenter motion present in the GPS data but not in the GRACE data.
[24] We can use this comparison to assess the annual component of a GRACE-derived deformation field, and therefore any method of filtering that is used in the reconstruction of this deformation. In particular, we have assessed the effects of destriping, isotropic and anisotropic Gaussian smoothing, and the statistical filtering technique. We have also performed the assessment on a range of smoothing radii. Figure 6 shows several of the reconstructed radial deformation fields for the study region. Increasing the smoothing radius has a significant effect on the Gaussian smoothed fields, as expected. Destriping also has a significant effect, and the stripes are almost absent from the fields smoothed with a radius of 250 km. The anisotropic smoothing appears at first glance to have an insignificant effect, but in fact this radius refers to the larger smoothing radius. We adopted a value for the smaller radius that was 50% of the larger. The anisotropic Gaussian smoothing filter thus has an averaging area $50% that of the isotropic Gaussian smoothing filter with the same labeled smoothing radius.
[25] As noted earlier (see Figure 3), the destriping technique does not produce much difference in the statistically filtered annual components (Figure 6, right), despite con- siderable differences in the resulting filters themselves ( Figure 5). There are, however, observable differences between the statistically filtered fields and the smoothed fields, especially for the region south of latitude 20°S.
[26] To assess these fields, we used estimates of in-phase and out-of-phase annual radial deformation from GPS data acquired at the sites shown in Figure 6. For each field, we calculated the reduced c 2 residual between the GPS estimates and the GRACE-determined values after estimation of the geocenter correction [Davis et al., 2004]. (King et al. [2006] employed a similar method using a global network of GPS sites.) The GPS uncertainties were the only contribution to the variance used in the calculation of the c 2 residual, as they are assumed to be much larger than the GRACE uncertainties. The smoothing radii (the r 1 values for the anisotropic Gaussian smoothing filter) were varied between 0 and 1000 km. The values r = r 1 /2 and m 1 = 15 were used for the anisotropic Gaussian smoothing filter. The results are shown in Figure 7.
[27] The undestriped anisotropic Gaussian smoothing filter (solid red line in Figure 7) produces a field that agree with the GPS measurements better than the isotropic Gaussian smoothing filter (solid blue line) does at radii greater than $250 km. At radii below 250 km, there is little difference between isotropic and anisotropic Gaussian smoothing. The destriped smoothed fields (red and blue dashed lines) agree with GPS better than the undestriped smoothed fields (red and blue solid lines).
[28] The destriped anisotropically smoothed fields (red lines) show little variation in the c 2 assessment for a large range of radii, $300 km to $700 km. However, as Swenson and Wahr [2006] point out, the destriping process affects the underlying geophysical signals as well as reducing the correlated errors that lead to the stripes. If we had the complete field of GPS-determined radial deformation, we could apply destriping to that field for a more appropriate comparison. This limitation of the destriping method therefore also imposes a limitation on our assessment. Figure 6. Radial deformation amplitudes (left, in-phase; right, out-of-phase) for the South American region studied by Davis et al. [2004]. The location of the GPS sites are shown as black triangles. For each smoothing radius, four Gaussian-smoothed fields are calculated. Top to bottom they are isotropic, no destriping; anisotropic, no destriping; isotropic, destriping applied; anisotropic, destriping applied. For anisotropic smoothing, the radius refers to r 1 . The other parameters are r = r 1 /2 and m 1 = 15. At right are the statistically filtered fields with no destriping (top) and destriping applied (bottom).
[29] Because of these limitations, this comparison can provide an indication of the level of agreement with GPS of the statistical filter in relation to the smoothing process, but not of the accuracy. The undestriped statistically filtered fields (solid green line) perform equivalently when compared to the best agreeing anisotropically Gaussian smoothed fields, except when r 1 < 250 km, where the smoothed fields have poor agreement with GPS. (This result differs slightly when we use CSR RL01 fields. The reduced c 2 for the best agreeing smoothed field is $15% greater than undestriped statistically filtered field.) As might be inferred from the fields for the annual components themselves, destriping has only a small effect on the statistical filter. In contrast to the Gaussian smoothed fields, however, destriping worsens the agreement for the statistical filter. Given the ''pick and choose'' process that leads to the statistical filter, the reconstructed fields ( Figure 5) look remarkably smooth. The assessment using GPS data indicates that the (data-derived) statistically filtered fields perform as well or better (in terms of agreement with GPS) than fields that are smoothed with a filter based on prior assumptions regarding the GRACE errors. Indeed, the statistical filter for the annual components appears to inherently reduce errors associated with stripes.

Effect of Statistical Filter on Model Fields
[30] To assess the possible distortion that the statistical filter might have on observed fields, we generated a realistic simulation of the annual amplitude fields using a scaled version (see below) of the GLDAS/Noah soil moisture fields [Rodell et al., 2004]. We decomposed each monthly GLDAS field into spherical harmonic coefficients, and estimated annual amplitudes and trends using the same computer code used to analyze the GRACE data. The resulting annual amplitude fields (Figure 8a) were used as the model for our simulation of the effects of filtering. The similarity of the model and observed annual amplitude fields (e.g., Figure 3c) clearly indicates that hydrology is the dominant process at the annual period. Some ''ringing'' is observable in the model fields, associated with our truncation of the expansion at degree 60. We also observe a low-amplitude signal in the North Atlantic and Arctic oceans associated with our imposition of a self-consistent sea level in the land-ocean exchange of water, and with the omission of the degree 1 component.
[31] To simulate the response of our filter to the model annual fields, we assumed that the model field of Figure 8a represented ground truth and that the only variability in the monthly fields were the annual variation based on this model field and random noise from GRACE. Under these conditions, it can be shown that an approximate expression for the expected value of F 'm C (the F statistic of equation (5) for the GRACE cosine coefficient at a specific degree and order) for a large number n of GRACE fields is and similarly for F 'm S . The annual in-phase and out-of-phase amplitudes were defined in equation (3). We assume in equation (8) that the observation standard deviation s 'm for a GRACE coefficient is constant with time. Equation (8) intuitively takes the form of the squared ratio between the annual amplitude (the signal) and the measurement uncertainty (the noise), i.e., the signal-to-noise ratio (SNR). As discussed earlier, the statistic of equation (8) is Figure 7. Assessment of the reconstructed annual deformation in South America using the GPS sites of Davis et al. [2004] using the reduced c 2 . Blue indicates isotropic Gaussian smoothing filter. Red indicates anisotropic Gaussian smoothing filter [Han et al., 2005]. Green indicates statistical filter. Solid indicates undestriped. Dashed indicates destriped. For the anisotropic Gaussian smoothing filter, the abscissa refers to r 1 (see text).
F distributed with (2, n) degrees of freedom. This equation gives us a method for simulating the expected statistical filter response for the fields represented by Figure 8a and the GRACE observational system.
[32] To obtain realistic values for the s 'm we used the weighted root-mean-square (WRMS) residuals of our previous fit to the GRACE data. (We can also use the scaled standard deviations available in the CSR RL04 data set without significant change to our results.) Since the expected F value depends on the SNR, the amplitude of our model field must be commensurate with our observed GRACE field to perform an appropriate simulation. The GLDAS soil moisture field we used, however, accounts for only a part of the overall terrestrial water storage budget [e.g., Rodell et al., 2007]. Rather than attempt to model the other processes in detail, we adopted the expedient of multiplying the GLDAS fields by a factor of 1.5 to create model fields that approximated the appropriate amplitude.
[33] Figure 8b shows the simulated statistically filtered (noiseless) hydrology base fields. The filtered fields clearly show the effect of the reduced range of degrees that result from the filtering process. The sharpness of most features is reduced, as is the intensity, but otherwise most of the important features are reproduced. The largest exception is the small-scale negative amplitude feature in the in-phase field just north of the Amazon River mouth. (The locations of the largest differences are noted with a cross.) The statistical filtering field underpredicts the amplitude of this negative feature by $180 mm. For the Gaussian smoothed fields, the greatest difference from the original field is seen in the in-phase component in northern India. There are two spatially small, large-amplitude features that the smoothed fields underpredict by $215 mm. In both the statistical filtering and Gaussian smoothing cases, the largest differences occur when small-scale features are close to features having the opposite sign.
[34] It is useful to note how each of the filters responds to the low-level signal in the northern ocean areas mentioned above. The statistical filter (Figure 8b) tends to break this signal up. These can also be seen in the statistically filtered GRACE fields. The Gaussian filtered field (Figure 8c) preserves much of the intensity while smearing the signal out locally. When destriping is used (Figure 8d), the features stretch out into longitudinal filaments.

Statistical Filter for the Trend
[35] The degree-averaged statistical filter values for the trend (Figure 9) are distinct from those for the annual components ( Figure 4). They share the rejection of coefficients at low degrees, but many more are rejected for ' 10 for the trend. The annual component filter decreases rapidly between degrees 10 and 30, with only a small number of coefficients accepted beyond degree 30. In contrast, the trend filter decreases almost linearly between degrees 5 and 40, a much less rapid decrease than the Gaussian smoothing filter, and tails off slowly beyond degree 40.
[36] The significance of this result can be seen in the reconstructed field. In Figure 10 we show six such reconstructions. We choose to show the trend for the geoid instead of surface-mass-thickness to stress that some of the gravity trend is due to GIA. (Inferences of surfacemass-thickness or deformation assume that the observed gravity changes are associated only with a change of surface mass on an elastic Earth, whereas the gravity changes associated with GIA are due to viscoelastic deformation.) Figure 10a shows, for reference, the undestriped, unfiltered rate field, which is clearly dominated by the striping artifact. The trend patterns in the destriped fields ( Figure 10b) are much clearer, even with no Gaussian smoothing, with the exception of the equatorial region (consistent with the findings of Swenson and Wahr [2006]).
[37] Smoothing either with an isotropic Gaussian filter (Figure 10c) or with a smaller-footprint anisotropic Gaussian filter (Figure 10d) leads to a qualitative improvement in the maps, with the isotropic Gaussian smoothing filter yielding more rounded shapes. Some of the peak trends (in Alaska, Laurentia, southern Greenland, and Antarctica, for example) are reduced with the isotropic Gaussian smoothing filter.
[38] Compared to the negligible effect of destriping for the annual components, the destriping process has an observable effect on the statistical filter for the trend, although the statistical filter by itself removes most of the stripes (Figure 10e). The remnant stripes (Figure 10e) are located at the latitude where the stripes in the unfiltered field are most prominent, and appear to be associated with the (spherical harmonic) degree 60 term. When we previously used the CSR RL01 GRACE fields, which includes degrees 61-70 (and greater for some fields), we obtained a similar result except the ''offending'' components were above degree 60.
[39] All the major features in the Gaussian smoothed geoid trend fields (Figures 10c and 10d) are present in the statistically filtered destriped field (Figure 10f), although the detailed shapes of the features can be significantly different. Some minor features that appear in the Gaussian smoothed trend maps, like the signal in India or in the southeast Atlantic Ocean, are not present in the statistically filtered maps. Furthermore, the positive geoid trends that appear between ±60°longitude in the Antarctic do not appear in the Gaussian smoothed maps. The trend in Europe is much broader in the statistically filtered field.
[40] The statistically filtered trend maps (destriped and undestriped) have a somewhat ''mottled'' appearance in some areas. The reason for this can be seen from Figure 11, which shows the spatial representation of the statistical filters. In addition to the main Gaussian-like central lobe of the filter with elongated east-west axis, there are sidelobes out as far as $15°. As with the annual filter, the first sidelobes on the north-south axis are negative. We do not yet understand the features of the GRACE data sets that give rise to these sidelobes.

Discussion
[41] We have developed a new approach for filtering GRACE gravity data sets. Our statistical filter is distinct from Gaussian smoothing in that it uses the data themselves to test for specific components of the time-varying gravity field. The statistical filter therefore adapts itself to the data set (and may be termed an adaptive filter). The main drawback to the method is that it requires a model for the temporal evolution of the gravity coefficients, and parameters of that model become ''targets'' for recovery by the filter. As the GRACE time series becomes longer, the simple annual-plus-trend model that has commonly been used may require additional complexity. Added complexity would tend to weaken the statistical test on which the filter is based, even as the lengthened time series strengthened it.
[42] The statistical filter appears inherently to remove the stripes, although destriping the fields prior to filtering seems to help the trend recovery due to leakage of the stripe error at the higher degrees. To the extent that the stripes are noise in the temporal domain, the statistical filter should indeed remove the stripes. The stripes could decrease the filter's effectiveness, however, by ''smearing'' the rate signal.
[43] To a great extent, the statistical filter produces maps similar to those of the anisotropic Gaussian smoothing filter because it has a spatial representation qualitatively similar to Gaussian smoothing. This result therefore explains the relative success of those filters. Unlike the Gaussian smoothing filters, however, the statistical filter has significant sidelobes, including large negative sidelobes on the north-south axis. We assume that the details of the statistical filter, including the sidelobes, reveal information on the errors, and the correlations among the errors, for the GRACE coefficients. In this case, it may be useful to compare the statistical filters derived from GRACE coefficients produced by different processing centers.
[44] Compared to the previous CSR release (RL01) of the GRACE fields, the RL04 fields require much less in the way of smoothing or filtering. We used Gaussian smoothing radii that were $50% of the radii that we typically used with the RL01 fields. The destriped rate fields are nearly to the point at which no or little filtering is required. This difference is due mainly to improvements in the analysis (including models), as the data time span is only $10% longer for RL04.
[45] The statistical filter yields features that we assessed using ground-based GPS data in South America by assuming that the annual component of the gravity signal is associated only with groundwater storage [after Davis et al., 2004]. The undestriped, statistically filtered field had a c 2 value relative to the GPS data consistent with the best result from smoothing. We did not attempt to test the trend field, for several reasons: it is not possible to associate the trends with a single process, such as loading, that can relate gravity variations to surface deformation; GPS determinations of the vertical (radial) rate may be highly uncertain due to a range of systematic errors; and models for various relevant processes (GIA, hydrology) are at least as uncertain as the differences between the statistically and Gaussian smoothed trend fields. However, in the future it may be possible to isolate a known signal in a specific area and to perform this comparison.
[46] Acknowledgments. This research was supported by NASA grants NNG04GF09G, NNG04GL69G, and NAG5-13748 and by the Natural Sciences and Engineering Research Council of Canada (J.X.M.). Some figures were generated using Generic Mapping Tools version 4 [Wessel and Smith, 1998]. We appreciate the information provided by C. K. Shum on the anisotropic Gaussian smoothing filter. We thank S. Hagen for his opinion on certain points of grammar. The manuscript benefited from reviews by P. Tregoning, R. Guillaume, and an anonymous referee.