Summary

We present a new non-isotropic Gaussian filter for smoothing mass changes computed from Gravity Recovery and Climate Experiment (GRACE) L2 products and a new method to reduce land–ocean signal leakage caused by Gaussian smoothing. The kernel of our non-isotropic filter is the product of two Gaussian functions with distinct latitudinal and longitudinal smoothing radii. When expressed as number of kilometres at the Earth's surface, the longitudinal smoothing radius, defined as a fixed longitude interval, is longer at the equator, and shorter at higher latitude. This is principally in accordance with the resolution of the GRACE data, and permits us to produce homogeneously smoothed results without excessive smoothing in latitudinal direction. This filter is not applicable in polar regions, where we choose to use the classical isotropic Gaussian filter. A smoothing radius choice scheme is proposed for the two filters to mesh seamlessly. In our leakage reduction method, the inputs are the mass change data after smoothing using a Gaussian filter. Along coasts where mass change signal on land is far larger than that over ocean (or signal over ocean is reduced to a very small magnitude by removing a model beforehand, and adding the model back afterwards), our method approximately recovers a smoothed mass change signal over both land and ocean sides as if a regional Gaussian filter with the same smoothing radius were applied over land and ocean separately, in which no signal leakages appear. The side lobe problem does not appear in our approach. Our leakage reduction method could also be used to study mass changes within a region where signal is far larger than that in surrounding regions, or where signal in surrounding regions could be reduced to very low magnitude by removing a model.

Introduction

The final results of the Gravity Recovery and Climate Experiment (GRACE) satellite observations are customarily provided as monthly gravitational field models in terms of spherical harmonic coefficients (SHCs), known as Level 2 (L2) products, which are widely used to study mass changes at the Earth's surface. However, large errors in the form of heavy south–north stripes are present in mass changes obtained. In practical applications, the stripes are customarily removed by smoothing, that is, the values of mass changes at the Earth's surface are averaged with a certain prescribed weighting function (e.g. Jekeli 1981; Wahr et al. 1998; Swenson & Wahr 2002; Han et al. 2005). Since the discovery that a large part of the stripes is due to correlated errors among the SHCs (Swenson & Wahr 2006), many authors apply a high pass filter to the SHCs for even and odd degrees separately before smoothing for reducing correlated errors in the SHCs (e.g. Chambers 2006; Swenson & Wahr 2006; Chen et al. 2007b; Duan et al. 2009), that we refer to as decorrelation here.

In this work, we follow the decorrelation-smoothing approach to compute mass changes from GRACE L2 data, focussing on the smoothing step. However, other methodologies have been developed: Kusche (2007) and Klees et al. (2008) made use of GRACE solution error covariances and a priori signal covariances to design optimal filter. Schrama et al. (2007) and Wouters & Schrama (2007) used empirical orthogonal function approaches. Davis et al. (2008) presented a statistical approach to extract trend and annual signals. A summary in more detail could be found in Duan et al. (2009).

Smoothing is customarily performed using an isotropic Gaussian filter (e.g. Jekeli 1981; Wahr et al. 1998), that is, when computing the smoothed value of mass change at a certain point, the Gaussian function depending only on the distance to and having maximum at that point is chosen as the averaging weighting function, often called smoothing kernel. A smoothing radius—the distance to that point where the value of the Gaussian function is half of its maximum—should be prescribed empirically (Chen et al. 2005; Seo & Wilson 2005). Although the isotropic Gaussian filter is applied in the spatial domain consisting of the values of mass changes at the Earth's surface, it is quite simple to do the computation in the spectral domain consisting of the values of the SHCs computed from the mass changes: each of the SHCs Clm and Slm is multiplied with a degree-only dependent factor Wl computed from the smoothing radius. The values of Wl decrease as l increases. Based on the fact that the SHCs have larger errors when the order m is larger, Han et al. (2005) generalized the isotropic Gaussian filter in the spectral domain by replacing the degree-only dependent factor Wl with a degree-and-order dependent factor Wlm. The values of Wlm depend on l in the same form as Wl, but computed using larger smoothing radius for larger m up to a certain limit m1. Hence this generalization smoothes different order spherical harmonic expansion terms with isotropic Gaussian filters of different smoothing radii. The authors referred to this generalization as a non-isotropic filter. It requires two prescribed smoothing radii: One for m= 0 that the authors referred to as south–north smoothing radius, and one for m=m1 that the authors referred to as east–west smoothing radius. The smoothing radii for 0 < m < m1 are obtained by linear interpolation. Isotropic smoothing algorithms have also been developed for series expansions using other base functions (Fengler et al. 2006; Schmidt et al. 2006). Furthermore, non-Gaussian filters have also been applied (Wu et al. 2009).

The non-isotropic filter we present in this work is a different generalization of the isotropic Gaussian filter. We recognize the fact that mass changes obtained from GRACE data have better precision at higher latitude. This may be due to the fact that GRACE satellites have nearly polar orbits, that is, the numbers of ground tracks across different parallels are equal except very close to the poles, resulting in denser ground tracks (number of track per kilometre along parallels) at higher latitudes. Due to this characteristic of orbits, the east–west resolution wave length of the GRACE measurements may be better represented using a longitudinal interval rather than a latitude-independent length (i.e. number of kilometres) at the Earth's surface (the same longitudinal interval Δλ corresponds to a surface length of R sin θΔλ that is shorter at higher latitude, where R is the radius of the Earth, and θ is the colatitude). In our non-isotropic Gaussian filter, the weighting function is the product of two Gaussian functions depending on the colatitude θ and the longitude λ, respectively, with different prescribed empirical smoothing radii. However, in polar regions, our non-isotropic filter is not applicable as the longitudinal smoothing radius expressed in terms of length at the Earth's surface becomes too small, and thus an isotropic Gaussian filter that meshes seamlessly with our non-isotropic filter is applied.

We also mention that, for explaining our leakage reduction method, we will devise computation algorithms for the isotropic and our non-isotropic Gaussian filters based on data over a latitude–longitude grid. These algorithms may be easily modified to freely define new filters with position-dependent smoothing radius or radii.

