1. Introduction
The Gravity Recovery and Climate Experiment (GRACE) satellites, launched on March 17, 2002, provided global gravity solutions from April 2002 to June 2017. Temporal gravity variations from GRACE have been used to measure water and ice mass load redistribution after removal of effects due to atmospheric mass, ocean dynamics [
1] and glacial isostatic adjustment in Earth’s mantle [
2]. GRACE has provided the leading data type to understand contemporary ice mass loss in polar regions (e.g., [
3]) and subsequent sea level rise (e.g., [
4]), variations of total water storage (e.g., [
5]) and hydrological components, such as river discharge and evapotranspiration (e.g., [
6,
7]).
GRACE gravity solutions are provided as changes in spherical harmonic (SH) coefficients and subsequently as changes in mass concentrations (mascon) associated with defined geographical regions. GRACE SH coefficients are contaminated by correlated errors associated with satellite orbit resonance [
8] and other measurement noise, which are suppressed using various filtering steps. Correlated errors produce a north-south striping pattern, which is suppressed by a decorrelation filter [
9]. Other noise, especially at high SH degree and order, is normally suppressed by Gaussian smoothing [
10] or special optimum filters [
11]. These filters are effective but also cause reduced spatial resolution. The resulting smoothing is particularly problematic in coastal areas because larger land signals contaminate nearby ocean signals. This spatial leakage from land to oceans causes signal amplitude on land to be reduced, while at the same time masking residual ocean dynamic and self-attraction and loading signals [
12].
Forward Modeling (FM) [
13] is an iterative algorithm that estimates leakage-corrected mass loads by empirically updating terrestrial mass loads until the smoothed updated mass load converges to the smoothed GRACE signal. The FM method can correct leakage into the oceans to estimate average terrestrial mass changes but the spatial distribution over land is not necessarily correct unless surface mass load geography is known in detail. This is often the case for ice mass loss at glacial outlets [
14] but for terrestrial water storage changes, change locations are not well-known in advance.
In this study we develop a linear least squares method to correct for spatial leakage that, like FM, works directly with smoothed SH solutions. Spatial regularization constraints are used to preserve general spatial patterns of mass change found in the SH solutions. The method is tested using synthetic GRACE data and applied to real GRACE data to estimate terrestrial mass loads in Australia. Estimated surface mass loads show effective suppression of leakage into the oceans, as with FM, while retaining spatial distribution of surface mass load indicated in the original solutions.
2. Method
A surface mass load estimate is
, a vector containing mass changes for locations on Earth’s surface, normally associated with defined geographical regions, for example 1 × 1 degree pixels. Observed values from GRACE are
. They are linearly related by some matrix G
with
a vector of prediction misfits. For
, previous studies have used gravitational potential differences [
15] or gravitational acceleration [
16]. Potential differences were estimated from range-rate perturbations of GRACE satellites, while gravitational acceleration is obtained from GRACE SH coefficients. In this study we use smoothed surface mass loads in all pixels as
. In this case,
is approximately the convolution operator associated with Gaussian smoothing. Using the continent of Australia as an example, we divide it into 1 × 1 degree pixels, including regions that extend 400 km into the oceans.
Figure 1a shows grid points for both
and
and black and gray pixels represent land and ocean grid points, respectively. We build the convolution operator
G based on the Green function representing 400 km Gaussian smoothing for a unit mass (
Figure 1b). Thus, elements of
G,
, are values shown in
Figure 1b determined by the angular difference between
ith and
jth grid points.
can be estimated by minimizing the squared misfit:
but generally, this estimation problem is ill-posed because it is a downward continuation of the smoother gravity field at satellite altitude to Earth’s surface where the sources are located. To stabilize the solution, a regularization parameter
[
17] can be used:
where
is an identity matrix. The regularization parameter can be determined empirically or objectively with the L-curve method (e.g., [
18]).
In this study we include a penalty that seeks to minimize spatial gradient (steepness) in the cost function [
19]. This regularizes the problem and suppresses noise amplified as a result of the deconvolution. The steepness can be represented using the gradient of two surface mass loads. For example, for steepness between the first surface mass load,
, and others (
,
matrix can be constructed:
in which
is the distance between the first and
nth mass load. Therefore, steepness between the first mass load and all others is
. Similarly, a steepness matrix (
) between the second and all other mass loads can be constructed and similarly for loads in all pixels. We define matrix
considering all:
Including the steepness, the least squares solution is:
where
and
is a regularization parameter (different from
in Equation (3)). Larger value of
will suppress more noise but the estimate will be smoother. We empirically determine
in the next section.
Estimated
over oceans (
via Equation (6) includes three sources: ocean mass signal determined by self-attraction and loading (SAL) effects [
12], residual ocean dynamic effects after applying atmospheric and ocean de-aliasing (AOD) fields from numerical models [
20] and spurious leakage from terrestrial water mass load (
). Terrestrial water mass loads cause leakage into oceans and at the same time, leakage into land by oceanic mass loads is also present. In general, the leakage into the oceans is much larger than leakage onto land. Considering both leakage effects is complicated. In this study, we only attempt to suppress spatial leakage from land to oceans by assuming that
is zero. This constraint effectively ignores local ocean signals. It can be realized via an additional linear relationship:
where
is ones for ocean elements in
and zeros elsewhere. This linear relationship can adjust
by constraining
to be zeros. The result is estimation of leakage corrected terrestrial mass loads.
Equation (6) can be modified by including Equation (7) in the cost function and the two equations must be solved simultaneously by combining them into the matrix form (see Section 3.10. in [
21]):
where
is a null parameter. We finally obtain the equation for
, which refers to an inversion solution, representing a constrained deconvolution of the GRACE data. We have yet to determine a useful value for the regularization parameter:
3. Synthetic Experiment in Australia
We apply Equation (9) to estimate surface mass loads over Australia using synthetic GRACE data. For surface mass load, surface water storage from Global Land Data Assimilation System (GLDAS) is used [
22]. Surface water storage includes components of soil moisture, snow, and canopy water. Ocean dynamic effects are nominally corrected in GRACE data processing but residual ocean bottom pressure (OBP) signals, which may cause leakage contamination from oceans to land, may remain. We consider such effects using average ocean bottom pressure from CSR [
23], GSFC [
24] and JPL [
25] mascons. Variations of total surface water storage and residual ocean dynamic effects are converted to SH coefficients up to degree and order 60. We also include GRACE noise, taking it to be the difference between GRACE data and smoothed GRACE data [
6]. This difference field includes GRACE noise, which likely dominates spatially high frequency components, but it may also include some GRACE signal. The residual signal can be additionally separated from noise using Empirical Orthogonal Functions (EOF) [
6]. The synthetic GRACE data is reduced by conventional GRACE data reduction methods, applying a de-correlation filter [
9] and 400 km Gaussian smoothing.
Figure 2a shows GLDAS surface mass loads over Australia and residual OBP for September 2005 and
Figure 2b shows synthetic GRACE data, the sum of surface mass load shown in
Figure 2a and estimated GRACE noise after applying decorrelation filter. GRACE noise is predominant in the synthetic GRACE data.
Figure 2c shows ‘reduced’ surface mass load from synthetic GRACE data after applying a de-correlation filter and 400 km Gaussian smoothing. This is surface water mass load as in a conventional GRACE estimate. The mass load field is global but we only consider Australia and the surrounding 400 km of oceans. Spatial leakage of terrestrial signals into the oceans is evident, leading to reduced signal strength over land along the coast.
Figure 2d shows
, in Equation (9) estimated from
Figure 2c as
for a particular
, which was chosen by varying λ, and seeking the best agreement based on the regression coefficients between synthetic surface mass load (modeled by GLDAS) and estimated surface mass loads,
, and the root mean squared (RMS) difference between the two.
Table 1 shows mean RMS and regression coefficients from Equation (9) with varying λ and different amounts of Gaussian smoothing during the entire period (January 2003–August 2016) of synthetic data. We determine
to estimate
because the regression coefficient is larger than any value from Gaussian smoothing while the RMS difference is close to the case of 300 km Gaussian smoothing, which provides the smallest RMS among Gaussian smoothing. Smaller λ provides a larger regression coefficient but simultaneously larger RMS difference. Estimated loads from Equation (9) with
(
Figure 2d) are closer to GLDAS fields relative to those in
Figure 2c. For example, large leakage into the oceans is found in South Western and North Eastern Australia in
Figure 2c and signal leaking to the oceans is moved back to the continent in
Figure 2d. Estimated loads in Tasmania are evidently different from GLDAS loads due to its small spatial scale and the presence of large GRACE noise (shown in
Figure 2b).
Figure 2e shows
from Equation (6) without considering the constraint that
are zeros, which is realized in Equation (7). The estimated surface mass loads shown in
Figure 2e are similar to
Figure 2c, the case of 400 km Gaussian smoothing. The comparison between
Figure 2d,e clearly exhibits effect of the leakage correction in Equation (9), incorporating the constraint.
Figure 3 shows scatter plots between GLDAS and estimated surface mass loads from 400 km Gaussian smoothing (blue) and inversion (red) shown in
Figure 2c,d, respectively. The linear regression coefficient from inversion (0.69, r
2 = 0.68) is much higher than that from Gaussian smoothing (0.43, r
2 = 0.65), with higher significance determined by the coefficient of determination, r
2. The RMS value from inversion (28.56 mm) is also much lower than that from 400 km Gaussian smoothing (32.78 mm).
Because Gaussian smoothing with a 300 km length scale is superior to other choices shown in
Table 1, we compare terrestrial water storage (TWS) estimates from Equation (9) for the 300 km case. Gray, black and red lines in
Figure 4a show mean terrestrial water storage (TWS) time series over Australia associated with GLDAS, 300 km Gaussian smoothing and Equation (9) respectively. Both red (inversion) and black lines (300 km Gaussian smoothing) are close to GLDAS (gray) and similar to one another. As shown in
Figure 4b, there are similar differences between estimates and GLDAS loads. This occurs because the sign of leakage differs across the continent. The sign is negative in North Australia and positive in the south as shown in
Figure 2c. They tend to cancel in a continental scale average in this synthetic experiment. Cancellation would not be as important for real GRACE data as discussed in the next section.
The synthetic data experiment shows that constrained deconvolution, via Equation (9), with proper choice of λ, is superior to ordinary SH estimates subjected to Gaussian smoothing. TWS changes are closer to GLDAS values by suppressing leakage into the oceans, but they preserve TWS spatial patterns evident in the original SH solutions.
4. TWS in Australia
We now use observed GRACE data (CSR RL06 from January 2003 to August 2016) with Equation (9) to estimate TWS change in Australia. Initially we used the λ value determined from the synthetic data experiment but subsequently modified it as discussed shortly. Glacial isostatic adjustment is corrected by a model prediction based on ICE-6G_C ice thickness history and VM5a radial viscosity profile [
26], although its effect is very small in Australia. North-south stripes are removed by a decorrelation filter [
9].
Figure 5a shows surface mass loads for September 2005, the same month as the synthetic data experiment (
Figure 2), including 300 km Gaussian smoothing. Spatial leakage into the oceans is evident. Pixels are 1 × 1 degrees over Australia out to 400 km from the coast, as in the synthetic experiment.
Comparing GLDAS-derived (
Figure 2c) and GRACE fields (
Figure 5a) shows that GLDAS has likely underestimated the magnitude of Australian TWS change. A larger TWS signal relative to GLDAS is evident also in
Figure 6a, where continent-wide variations from GRACE are larger than GLDAS values.
The value for
is chosen to suppress noise but it also suppresses signal, so must depend upon the signal to noise ratio. From the synthetic data experiment
but if the signal level is higher, relative to GLDAS used in the synthetic experiment, then a smaller value may be appropriate. The RMS value from the black line in
Figure 6a is about 1.3 times larger than that of the black line in
Figure 4a. The discrepancy reflects lower signal strength in GLDAS relative to GRACE. Therefore, we returned to the synthetic experiments to adjust GLDAS variations by a factor of 1.3 to compensate for the underestimate of GLDAS. From this, we selected an optimum
.
Figure 5b shows estimated TWS change from Equation (9). Large TWS signals appear along coastal areas relative to conventional Gaussian smoothed estimates (
Figure 5a). The red line in
Figure 6a shows resulting continent-wide TWS variations. Variability in black and red lines in
Figure 6a is very close to each other (correlation coefficient, 0.99) while peak-to-peak changes in the red line are much larger than the black line. An additional comparison is shown in
Figure 6a with CSR mascon solutions [
23] (blue line). In general, the CSR mascon solutions show smaller peak-to-peak variations relative to the Equation (9) solution and they are rather close to the case of Gaussian smoothing. We also compare GSFC [
24] and JPL [
25] mascon solutions with inversion (
Figure 6b). Smoothing effects in JPL mascon solutions were corrected by applying gain factors [
27]. Similar to CSR mascon solutions, GSFC and JPL mascon solutions also show smaller TWS variations than the Equation (9) solution.
Figure 5c–e show TWS from CSR, JPL and GSFC mascons respectively, for the same month as in
Figure 5a,b. TWS spatial patterns are different from one another while the CSR mascon solution (
Figure 5c) is close to the inversion estimate (
Figure 5b).
We further compare TWS from Gaussian smoothing, mascons and inversion at basin scales to examine their uncertainty. First, we choose the Timor Sea Basin, located on the north coast of Australia, showing large TWS variations. It is a good choice for studying leakage effect and corrections among the three methods.
Figure 7a shows TWS variations from the three methods. Both CSR mascon (blue) and inversion (red) provide similar TWS variations while their peak-to-peak variations are much larger than those from Gaussian smoothing (black). This indicates that the mascon and Equation (9) solutions suppress leakage effects in basins near the oceans. Similar results are also found comparing JPL mascon (blue) and inversion (red) in
Figure 7b. The GSFC mascon solution (black), however, evidently underestimates TWS variations relative to inversion and other mascons.
Figure 7c shows a synthetic experiment for the Timor Sea Basin, where inversion (Equation (9)) TWS (red) is much closer to GLDAS (gray) than the Gaussian estimate (black). Some discrepancies between red and gray curves are evident in
Figure 7c, indicating leakage may not be completely corrected using Equation (9). Despite this likelihood in the case of real GRACE data shown in
Figure 7a, Equation (9) estimates removed much of the leakage error associated with Gaussian smoothing. A similar limitation appears in mascon solutions, which show very similar TWS variations in
Figure 7a,b.
We apply the same analysis to the Murray Darling Basin, located in the continental interior. It differs from the Timor Sea basin in distance from the ocean, with small TWS variations over a basin of much larger size. The Murray Darling and Timor Sea Basin areas are 1.15 × 10
6 km
2 and 0.62 × 10
6 km
2 respectively.
Figure 8a shows TWS variations over the basin from Gaussian smoothing (black), CSR mascon (blue) and Equation (9) inversion (red). Inversion and Gaussian smoothing provide very similar TWS variations because the effect of smoothing should be small in a large interior basin. The same conclusion can be found in the synthetic experiment (
Figure 8c): both Gaussian smoothing and Equation (9) show similar variations as GLDAS. However, mascons (blue line in
Figure 8a and black and blue lines in
Figure 8b) are distinctly different. The cause of this discrepancy is unknown but it suggests that CSR, GSFC and JPL mascon solutions may not perform well in certain situations. Similar results are found in other basins in Australia, such as the Lake Eyre Basin.
5. Conclusions
We have developed a constrained deconvolution method to estimate surface mass loads from GRACE data. Noise is suppressed by minimizing steepness of the estimated surface mass loads in the cost function. Signal leakage from land to oceans is suppressed by constraining ocean signals to be approximately zero. Testing this method using synthetic GRACE data shows results that are superior to Gaussian smoothing.
Using this method, we estimate surface mass loads over Australia from GRACE data. As in the synthetic case, estimates show higher signal strength near the coast relative to Gaussian smoothing but retain high correlation with Gaussian smoothing estimates. TWS variations are also compared with those from CSR, GSFC and JPL mascon solutions. In the examples studied here, mascon solutions suppress leakage into the oceans but show smaller peak-to-peak variation and may be problematic in interior basins with small signals.
The constrained deconvolution method should be useful in geographical regions adjacent to oceans or lakes to suppress leakage error. Calibration of λ for particular problems is necessary. We show that it likely depends both on signal strength and probably also, to some extent, on pixel size. However, the process should be straightforward, using methods similar to those in the synthetic data experiment. Additional estimates incorporating other data types, such as GPS loading observations, would be particularly useful to enhance spatial resolution while reconciling leakage problems.