## 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 *C _{lm}* and

*S*is multiplied with a degree-only dependent factor

_{lm}*W*computed from the smoothing radius. The values of

_{l}*W*decrease as

_{l}*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

*W*with a degree-and-order dependent factor

_{l}*W*. The values of

_{lm}*W*depend on

_{lm}*l*in the same form as

*W*, but computed using larger smoothing radius for larger

_{l}*m*up to a certain limit

*m*

_{1}. 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*=

*m*

_{1}that the authors referred to as east–west smoothing radius. The smoothing radii for 0 <

*m*<

*m*

_{1}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:

where the standard deviation σ_{x}is related to the smoothing radius

*r*customarily used in literature (Jekeli 1981; Wahr

_{x}*et al.*1998), that is, the value of

*x*where

*G*(

*x*) =

*G*(0)/2, in the form of We have used the subscript

*x*in σ

_{x}and

*r*to specify that the variable is

_{x}*x*. In the following,

*x*will be replaced by the geocentric angle ψ defined below, and accordingly, σ

_{x}and

*r*will be replaced by σ

_{x}_{ψ}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 is defined as

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, *r _{d}*, that is related to

*r*

_{ψ}in the form of

Practical computation formulae using spherical harmonic expansion were developed by Jekeli (1981) who obtained the factor *W _{l}* 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 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 and by σ_{θ} and σ_{λ}, and *r _{x}* 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 to be

*W*(θ) is defined as For evident reason, it suffices to define the domain of integration to be for λ, and 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

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

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

*c*

_{0}= 2, and

*c*= 1 if 0 <

_{n}*n*≤

*N*. 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

where we used the relationIt 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:

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 |*B _{c}*|, the smoothing radii in terms of length at the Earth's surface along latitude,

*Rr*

_{θ}, and that along longitude,

*Rr*

_{λ}cos

*B*, are equal, and the non-isotropic filter become isotropic. Thus, we have cos

_{c}*B*=

_{c}*r*

_{θ}/

*r*

_{λ}. The classical isotropic Gaussian filter is then applied for |

*B*| > |

*B*|. For example, if we assign

_{c}*r*

_{λ}= 6°,

*r*

_{θ}= 3°, we have |

*B*| = 60°. In this way, the results of the two filters mesh seamlessly at |

_{c}*B*|.

_{c}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 |*B _{c}*| = 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.

## 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 and as the ocean area and land area, respectively, in the rectangle and . For mass changes *f*(θ, λ) at the Earth's surface, we modify (7) and (8) to define its smoothed value at the point to be

*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

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 *f*(θ_{j}, λ_{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: 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.

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

The values of*f*(θ

^{s}

_{α}, λ

^{s}

_{β}) are obtained by interpolation from the values of

*f*(θ

_{j}, λ

_{k}) (Guo & Shum 2009).

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

*m*and

*n*are integers, which we choose to be odd so that the grid cell centre of interest, , could be most likely at the centre of integral domain in (23) and (24). In fact, could be truly at the centre of integral domain only if the rectangle and 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 such that it includes only the continent or island where 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:

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

*f*(θ

_{s}_{j}, λ

_{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*_{3ψ} and *O*_{3ψ}, respectively. We then have

^{αβ}

_{jk}is computed using (4) with and (θ, λ) = (θ

^{s}

_{α}, λ

^{s}

_{β}), and

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.

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

where the GRACE error*e*is responsible for stripes. We use (where stands for ‘global’) to denote a filter without distinguishing land and ocean, and , the corresponding filter distinguishing land and ocean. The purpose of filtering is to reduce the influence of the error term

_{g}*e*. After smoothing without distinguishing land and ocean, the error term

_{g}*e*is very much reduced, and we have When the smoothing radius is large enough, we assume that the ‘smoothed error’

_{g}*e*

_{sg}is negligibly small, that is, the filter 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 . Our objective of this section is to find a way to recover from , which is assumed equivalent to as expressed in (35) where

*e*

_{sg}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, can be considered equivalent to , 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 hereafter, keeping in mind that it is actually obtained by computing . Here we state our objective more concretely: to recover from 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 on land, that will then be applied to recover 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, . This is quite straightforward. We use *s*_{glb}(θ_{j}, λ_{k}) and *s*_{dlo}(θ_{j}, λ_{k}) to express the values of and at the gridpoint (θ_{j}, λ_{k}). We can express the result as

*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 from , we suppose (θ_{j}, λ_{k}) is on land, and rewrite the formulae of the filter distinguishing land and ocean, that is, , in the form of

*e*

_{ol}, 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

*s*

_{glb}(θ

_{j}, λ

_{k}) by the smoothed result of GRACE data,

*f*

_{glb}(θ

_{j}, λ

_{k}), to compute

*s*

_{dlo}(θ

_{j}, λ

_{k}).

The error *e*_{ol} may be obtained by comparing (36) and (38) to (41)

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

We then smooth*s*using the isotropic Gaussian filter without distinguishing land and ocean, and use

^{L}*s*

^{L}_{glb}(θ

_{j}, λ

_{k}) to denote the result at the point (θ

_{j}, λ

_{k}). The values of

*s*

^{L}_{glb}(θ

_{j}, λ

_{k}) over ocean are then the leakage from land. The smoothed signal over ocean after removing leakage, , is then

*s*

_{glb}(θ

_{j}, λ

_{k}) −

*s*

^{L}_{glb}(θ

_{j}, λ

_{k}). Here we mention again that

*s*

_{glb}(θ

_{j}, λ

_{k}) is to be replaced by

*f*

_{glb}(θ

_{j}, λ

_{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 , 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.

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*| < =*B _{c}*, and the classical isotropic Gaussian filter when |

*B*| >

*B*. The smoothing radii of the filters are suitably chosen so that the results of the two filters mesh seamlessly at |

_{c}*B*| =

*B*.

_{c}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*(θ, λ)

*f*(θ, λ) on a latitude–longitude grid of constant latitude and longitude intervals to compute the cos-Fourier coefficients

*A*and

^{m}_{n}*B*. The discrete form of (A1) can then be written as an inverse Fourier transform and an inverse Chebyshev transform of the spatial variable

^{m}_{n}*x*= cos θ where θ

_{j}and λ

_{k}are the nodal values of colatitude and longitude.

The cos-Fourier coefficients *A ^{m}_{n}* and

*B*are computed in two steps. The first step is a discrete Fourier transform:

^{m}_{n}*x*= cos θ. In this work, we choose to use the nodes over θ to be where

*c*

_{0}= 2, and

*c*= 1 if 0 <

_{n}*n*≤

*N*.

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 *B ^{M}_{n}* should not be included in the cos-Fourier expansions. Hence, (A6) should be replaced by

*L*= 2

*M*− 1, and the factor 2/(

*L*+ 1) in the formula of

*A*in (A4) should be replaced by 1/(

^{M}_{n}*L*+ 1).

## References

*Spectral Methods in Fluid Dynamics*

*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.