The smoothed value of mass changes at a point is a weighted average of the unsmoothed values around that point. Hence, at a point in coastal region, whether over land or ocean, unsmoothed values over both land and ocean are used to compute the smoothed value. Supposing the signal of mass changes on land is larger than that over ocean, smoothing would then attenuate the amplitude on land (e.g. Chen et al. 2005, 2007a). This attenuation represents a loss of a part of the signal on land that is added to the smoothed value over ocean. In other words, we say that the signal on land is leaked into the coastal ocean. Wahr et al. (1998) presented an iterative algorithm to alleviate the effect of the signal leakage. The essence of their algorithm to compute the smoothed mass changes over ocean is (1) filtering the GRACE results with a small but large enough smoothing radius so that most of errors/stripes are removed, (2) removing the smoothed mass changes on land from the unsmoothed GRACE results and then (3) filtering to obtain the smoothed mass changes over ocean. The algorithm for computing the smoothed mass changes on land is similar. Some authors used an empirical amplitude scaling factor to correct the leakage while estimating regional mass changes (e.g. Fenoglio-Marc et al. 2006; Velicogna & Wahr 2006; Baur et al. 2009). Swenson & Wahr (2002) devised two smoothing algorithms in addition to the Gaussian filtering for estimating regional mass changes that may suffer from signal leakage. One algorithm is to minimize the sum of satellite data and signal leakage errors. The other is to minimize the leakage error while preserving a prescribed satellite data error. Chen et al. (2006) used a trial-and-error forward modelling approach to avoid leakage effect in their estimate of glacier melt.

In this work, we present an algorithm to reduce signal leakage by, loosely speaking, reverting the Gaussian filter. A flowchart of the method could be found in fig. 2 of Shum et al. (2010). The input is the mass change signal computed from GRACE L2 data after decorrelation-smoothing. The choice of decorelation filter and smoothing radius is actually empirical, for which a guideline is provided by Duan et al. (2009). The algorithm is based on the assumption that errors causing stripes are considerably removed by decorrelation-smoothing, so that the mass change signal obtained has little error, and can be considered practically the same as that obtained by smoothing purely the mass change signal alone. Assuming the mass change signal on land is far larger than that over coastal ocean (or signal over ocean is reduced to very small magnitude by removing a ocean model from GRACE data before smoothing), so that the contribution of the signal over ocean can be neglected in the results of the smoothed mass changes on land, the leakage-reduced smoothed mass changes on land are computed by re-scaling the output of the Gaussian filter so that the final results are equal to that obtained by a weighted averaging of the land signal alone. The result of this approach is in fact equivalent to that of a regional filter for land alone while avoiding the side lobes that Jekeli (1981) and Wahr et al. (1998) mentioned. To recover the leakage-reduced smoothed mass changes over ocean, we construct a case of hypothetical global mass changes which are equal to the leakage-reduced smoothed mass changes on land, and vanish over ocean. We then smooth it using the same Gaussian filter. The results over ocean then represent the leakage from land caused by Gaussian smoothing. Hence, the leakage-reduced smoothed mass changes over ocean are obtained by subtracting the smoothed results of the hypothetical case from the smoothed GRACE results. If the mass change signal is strictly zero over ocean, our algorithm removes the leakage completely on land, while the iterative algorithm of Wahr et al. (1998) does not. However, the assumption that the mass change signal on land is far larger than that over ocean is not true everywhere. Hence, the algorithm should be applied on a regional basis, that is, every coast should be carefully examined if the assumption is satisfied, or could be made satisfied by removing a model of mass changes over ocean (model should be added back afterwards). This algorithm could also be used to study regional mass change signal like in the case of Fenoglio-Marc et al. (2006) who studied mass changes over Mediterranean Sea by first removing mass change signal on surrounding land using data from a hydrological model.

This paper is arranged as follows. In Section 2, we briefly summarize the classical Gaussian filter as reference. In Section 3, we define our non-isotropic filter, and derive more efficient computation formulae based on a cos-Fourier expansion (Guo & Shum 2009). In Section 4, we define our non-isotropic and the classical isotropic Gaussian filters for smoothing the mass changes over land and ocean separately as background for the leakage reduction method presented in Section 5. The results are summarized in the Section 6.

The Isotropic Gaussian Filters

We begin by summarizing the essence of the classical isotropic Gaussian filter as it serves as a background for our non-isotropic Gaussian filter to be defined in the next section. Furthermore, we will also deal with signal leakage reduction caused by the classical isotropic Gaussian filtering.

In this work, we write the Gaussian function in the form of:  

1
formula
where the standard deviation σx is related to the smoothing radius rx customarily used in literature (Jekeli 1981; Wahr et al. 1998), that is, the value of x where G(x) =G(0)/2, in the form of  
2
formula
We have used the subscript x in σx and rx to specify that the variable is x. In the following, x will be replaced by the geocentric angle ψ defined below, and accordingly, σx and rx will be replaced by σψ and rψ, respectively.

We assume the Earth to be a sphere of radius R, and use the colatitude θ and longitude λ to represent positions at the Earth's surface. We use f(θ, λ) to denote the mass changes at the Earth's surface inferred from GRACE observations, which exhibit stripes, and are to be smoothed. For the classical isotropic Gaussian filter, the smoothed value of f(θ, λ) at the point graphic is defined as  

3
formula
where ψ is the geocentric angle between the points (θ, λ) and graphic given by the equation  
4
formula
and  
5
formula
that is in fact independent of graphic.

The smoothing radius rψ is an angular distance on sphere. In the literature, it is most often expressed in terms of length at the Earth's surface, rd, that is related to rψ in the form of  

6
formula

Practical computation formulae using spherical harmonic expansion were developed by Jekeli (1981) who obtained the factor Wl mentioned in the Introduction, and not to be reproduced here. The result was cited by many others (e.g. Wahr et al. 1998; Han et al. 2005). The formulae (3)–(5) will be used later in this work to derive an computation algorithm using the values of f(θ, λ) over a latitude–longitude grid.

A New Non-Isotropic Gaussian Filter

In the isotropic Gaussian filter defined in (3)–(5), the smoothing kernel is G(ψ), and the curve defined by kernel = 1/2 is a circle of radius rψ on the unit sphere centred on the point graphic where the filtered value is computed. In the non-isotropic filter we will define below, the curve defined by kernel = 1/2 will be an ellipse with distinct semi-axes, or radii, rθ and rλ along latitude and longitude, respectively.

In our non-isotropic Gaussian filter, the smoothing radii along latitude and longitude are based on θ and λ, and are expressed in the form of rθ=δθ and rλ=δλ, where δθ and δλ may be different. The smoothing radius in terms of number of kilometres at the the Earth's surface along latitude is then Rrθ, and that along longitude is R sin θrλ. Our non-isotropic filter is different from that of Han et al. (2005) in two ways: (1) we define explicitly our filter in the spatial domain; (2) our longitudinal smoothing radius expressed in terms of length at the Earth's surface, R sin θrλ, varies as function of latitude. For example, the value of R sin θrλ at latitude |B| = 60° is half of its value at the equator.

Like the isotropic Gaussian filter, our non-isotropic Gaussian filter is also based on the Gaussian function defined in (1) and (2). In the following text, we will replace x by graphic and graphic by σθ and σλ, and rx by rθ and rλ to define our non-isotropic filter. The meaning of these symbols are self-explanatory.

For our non-isotropic Gaussian filter, we defined the smoothed value of f(θ, λ) at the point graphic to be  

7
formula
 
8
formula
where graphic is the smoothing kernel, and the function W(θ) is defined as  
9
formula
For evident reason, it suffices to define the domain of integration to be graphic for λ, and graphic for θ. We have extended the domain of integration to the whole surface of the Earth for deriving an easy computation algorithm. It can be inferred from the formulae of definition that the poles are singular points, as the domain of integration for θ cannot be extended beyond the poles. However, in polar regions, because the smoothing radius along longitude is too short when expressed as length at the Earth's surface, the non-isotropic filter thus defined is not appropriate as will be seen from numerical examples, and the classical isotropic Gaussian filter will be applied. Hence this singularity at the poles is not a problem.

The definition formulae of the filter, (7)–(9), should be further developed for practical computation. Unfortunately, we do not see any possibility to use the spherical harmonics as in the case of isotropic Gaussian filter and its variant of Han et al. (2005).

Here, to derive a fast computation algorithm of the filter defined in (7)–(9), we will use the cos-Fourier expansion (Canuto et al. 1987; Guo & Shum 2009) summarized in the Appendix to represent f(θ, λ). The values of f(θ, λ) at the gridpoints (θj, λk), which are considered as the grid cell centres of a constant interval regular latitude–longitude grid, are first computed according to the formulae given by Wahr et al. (1998) using the GRACE determined geopotential coefficients, that is, the SHCs of the L2 products. The cos-Fourier expansion coefficients are then computed using the values at the gridpoints. According to (A2), we express (7) and (8) in the following forms for further treatment  

10
formula
where  
11
formula
 
12
formula
 
13
formula
 
14
formula

In (11) and (12), the practical value of σλ is very small compared to π, so that graphic can be regarded as zero for graphic, and extending the domain of integration of λ from graphic to [−∞, ∞] introduces only a negligible error. With this approximation we obtain, after some mathematical manipulation,  

15
formula
 
16
formula

To evaluate the integrals in (13) and (14), we express graphic in the form of a cosine expansion (Guo & Shum 2009; Canuto et al. 1987)  

17
formula
where the coefficients graphic can be computed using the values of graphic at the gridpoints of θ defined in (A7) and reproduced in (21):  
18
formula
where c0= 2, and cn= 1 if 0 < nN. The formulae (17) and (18) are a pair discrete Chebychev transform and its inverse, that are similar to the pair of (A3) and (A8).

Substituting (17) into (13) and (14), we obtain  

19
formula
where we used the relation  
20
formula

It is straightforward to write out the formulae for computing the filtered values at the gridpoints given in (A5) and (A7), that we rewrite here using grid intervals:  

21
formula
 
22
formula
where the formula for λk is rewritten to include the case when the number of nodes is even that is discussed in the last paragraph of the Appendix. In fact, (21) has been used in (18), and both (21) and (22) will be used later in the next section.

As already mentioned, the filter thus designed is not applicable in polar regions where we choose to use the classical isotropic Gaussian filter. For GRACE data presently available, we always assign a larger value to rλ than to rθ. At a certain critical latitude |Bc|, the smoothing radii in terms of length at the Earth's surface along latitude, Rrθ, and that along longitude, Rrλ cos Bc, are equal, and the non-isotropic filter become isotropic. Thus, we have cos Bc=rθ/rλ. The classical isotropic Gaussian filter is then applied for |B| > |Bc|. For example, if we assign rλ= 6°, rθ= 3°, we have |Bc| = 60°. In this way, the results of the two filters mesh seamlessly at |Bc|.

In this work, we choose to use the JPL GRACE Level 2 (L2) product RL04 truncated at degree and order 70 for illustrating our method, as higher degee/order coefficients have little influence on the results after filtering. As our purpose here is limited to illustrate the characteristic of our filter (and leakage reduction method later), we choose not to apply any postglacial rebound (PGR) correction, and skip the decorrelation step. In our numerical example, the trend obtained by fitting the geoidal height changes with an offset, a linear term that represents the trend, an annual term, a semi-annual term and a biennial term is used. The mass changes are computed from the geoidal height changes using the load Love numbers that Guo et al. (2004) computed by solving the differential equations for each degree.

Some results are shown in Fig. 1 with W(θ) = sin θ. In fact, the results obtained with W(θ) = sin θ or W(θ) = 1 do not show visible difference. We have chosen rλ= 6° and rθ= 3°, hence |Bc| = 60°. We see that the results of our non-isotropic filter are not appropriate in polar regions where the isotropic Gaussian filter is applied. However, in the latitude range |B| < =60°, our non-isotropic filter provides far smoother results than the classical Gaussian filter. This is attributed to our longer smoothing radius along longitude alone. The longitudinal smoothing radius of our non-isotropic filter, when expressed in terms of distance at the Earth's surface, is R· sin 90°· [π· (6/180)] = 668 km at the equator, and decreases to R· sin 30°· [π· (6/180)] = 334 km at |B| = 60°, while the latitudinal smoothing radius R· [π· (3/180)] = 334 km being the same everywhere. As compared to the isotropic filter, the choice of the non-isotropic filter permits us to obtain smoother results without excessive filtering in the latitudinal direction. The mixed filtering takes advantage of both of the classical isotropic and our non-isotropic Gaussian filters.

Figure 1

Trends of geoidal height changes (left-hand panel) and surface mass changes (right-hand panel) computed from the JPL L2 product RL04 smoothed using different filters. For our illustration purpose, postglacial rebound correction is not made. The first row shows the result without any smoothing applied. The second row shows the results of Gaussian filter with a smoothing radius corresponding to 3° at the equator: R· [π· (3/180)] = 334 km. The third row shows the results using our non-isotropic filter with rθ= 3° and rλ= 6°. The last row shows the results of mixed filter: our non-isotropic filter is used in |B| < =60°, and the isotropic Gaussian filter is used in |B| > 60°. For the values of rθ and rλ chosen, our non-isotropic filter is isotropic at |B| = 60°. Hence the results of the two different filters mesh seamlessly.

Figure 1

Trends of geoidal height changes (left-hand panel) and surface mass changes (right-hand panel) computed from the JPL L2 product RL04 smoothed using different filters. For our illustration purpose, postglacial rebound correction is not made. The first row shows the result without any smoothing applied. The second row shows the results of Gaussian filter with a smoothing radius corresponding to 3° at the equator: R· [π· (3/180)] = 334 km. The third row shows the results using our non-isotropic filter with rθ= 3° and rλ= 6°. The last row shows the results of mixed filter: our non-isotropic filter is used in |B| < =60°, and the isotropic Gaussian filter is used in |B| > 60°. For the values of rθ and rλ chosen, our non-isotropic filter is isotropic at |B| = 60°. Hence the results of the two different filters mesh seamlessly.

A Naive Filter Distinguishing Land and Ocean

Suppose we have data of mass changes all over the Earth's surface on a latitude–longitude grid. Let's say, the mass changes are related to water storage including underground water, snow and ice, etc. on land, and to the mass redistribution over ocean. We suppose that these data do not have significant errors as the GRACE observed mass changes do, and we still want to smooth the data using a Gaussian filter. However, as the data represent completely different phenomena over land and ocean, we do not want to merge the data over land and ocean to produce leakage as what a usual Gaussian filter may result in.

Here we define a filter for such a purpose: only the data on land are used to computed the smoothed values at the gridpoints on land, and only the data over ocean are used to computed the smoothed values at the gridpoints over ocean. Hence, the filter does not produce leakage between land and ocean. However, it can be imagined that at a gridpoint very close to a coast, not all data within a certain distance to the gridpoint are used to compute the smoothed value as in the case of usual Gaussian filters. Instead, only the data at the land or ocean side are used, depending on which side the gridpoint is located. This may be understood as a sacrifice of the quality of the smoothed data in favour of avoiding leakage. Due to this fact, if the original data have large errors like the stripes in the GRACE observed mass changes, or if the data have been expressed as spherical harmonics series which include truncation errors, side lobes along coasts at both the land and ocean sides may appear in the smoothed data (Jekeli 1981; Wahr et al. 1998). For this reason, we call this filter ‘naive’. Nevertheless, this filter represents a compromise between the quality of the smoothed data and the removal of leakage if the original data do not have significant errors. We will use this filter as a reference in the next section where a method for reducing leakage caused by usual Gaussian filters, for example, the classical isotropic Gaussian filter and the non-isotropic one presented in the last section, will be presented.

We first modify the non-isotropic Gaussian filter presented in the last section to filter land and ocean separately. We then modify the classical isotropic Gaussian filter, that is quite straightforward.

We define graphic and graphic as the ocean area and land area, respectively, in the rectangle graphic and graphic. For mass changes f(θ, λ) at the Earth's surface, we modify (7) and (8) to define its smoothed value at the point graphic to be  

23
formula
 
24
formula
where  
25
formula
and the Gaussian function G and the function W are defined in (1) and (9), respectively.

The parameters σθ, σλ and a, b should be prescribed. As already mentioned, the smoothing radius is often used in the literature. In our case, the smoothing radii on latitude rθ and on longitude rλ are different. They are related to σθ and σλ, respectively, in the same form as (2). For evident reason, we define a and b according to  

26
formula

In practice, computations are done using values of f(θ, λ) over a regular latitude–longitude grid. Let's suppose the grid intervals are Δθ and Δλ, respectively, and the gridpoints (θj, λk), which are understood as the grid cell centres where the values of fj, λk) are provided, are defined in (21) and (22).

We do not further divide a grid cell into ocean and land, i.e. whole of the cell is considered either ocean or land. Filtered values are to be computed for all the gridpoints: graphic defined in (21) and (22). Fig. 2 shows a portion of a grid, where thin solid lines separate grid cells, and large dots mark the gridpoints (θj, λk) located at grid cell centres. Land grid cells are grey, and ocean grid cells are white.

Figure 2

A portion of a grid used in the evaluation of the integrals of the non-isotropic Gaussian filter separating land and ocean. The original grid cells are separated by thin solid lines, and their centres are marked with large dots. The land is marked with gray. The domains of integrations of the convolution for different points, the circled large dots, are marked using thick solid lines. The subdivided grid cells for evaluating the integrals are separated by dashed and thin solid lines, and their centres are marked with small and large dots.

Figure 2

A portion of a grid used in the evaluation of the integrals of the non-isotropic Gaussian filter separating land and ocean. The original grid cells are separated by thin solid lines, and their centres are marked with large dots. The land is marked with gray. The domains of integrations of the convolution for different points, the circled large dots, are marked using thick solid lines. The subdivided grid cells for evaluating the integrals are separated by dashed and thin solid lines, and their centres are marked with small and large dots.

For more accurately evaluating the integrals in the formulae of the filter, we subdivide the grid to grid intervals Δθ/μ and Δλ/μ, where μ is an integer which we choose to be odd, so that the original grid cell centres are still the centres of certain subdivided grid cells. We show an example for μ= 3 in Fig. 2, where dashed and solid thin lines separate subdivided grid cells, and small and large dots mark the centres of these smaller cells. The coordinates of the subdivided grid cell centres are given by  

27
formula
The values of fsα, λsβ) are obtained by interpolation from the values of fj, λk) (Guo & Shum 2009).

In discrete case, we round a and b in the form of  

28
formula
where Δθ/μ and Δλ/μ are grid intervals of the subdivided grid, and m and n are integers, which we choose to be odd so that the grid cell centre of interest, graphic, could be most likely at the centre of integral domain in (23) and (24). In fact, graphic could be truly at the centre of integral domain only if the rectangle graphic and graphic is purely ocean or land. We show in the Fig. 2 some integral domains of the convolution in (23) and (24) using thick solid lines assuming m=n= 3μ= 9, where circled large dots inside the domains are the grid cell centres where the filtered values are to be computed.

Different continents and islands may be considered separately in the filter. Here we refine the definition of graphic such that it includes only the continent or island where graphic is located. Our guidelines in dividing continents and islands are: (1) two land grid cells having a common side or corner belong to the same continent or island; (2) two land grid cells separated by ocean grid cells belong to different continents or islands.

Based on the subdivided grid, we can write the filtered values for all the original grid cell centres defined in (21) and (22) according to (23) and (24) in the form of summations:  

29
formula
 
30
formula
where S is the domain of integration as illustrated in Fig. 2. A common factor (Δθ/μ)(Δλ/μ) is omitted in both formulae, that does not influence the final results of fsj, λk).

Similar to the case without distinguishing land and ocean, our non-isotropic filter is not applicable in polar regions, where we will use a modified version of the isotropic Gaussian filter. With the grids defined for our non-isotropic Gaussian filter, it is straightforward to write out the formulae of the isotropic Gaussian filter that distinguishes land and ocean. When evaluating the integrals in (3) and (5) to compute the smoothed value at the grid cell centre (θj, λk), we need to include only the subdivided grid cells within a distance of 3σψ that belong to the same continent, island or ocean as (θj, λk), which we denote using L and O, respectively. We then have  

31
formula
 
32
formula
where ψαβjk is computed using (4) with graphic and (θ, λ) = (θsα, λsβ), and  
33
formula

If land and ocean are not distinguished (just consider the Earth's whole surface as ocean in the formulae), the filters defined in this section give the same results as the non-isotropic filter defined in the last section and the classical isotropic Gaussian filter, respectively, though requiring far more computer time. However, these discrete forms of the filters may be easily modified to freely define new filters with position-dependent smoothing radius or radii.

Filtering that distinguishes land and ocean may not be meaningful for the geoidal height changes that is continuous everywhere. In Fig. 3, we present the filtered results of the trends of mass changes computed from GLDAS data (we used the daily data provided by the Special Bureau for Hydrology of the IERS The Global Geophysical Fluids Center* to compute the monthly means) and the JPL GRACE L2 product RL04. The GLDAS trend is computed in a similar way as the GRACE trend, i.e. the GLDAS data are fitted with an offset, a linear term that represent the trend, an annual term, a semi-annual term and a biennial term. We see that there are side lobes along coasts in the case of GRACE results (we have chosen the JPL GRACE L2 product to truncate at degree 70 for the side lobes be very apparent). These side lobes may be eliminated by filtering the ocean function in a similar way as what Swenson & Wahr (2002) did for the regional margin for inferring regional mass changes from GRACE data. However, this may introduce leakages that need to be further treated. Instead of seeking to avoid the side lobes, in the next section, we will use the filters distinguishing land and ocean as references to reduce the leakages caused by the filters without distinguishing land and ocean.

Figure 3

Trends of mass changes at the Earth's surface computed from the GLDAS (left-hand panel) and the JPL GRACE L2 product RL04 (right-hand panel) smoothed using different algorithms distinguishing land and ocean. The first row is the result without any smoothing applied. The second row is the results of the isotropic Gaussian filter with a smoothing radius corresponding to 3° at the equator: R*3π/180 = 334 km. The third row is the results using our non-isotropic Gaussian filter with rθ= 3° and rλ= 6°. The last row is the results of mixed filter: our non-isotropic Gaussian filter is used in |B| < =60, and the isotropic Gaussian filter is used in |B| > 60. For the values of rθ and rλ chosen, our non-isotropic Gaussian filter is isotropic at |B| = 60°. Hence the results of the two different filters mesh seamlessly. We see the appearance of side lobes along coasts on the right-hand panel.

Figure 3

Trends of mass changes at the Earth's surface computed from the GLDAS (left-hand panel) and the JPL GRACE L2 product RL04 (right-hand panel) smoothed using different algorithms distinguishing land and ocean. The first row is the result without any smoothing applied. The second row is the results of the isotropic Gaussian filter with a smoothing radius corresponding to 3° at the equator: R*3π/180 = 334 km. The third row is the results using our non-isotropic Gaussian filter with rθ= 3° and rλ= 6°. The last row is the results of mixed filter: our non-isotropic Gaussian filter is used in |B| < =60, and the isotropic Gaussian filter is used in |B| > 60. For the values of rθ and rλ chosen, our non-isotropic Gaussian filter is isotropic at |B| = 60°. Hence the results of the two different filters mesh seamlessly. We see the appearance of side lobes along coasts on the right-hand panel.

Finally, we remark that GLDAS trend is reasonably smoothed without producing leakage.

A Leakage Reduction Method

We have defined two types of filters. One type does not distinguish land and ocean that can remove stripes in GRACE determined mass changes, but causes leakage between land and ocean. The other type distinguishes land and ocean, that does not cause leakage between land and ocean by definition, but creates side lobes along coasts when applied to GRACE determined mass changes. In this section, we present an algorithm to reduce leakage caused by filtering without distinguishing land and ocean.

For better explaining our algorithm, we divide the GRACE determined mass changes into signal and error parts  

34
formula
where the GRACE error eg is responsible for stripes. We use graphic (where graphic stands for ‘global’) to denote a filter without distinguishing land and ocean, and graphic, the corresponding filter distinguishing land and ocean. The purpose of filtering is to reduce the influence of the error term eg. After smoothing without distinguishing land and ocean, the error term eg is very much reduced, and we have  
35
formula
When the smoothing radius is large enough, we assume that the ‘smoothed error’esg is negligibly small, that is, the filter graphic removes stripes with only negligible error left, and produces a map of smoothed mass changes with signal leakage between land and ocean. Our ideal map of smoothed mass changes without leakage between land and ocean we expect to obtain is graphic. Our objective of this section is to find a way to recover graphic from graphic, which is assumed equivalent to graphic as expressed in (35) where esg is assumed negligible.

Before introducing our algorithm, we summarize here an algorithm for reducing leakage proposed by Wahr et al. (1998). Their procedure for estimating oceanic mass changes with reduced leakage is:

  • (i)

    Filter the global GRACE determined mass changes with a small but yet large enough smoothing radius so that most of the errors/stripes are removed, and then set the oceanic mass changes to zero. As result, the mass changes are equal to that on land, and zero over ocean. Here we call this result smoothed mass changes on land.

  • (ii)

    Compute the geopotential coefficients corresponding to the smoothed mass changes on land obtained in step (i), and then subtract the results from the original GRACE determined geopotential coefficients to obtain the geopotential coefficients corresponding to oceanic mass changes.

  • (iii)

    Compute the oceanic mass changes using the geopotential coefficients corresponding to oceanic mass changes obtained in step (ii) by smoothing appropriately using a filter.

Now we suppose that the mass changes on land are far larger than that over ocean. We can then see in step (i) that, along coasts, the smoothing would attenuate the amplitude of mass changes on land, accompanied by a leakage of signal from land into ocean. As those are the attenuated mass changes on land obtained in step (i) which are removed in step (ii), the result of oceanic mass changes obtained in step (iii) should still contain leakage from land, though reduced indeed. Our algorithm described bellow aims at further improving the reduction of leakage.

The leakage reduction method of Wahr et al. (1998) can be computed using spherical harmonics for saving computer time. However we don't see any similar way to reduce computer time for our method. So our algorithm will be based on the data over a latitude–longitude grid, and require far more computer time that is no longer a problem nowadays.

Based on our assumption that the error term in (35) could be neglected, graphic can be considered equivalent to graphic, or the stripes can be considered removed in the smoothed data after filtering. Hence there is no need to worry about the stripes when using the smoothed data to reduce the leakage. We will only refer to graphic hereafter, keeping in mind that it is actually obtained by computing graphic. Here we state our objective more concretely: to recover graphic from graphic under the assumption that the mass changes over land are far larger than that over ocean along coasts. A flowchart of the method was given in fig. 2 of Shum et al. (2010).

We will use similar forms of the discrete isotropic Gaussian filter defined in (31) and (32) to illustrate our algorithm, that can be easily extended to the case of our non-isotropic Gaussian filter. The advantage of the discrete forms is that they may be directly implemented for computation. We will first recover graphic on land, that will then be applied to recover graphic over ocean.

First, let's examine how the signal on land is attenuated in the classical isotropic Gaussian filter. For this purpose, we need to modify (31) and (32) to get the discrete form of the isotropic Gaussian filter without distinguishing land and ocean, that is, graphic. This is quite straightforward. We use sglbj, λk) and sdloj, λk) to express the values of graphic and graphic at the gridpoint (θj, λk). We can express the result as  

36
formula
 
37
formula
where we used s to replace f to indicate that we are dealing purely with signal without any error to cause the stripes. For understanding how signal on land is attenuated by the filter, we consider a fictitious simple situation: the signal is constant on land and vanishes over ocean, and the point (θj, λk) is right at the land side of the coastline. We see immediately that what we get from the above formulae is about half of the signal amplitude at the point (θj, λk).

For recovering graphic from graphic, we suppose (θj, λk) is on land, and rewrite the formulae of the filter distinguishing land and ocean, that is, graphic, in the form of  

38
formula
 
39
formula
As we supposed that the signal on land is far larger than that over ocean, we may neglect the oceanic contribution to the summation in (36), that is,  
40
formula
Hence, we obtain from (36), (38) and (40) the desired smoothed mass changes over land with reduced leakage  
41
formula
The error due to leakage from ocean, eol, represents the error in the approximation (40). It can be neglected in practical computation based on our assumption made above. Due to (35) where the error term is assumed negligible, we can replace sglbj, λk) by the smoothed result of GRACE data, fglbj, λk), to compute sdloj, λk).

The error eol may be obtained by comparing (36) and (38) to (41) 

42
formula
This is the error due to non-zero signal over ocean (or remanent signal if signal over ocean was removed based on available data). The relative error is  
43
formula
The order of magnitude of the relative error at a gridpoint is approximately the product of two factor: (1) the ratio between the amplitudes of signals over ocean and land, (2) the ratio between the areas of ocean and land within the distance of smoothing radius.

Now we determined the smoothed signal over ocean with reduced leakage. Our idea is to first determine the signal leaked into ocean from land, and then remove the leakage from the smoothed signal over ocean. For this purpose, we construct the following hypothetical mass changes  

44
formula
We then smooth sL using the isotropic Gaussian filter without distinguishing land and ocean, and use sLglbj, λk) to denote the result at the point (θj, λk). The values of sLglbj, λk) over ocean are then the leakage from land. The smoothed signal over ocean after removing leakage, graphic, is then sglbj, λk) −sLglbj, λk). Here we mention again that sglbj, λk) is to be replaced by fglbj, λk) in practical computation.

While we used the discrete form of the filters to demonstrate the principle, the formulae in spectral domain presented in Jekeli (1981) and Section 3 of this paper may be used to save computer time whenever applying the filter graphic, that is, filtering without distinguishing land and ocean.

For the quite idealized hypothetical special case when the signal over ocean vanishes, our algorithm removes leakage completely on land, while the original algorithm presented in Wahr et al. (1998) does not. This can be easily inferred from the characteristics of the methods. However, as the signal on land is smoothed twice when determining the signal leaked into ocean, the leakage over ocean is not completely removed.

For easy reference, we summarize here the procedure of computation:

  • (i)

    Step 1: Smooth the GRACE data using a Gaussian type filter without distinguishing land and ocean.

  • (ii)

    Step 2: Rescale the smoothed data on land using the formula (41) neglecting the error term.

  • (iii)

    Step 3: Create a hypothetical model that has the same data as the result of Step 2 on land, and zero over ocean, and then smooth this model in the same way as STEP 1.

  • (iv)

    Step 4: Obtain the final result that has the same data as the result of Step 2 on land, and the difference between the results of Steps 1 and 4 over ocean.

As numerical examples, we still used the trend of mass changes determined from GLDAS and GRACE L2 product as in Fig. 3 in the last section. In our computation, we assume that the mass change signal on land is everywhere far larger than oceanic signal. This is certainly not true. However, as leakage is local, regions where the assumption may be considered true are not influenced by other regions. Let us look at the South America, particularly the Patagognia region, and the neighboring ocean area where our assumption is believed approximately true for the trend. We see that our leakage reduction algorithm intensifies the signal on land, and drastically removes the leakage over ocean.

For the GLDAS data, the result after reducing leakage shown in Fig. 4 (second row left-hand side) is identical to that of filtering distinguishing land and ocean shown in Fig. 3 (last row left-hand side) over land.

Figure 4

Trends of mass changes at the Earth's surface computed from the GLDAS (left-hand panel) and the JPL GRACE L2 product RL04 (right-hand panel). The first row presents the results of our non-isotropic Gaussian filter (rθ= 3°, rλ= 6°) without distinguishing land and ocean, thus with leakage present. The second row presents the results after removing leakage according to our algorithm. The last row presents the difference between the two rows above, that is the leakage of the filter without distinguishing land and ocean. We see that the removal of leakage only affects the nearby area along coastal lines.

Figure 4

Trends of mass changes at the Earth's surface computed from the GLDAS (left-hand panel) and the JPL GRACE L2 product RL04 (right-hand panel). The first row presents the results of our non-isotropic Gaussian filter (rθ= 3°, rλ= 6°) without distinguishing land and ocean, thus with leakage present. The second row presents the results after removing leakage according to our algorithm. The last row presents the difference between the two rows above, that is the leakage of the filter without distinguishing land and ocean. We see that the removal of leakage only affects the nearby area along coastal lines.

As already mentioned in the last paragraph, the assumption that signal on land is far larger than that over ocean is not true along every fraction of coasts in the world. Hence, our algorithm of reducing leakage could be applied only to specific regions like South America, Greenland and Antarctica where the assumption is supposed to be true. Similar algorithm may be developed if signal over ocean is far larger than that on land, such as the case of Fenoglio-Marc et al. (2006) who studied the mass changes over the Mediterranean Sea by removing the signal on surrounding land beforehand using hydrological data. To remove leakage in the worldwide map of mass changes determined from GRACE observations, every fraction of coasts should be considered specifically. Special care should be especially taken when applying the algorithm to small islands surrounded by ocean where errors may be quite large as shown in (43).

Concluding Remarks

We presented a new non-isotropic Gaussian filter and a new leakage reduction method for determining mass changes using the GRACE L2 products. These methods were proven with real observation data to be valid alternatives to the methods already available. In practical applications, most people would apply a Swenson & Wahr (2006) type method to remove correlated errors in the SHCs of the GRACE L2 product before using them to compute mass changes. However, as this algorithm does not remove all errors in the SHCs, the mass changes thus computed still need to be smoothed. Our non-isotropic filter is an alternative choice to smooth the mass change at this step besides the classical isotropic Gaussian filter and its variants like the one of Han et al. (2005). Our leakage reduction method is an alternative to the one described by Wahr et al. (1998).

In our non-isotropic filter, we defined the kernel, i.e. the averaging weighting function, as the product of two Gaussian functions depending on latitude and longitude, respectively, with distinct prescribed smoothing radii in terms of intervals of latitude and longitude. In this way, the smoothing radius in the longitudinal direction, that is equal to a fixed longitude interval, is shorter at higher latitude when expressed in terms of length at the Earth's surface. This feature makes the filter particularly suitable for smoothing GRACE data, which has better longitudinal resolution (in terms of distance at the Earth's surface) at higher latitude. However, at the polar regions, as our non-isotropic filter is not applicable, the classical isotropic Gaussian filter is applied. Hence, the filter finally used is a combination of two: Our non-isotropic filter when |B| < =Bc, and the classical isotropic Gaussian filter when |B| > Bc. The smoothing radii of the filters are suitably chosen so that the results of the two filters mesh seamlessly at |B| =Bc.

We tested our filter using the trends of geoidal height changes and mass changes obtained from the GRACE L2 product. The results show that, as compared to the isotropic filter, our non-isotropic filter may produce reasonably homogeneously smoothed results without excessive filtering in the latitudinal direction. However, one should bear in mind the non-isotropy of the filter when interpreting the results.

Our leakage reduction method is devised to reduce signal leakage of mass changes between land and ocean caused by Gaussian smoothing under the assumption that signal on land is far larger than that over ocean in coastal regions. The results are approximately equal to those of smoothing using a regional Gaussian filter with the same smoothing radius designed for land and ocean separately, that is, with no leakage present, while avoiding side lobes that may appear when directly smoothing GRACE mass change data regionally. However, the assumption is not necessarily true along every fraction of coasts. Hence the method should be used on a regional basis, as whether the basic assumption holds beyond a distance of the smoothing radius of the Gaussian filter would not significantly influence the results in the region being studied.

The efficiency of the leakage reduction method was shown using the trends of mass changes computed from both GLDAS data and GRACE L2 products. For GLDAS data, the results on land are identical to those obtained by filtering for land and ocean separately, that is, the leakage is completely removed. For GRACE data, we see clear signal reinforcement on land and drastic leakage removal over ocean along several coasts, particularly in Patagognia and the surrounding ocean.

Based on the formulation, our leakage reduction method could be modified to study mass changes in a selected region, either continental or oceanic, provided that the signal in the selected region can be considered predominant after removing the signal outside the selected region based on models. One such example is the case of Fenoglio-Marc et al. (2006) who studied mass changes over Mediterranean Sea by first removing the signal on surrounding land using data from a hydrological model.

However, if the signals in the region of interest (e.g. land in the paper) and in immediate surrounding regions (e.g. coastal oceans in the paper) are known to have comparable magnitudes, and there is no a priori information like a model to reduce the signal in surrounding regions to considerately lower magnitude, the approach of Wahr et al. (1998) may still be used to alleviate the leakage. We mention that sea level signals due to remote enough sources could be inferred from GRACE data (e.g. Baur et al. 2009).

GRACE L2 data are given in the form of truncated spherical harmonic series. Truncation introduces errors. In the scale factor method (Fenoglio-Marc et al. 2006; Velicogna & Wahr 2006; Baur et al. 2009) where the scale factor is obtained based on forward modelling, the truncation error is taken into account in meanwhile with leakage. The same is in the trial-and-error forward modelling method (Chen et al. 2006b). Our leakage reduction method for computing regionally smoothed mass changes does not reduce truncation error. Hence, if the truncation error is supposed to be significant, an extra step is required to reduce it. The most straightforward way to us is to use a scale factor similar to the methods mentioned above. In our case the scale factor is position dependent, and should be computed based on forward modelling using the mass change signal obtained after leakage reduction. The same as the methods mentioned above, the scale factor approach is empirical and case dependent, thus being not explored further here.

Finally, we mention that the discrete forms of the isotropic and our non-isotropic filters may be easily modified to freely define new filters with position-dependent smoothing radius or radii.

Acknowledgments

This research is partially supported by grants from NASA (NNG04GN19G, JPL-1356632 and JPL-1384376), and the Ohio State University Climate, Water and Carbon Program (http://cwc.osu.edu). XJD is partially supported by the Natural Science Foundation of China (60974124).

Appendix

Appendix A: The Cos-Fourier Expansion

Here we concisely summarize the formulae required in this work. The detail may be found in Canuto et al. (1987) and Guo & Shum (2009).

In our smoothing algorithm, we used the following cos-Fourier expansion to express f(θ, λ) 

(A1)
formula
We will use the values of f(θ, λ) on a latitude–longitude grid of constant latitude and longitude intervals to compute the cos-Fourier coefficients Amn and Bmn. The discrete form of (A1) can then be written as an inverse Fourier transform and an inverse Chebyshev transform of the spatial variable x= cos θ 
(A2)
formula
 
(A3)
formula
where θj and λk are the nodal values of colatitude and longitude.

The cos-Fourier coefficients Amn and Bmn are computed in two steps. The first step is a discrete Fourier transform:  

(A4)
formula
where  
(A5)
formula
and  
(A6)
formula
The second step is a discrete Chebyshev transform of the spatial variable x= cos θ. In this work, we choose to use the nodes over θ to be  
(A7)
formula

The transform is then  

(A8)
formula
where c0= 2, and cn= 1 if 0 < nN.

We see from (A5) and (A6) that the number of nodes in longitude is odd in the above formulation. However, it is mostly even in practical applications. For the number of nodes in longitude to be even, the terms including BMn should not be included in the cos-Fourier expansions. Hence, (A6) should be replaced by L= 2M− 1, and the factor 2/(L+ 1) in the formula of AMn in (A4) should be replaced by 1/(L+ 1).

References

Baur
O.
Kuhn
M.
Featherstone
W.E.
,
2009
.
GRACE-derived ice-mass variations over Greenland by accounting for leakage effects
,
J. geophys. Res.
 ,
114
,
B06407
, doi:10.1029/2008JB006239.
Canuto
C.
Hussaini
M.Y.
Qarteroni
A.
Zhang
T.A.
,
1987
.
Spectral Methods in Fluid Dynamics
 ,
Springer-Verlag
, New York.
Chambers
D.P.
,
2006
.
Evaluation of New GRACE Time-Variable Gravity Data over the Ocean
,
Geophys. Res. Lett.
 ,
33
,
L17603
, doi:10.1029/2006GL027296.
Chen
J.L.
Wilson
C.R.
Famiglietti
J.S.
Rodell
M.
,
2005
.
Spatial sensitivity of the gravity recovery and climate experiment (GRACE) time-variable gravity observations
,
J. geophys. Res.
 ,
110
,
B08408
, doi:10.1029/2004JB003536.
Chen
J.L.
Wilson
C.R.
Tapley
B.D.
,
2006
.
Satellite gravity measurements confirm accelerated melting of greenland ice sheet
,
Science
 ,
313
,
1958
1960
, doi:10.1126/science.1129007.
Chen
J.L.
Wilson
C.R.
Famiglietti
J.S.
Rodell
M.
,
2007
.
Attenuation effects on seasonal basin-scale water storage change from GRACE time-variable gravity
,
J. Geod.
 ,
81
(
4
),
237
245
, 10.1007/s00190-006-0104-2.
Chen
J.L.
Wilson
C.R.
Tapley
B.D.
Grand
S.
,
2007
.
GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earthquake
,
Geophys. Res. Lett.
 ,
34
,
L13302
, doi:10.1029/2007GL030356.
Davis
J.L.
Tamisie
M.E.
Elósegui
P.
Mitrovica
J.X.
Hill
E.M.
,
2008
.
A statistical filtering approach for Gravity Recovery and Climate Experiment (GRACE) gravity data
,
J. geophys. Res.
 ,
113
,
B04410
, doi:10.1029/2007JB005043.
Duan
X.J.
Guo
J.Y.
Shum
C.K.
Van Der Wal
W.
,
2009
.
On the postprocessing removal of correlated errors in GRACE temporal gravity field solutions
,
J. Geod.
 ,
83
,
1095
1106
, doi:10.1007/s00190-009-0327-0.
Fengler
M.J.
Freeden
W.
Kohlhaas
A.
Michel
V.
Peters
T.
,
2006
.
Wavelet modelling of regional and temporal variations of the Earths gravitational potential observed by GRACE
,
J. Geod.
 ,
81
,
5
15
, doi:10.1007/s00190-006-0040-1.
Fenoglio-Marc
L.
Kusche
J.
Becker
M.
,
2006
.
Mass variation in the Mediterranean Sea from GRACE and its validation by altimetry, steric and hydrology fields
,
Geophys. Res. Lett.
 ,
33
,
L19606
, doi:10.1029/2006GL026851.
Guo
J.Y.
Shum
C.K.
,
2009
.
Application of the cosine-Fourier series expansion in the transformation of data between latitude–longitude grids
,
Comput. Geosci.
 ,
35
,
1439
1444
, doi:10.1016/j.cageo.2008.09.010.
Guo
J.Y.
Li
Y.B.
Huang
Y.
Deng
H.T.
Xu
S.Q.
Ning
J.S.
,
2004
.
Green's function of the deformation of the Earth as a result of atmospheric loading
,
Geophys. J. Int.
 ,
159
,
53
68
, doi: 10.1111/j.1365-246X.2004.02410.x.
Han
S.C.
Shum
C.K.
Jekeli
C.
Kuo
C.Y.
Wilson
C.R.
Seo
K.W.
,
2005
.
Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement
,
Geophys. J. Int.
 ,
163
,
18
25
, doi:10.1111/j.1365-246X.2005.02756.x.
Jekeli
C.
,
1981
.
Alternative methods to smooth the Earth's gravity field
,
Rep. 327
 , Dep. of Geod. Sci. and Surv., Ohio State Univ., Columbus.
Klees
R.
Revtova
E.A.
Gunter
B.C.
Ditmar
P.
Oudman
E.
Winsemius
H.C.
Savenije
H.H.G.
,
2008
.
The design of an optimal filter for monthly GRACE gravity models
.
Geophys. J. Int.
 ,
175
,
417
432
, doi: 10.1111/j.1365-246X.2008.03922.x.
Kusche
J.
,
2007
.
Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models
,
J. Geod.
 ,
81
,
733
749
, doi:10.1007/s00190-007-0143-3.
Seo
K.W.
Wilson
C.R.
,
2005
.
Simulated estimation of hydrological loads from GRACE
,
J. Geod.
 ,
78
,
442
456
, doi:10.1007/s00190-004-0410-5.
Schmidt
M.
Fengler
M.
Mayer-Gürr
T.
Eicker
A.
Kusche
J.
Sánchez
L.
Han
S.C.
,
2006
.
Regional gravity modelling in terms of spherical base functions
,
J. Geod.
 ,
81
,
17
38
, doi:10.1007/s00190-006-0101-5.
Schrama
E.J.O
Wouters
B.
Lavallée
D.A.
,
2007
.
Signal and noise in Gravity Recovery and Climate Experiment (GRACE) observed surface mass variations
,
J. geophys. Res.
 ,
112
,
B08407
, doi:10.1029/2006JB004882.
Shum
C.K.
et al
,
2010
.
Inter-annual water storage changes in Asia from GRACE data
 , in Climate Change and Food Security in South Asia, eds Lal, R., Sivakumar, M., Faiz, S.M.A., Mustafizur Rahman, A.H.M. & Islam, K.R., Springer-Verlag, Dordrecht, in press.
Swenson
S.
Wahr
J.
,
2002
.
Methods for inferring regional surfacemass anomalies from gravity recovery and climate experiment (GRACE) measurements of time-variable gravity
,
J. geophys. Res.
 ,
107
(
B9
),
2193
, doi:10.1029/2001JB000576.
Swenson
S.
Wahr
J.
,
2006
.
Post-processing removal of correlated errors in GRACE data
,
Geophys. Res. Lett.
 ,
33
,
L08402
, doi:10.1029/2005GL025285.
Velicogna
I.
Wahr
J.
,
2006
.
Measurements of time-variable gravity show mass loss in Antarctica
,
Science
 ,
311
,
1754
1756
, doi:10.1126/science.1123785.
Wahr
J.
Molenaar
M.
Bryan
F.
,
1998
.
Time variablility of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE
,
J. geophys. Res.
 ,
103
(
B12
),
30205
30229
.
Wouters
B.
Schrama
E.J.O.
,
2007
.
Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spherical harmonics
,
Geophys. Res. Lett.
 ,
34
,
L23711
, doi:10.1029/2007GL032098.
Wu
X.
Blom
R.G.
Ivins
E.R.
Oyafuso
F.A.
Zhong
M.
,
2009
.
Improved inverse and probabilistic methods for geophysical applications of GRACE gravity data
,
Geophys. J. Int.
 ,
177
,
865
877
, doi:10.1111/j.1365-246X.2009.04141.x.