Observed changes in the Earth’s dynamic oblateness from GRACE data and geophysical models

A new methodology is proposed to estimate changes in the Earth’s dynamic oblateness (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {J_{2}}$$\end{document}ΔJ2 or equivalently, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\sqrt{5}\Delta {C_{20}}$$\end{document}-5ΔC20) on a monthly basis. The algorithm uses monthly Gravity Recovery and Climate Experiment (GRACE) gravity solutions, an ocean bottom pressure model and a glacial isostatic adjustment (GIA) model. The resulting time series agree remarkably well with a solution based on satellite laser ranging (SLR) data. Seasonal variations of the obtained time series show little sensitivity to the choice of GRACE solutions. Reducing signal leakage in coastal areas when dealing with GRACE data and accounting for self-attraction and loading effects when dealing with water redistribution in the ocean is crucial in achieving close agreement with the SLR-based solution in terms of de-trended solutions. The obtained trend estimates, on the other hand, may be less accurate due to their dependence on the GIA models, which still carry large uncertainties.


Introduction
Monthly Earth gravity field models based on data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission   (Liu et al. 2010). In spite of continuous improvements in data processing techniques, very low-degree spherical harmonic coefficients still cannot be determined with high accuracy. This is largely due to the mission design (low orbits, limited separation of the satellites, etc.) (Chen et al. 2005). In particular, this concerns variations of the C 20 coefficient ( C 20 , denoted as C 20 hereafter for simplicity), which describes changes of the Earth's dynamic oblateness J 2 (J 2 = − √ 5C 20 , where the factor √ 5 implicitly means that the C 20 is normalised). Estimations of this coefficient are corrupted by 161-day-period ocean tide aliases (Cheng et al. 2013). Therefore, the C 20 coefficient in GRACE gravity field models is recommended to be replaced with estimates from other techniques such as satellite laser ranging (SLR), which is likely to provide the most accurate C 20 information so far (Cheng and Tapley 2004).
An alternative source of information about variations of low-degree coefficients is surface mass loading inferred from the GPS-sensed solid Earth deformation, an approach known as the inversion method Gross et al. 2004;Wu et al. 2012). Swenson et al. (2008) developed a new method to determine the degree-1 coefficients by combining GRACE information with ocean bottom pressure (OBP) data, so that the usage of GPS data is not needed.
Here, we extend the methodology by Swenson et al. (2008) further to estimate the monthly C 20 coefficients from other GRACE gravity field model coefficients supported by the C 20 coefficients from an OBP model and a glacial isostatic adjustment (GIA) model. We validate our solutions against SLR-derived estimates. This study is motivated by the following considerations: (1) dense and evenly distributed measurements are used as the input. (2) The proposed procedure has better prospects regarding an increasing accuracy of future satellite gravity mission and related geophysical models. In addition, one will be able to use the proposed procedure for a mutual validation of the estimates based on GRACE data and on other techniques.

Methodology
Following equation (11) in Swenson et al. (2008), one can derive a similar equation for the determination of the C 20 coefficient: where C ocean 20 represents the oceanic component of C 20 . Integrals are defined over the entire globe, dΩ = sin θ dθ dφ is an element of solid angle. The summations exclude the estimated term C 20 . Indices l and m stand for spherical harmonic degree and order, respectively.P lm are normalised associated Legendre functions. θ is colatitude in spherical coordinates, φ is longitude, ϑ(θ, φ) denotes the ocean function, which equals 1 over ocean and 0 over land. C 20 , C lm and S lm denote the "mass coefficients" describing the surface mass change and are related to the dimensionless Stokes coefficients C 20 , C lm and S lm by in which a is the semi-major axis of the reference ellipsoid, ρ earth is the Earth's average density and k l denotes the degree-l load Love number (Wahr et al. 1998). Following Swenson et al. (2008), one can easily extend equation (1) to the case when four coefficients-C 10 , C 11 , S 11 , and C 20 -have to be simultaneously estimated, for which purpose a system of linear equations has to be solved: where the following notations have been used: and (similar for the other elements of vector G), in which the summations exclude the terms that are estimated. To solve the system of linear equations and obtain degree-1 and C 20 dimensionless Stokes coefficients, one needs (1) the oceanic component of degree-1 and C 20 , (2) higher-order Stokes coefficients and (3) GIA model coefficients. The input and output shown in the equation are mass coefficients, but they are directly related to the Stokes coefficients mentioned here through Eq. (2). The Stokes coefficients used in this study come directly from the GRACE level-2 data (also known as GSM), for which the oceanic and atmospheric mass variations are subtracted using the atmosphere and ocean de-aliasing level-1B (AOD1B) products (Flechtner et al. 2014). Monthly averages of the AOD1B product are available in Stokes coefficients stored in GAC and GAD files. GAC includes global oceanic and atmospheric effects, while GAD has the atmospheric contribution over land set to zero. To make the input coefficients compatible, the same oceanic and atmospheric effects need to be removed also from the oceanic coefficients, e.g., C ocean 20 .
Since the oceanic coefficients lack the contribution from atmosphere over continents, it is the GAD (rather than GAC) which should be subtracted. With this procedure, the output will also be GSM-like coefficients. If the full C 20 coefficients are needed, the contribution of GAC can be restored afterwards.
An alternative procedure requires that the AOD1B product is first added back to GSM coefficients and then full degree-1 and C 20 coefficients are estimated directly. Although the latter procedure is stated to be equivalent to the first one in Swenson et al. (2008), it is not favoured in this study for the reason outlined in Sect. 3.2.

Oceanic C 20
As has been mentioned above, the GAD contribution, denoted as C GAD 20 , needs to be removed from C ocean 20 coefficients. The GAD coefficients represent the OBP model that describes the pressure on the sea floor from both air and water column above. The water columns are output from the ocean model from circulation and tides (OMCT) (Thomas 2002). This ocean model applies the Boussinesq approximation and thus essentially conserves the ocean volume. A thin uniform layer of water is then added or removed to conserve the total ocean mass. can also be provided by GRACE data (Chambers and Schroeter 2011). In this study, we integrate the GSM coefficients over the continental areas to infer the total water mass variations (which are opposite to the water mass variations in the oceans, assuming mass conservation in the Earth system). Once the monthly water mass variation is known, the value of C exchange 20 can be obtained by assuming a certain spatial distribution of the exchanged water over the oceans. We implement two different approaches: (1) water redistributes as a uniform layer [eustatic approach, as in Swenson et al. (2008)]; (2) water redistributes accounting for Self-attraction and loading effects (SAL approach). SAL effects [or fingerprints, Mitrovica et al. (2001)] are computed by solving the sea-level equation , including the feedback from Earth rotation (Milne and Mitrovica 1998). It is worth noting that using GRACE to constrain total mass change over the continents requires the availability of a complete GRACE solution, which includes the coefficients being estimated through Eq. (3). Therefore, C exchange 20 needs to be determined through an iterative approach (starting from a GRACE solution where the four estimated coefficients are null, later updated with preliminary estimates of the same coefficients). Convergence is very quick, with the difference between subsequent solutions being smaller than 0.1 % in 3 or 4 iterations.
The degree-1 coefficients are estimated similarly, simultaneously with C 20 .

GRACE gravity field models
In this paper, we present results based on CSR RL05, GFZ RL05a and JPL RL05 time series in the period from Jan-uary 2003 to May 2013, all complete to or truncated at degree 60.
All the GRACE-based monthly gravity fields contain spatially correlated noise that reveals itself in the form of meridionally oriented stripes in the spatial domain. To solve the sea level equation and account for self-attraction and loading effects, we need to know the spatial distribution of the land load as accurately as possible. For this purpose, we use publicly available solutions that have been post-processed by means of the DDK4 filter ) (http://icgem. gfz-potsdam.de/ICGEM/). The DDK4 filter is a decorrelation filter making use of error covariance matrices, and an a priori signal covariance matrix in the spherical harmonic domain. In this way, the filtering ensures that a higher noise or/and lower signal level means harder damping and vice versa. Ultimately, the effect of this filter is somewhat similar to that of a combination of empirical destriping algorithm (Swenson and Wahr 2006) and Gaussian filter (Wahr et al. 1998).
When using Eq. (3), we need to deal with the limited spatial resolution of the GRACE gravity field models, which causes signals to spread over (or leak into) wider areas. The signal leakage is further increased by applying a filter, such as DDK4. As a result, the available observations cannot distinguish whether mass variations occurring in coastal areas are originating from the land or from the ocean. An attempt to define an ocean function without taking this fact into account may lead to a miscalculation of the total mass exchange between land and oceans as well as of the G vector. We correct for signal leakage by introducing a buffer zone around all land areas, similarly to what is done by Swenson et al. (2008) when computing the total ocean mass change. Differently from that study, we also consider the buffer zone to be part of the land areas when we define the ocean function ϑ(θ, φ), which means that we include the buffer zone in the definition of the G vector. We will show that such a buffer is crucial to obtain solutions close to SLR estimates. The use of a buffer introduces the risk that mass redistribution due to ocean dynamical processes in coastal areas is erroneously attributed to land processes. However, the problem is largely reduced using GSM coefficients in Eq. (3), under the assumption that the AOD products capture most of the ocean signal.

GIA models
The method discussed above and Eq. (2) imply that gravity field variations are solely due to a redistribution of mass at the Earth's surface. Solid Earth contributions such as those of tectonics and GIA should, therefore, be removed. Here, only GIA is accounted for as proposed by Swenson et al. (2008). The removed GIA signal is restored at the final data processing stage. Since GIA is characterised by a linear trend, the choice of a specific GIA model has no impact on seasonal and other short-term signals.
Considered that available GIA models are highly uncertain, we only show the resulting C20 trends for a few GIA realisations, based on different Earth rheologies and on two Antarctic ice histories. A full-scale sensitivity study is beyond the scope of this paper.
Four GIA models have been used in this study. All models are based on the ICE-5G ice history (Peltier 2004). Model-A, -B and -C are based on a simplified version of viscosity model VM2 (Peltier 2004), while Model-D assumes a lower mantle with a higher viscosity (10 22 Pa s) than VM2 (Mitrovica and Forte 1997). Model-A is taken from A et al. (2012)

Results
The following factors can affect the estimation of C 20 coefficients: (1) the choice of the input models (GRACE solutions, OBP and GIA models) and (2) implementation details (buffer zone width, the filter applied to GRACE solutions and whether or not accounting for self-attraction and loading effects). By trying different combinations of data processing parameters, we produced many variants of C 20 time series. Each of them was compared with the state-of-the-art C 20 time series based on SLR data from five geodetic satellites (LAGEOS-1 and 2, Starlette, Stella and Ajisai) (Cheng et al. 2013). Since all the results discussed are presented in the form of GSM-like coefficients, the AOD1B product (GAC coefficients) has also been removed from the reference SLR time series. We estimate bias, linear trend, acceleration, as well as annual and semi-annual periodic terms for each time series and make a comparison with corresponding parameters derived from the SLR-based time series.
We first compare de-trended (linear-trend removed) time series both visually and in terms of variance, where the percentage of the SLR variance explained is defined as R 2 = 1 − SLR − MODEL / SLR , where MODEL represents our estimation in this study and denotes the variance operator. We also compare annual amplitudes and phases against those of the SLR solution. Comparison of de-trended time series will lead to results invariant to the GIA model used. Later, we use one selected solution to compare the linear trend estimates resulting from different GIA models.

Seasonal variations
In Fig. 1, we show a few time series meant to illustrate the sensitivity of our GRACE-based solutions to implementation details and input models. The reference SLR solution is represented by a black solid line and by a grey band, indicating mean value and one standard deviation, respectively. In Table 1, we show statistics for the same models, as well as for a few additional experiments (different buffer widths, use of the DDK4 filter).
In Fig. 1a, we show the role of implementation details, namely of the use of a buffer zone and of the computation of SAL effects, based on GRACE CSR RL05 solutions. Not using any buffer and ignoring SAL effects (green line) largely underestimates the amplitude of the seasonal cycle. Nonetheless, most features of the SLR time series are already recognisable, such as the relative size of maxima and minima, as well as their phase. This solution explains about 59 % of the SLR variance, where the annual cycle is rather close in phase, but clearly smaller in amplitude (65 % of SLR). The addition of a 200 km buffer zone (blue line) largely improves the overall (explained variance) fit as well as the size of the peak amplitudes. The amplitude of the annual signal becomes statistically equivalent (within 2σ ) to the SLR solution. However, the improvement on the overall fit is moderate, where the new solution explains about 68 % of the SLR variance. Further increase in the width of the buffer zone to 250 and 300 km will begin to lower the explained variance slightly. When using a 300 km buffer width, the annual amplitude estimated becomes smaller. The analysis of the buffer zone width and more advanced way of handling signal leakage will be discussed in a separate paper.
Finally, accounting for SAL effects (red line) further improves the explained variance and at the same time significantly affects the amplitude of the estimated annual signal. The solution closest to SLR when including SAL effects (the explained variance is 71 %) makes use of a smaller buffer (150 km) than in the eustatic case. Note that the effects of feedback from the Earth rotation are accounted for during the computation of SAL. These effects on the estimated C 20 coefficients are negligible (not shown).
It is worth mentioning that the elimination of the buffer zone from the ocean functions prevents the accounting for SAL effects in the coastal regions. We have verified, however, that this has a little impact on the solution. We have considered the following two scenarios: (1) solving the sea level equation for the whole ocean; (2) solving the sea level equation for a slightly smaller ocean by reducing the ocean function 150 km along all boundaries while keeping the continental load unchanged (i.e., ignoring the mass variation inside the 150-km-wide buffer zone seen by GRACE). The resulting amplitude of the annual signal in the second sce-  nario increases, compared to the first one, by only about 2 %, which is less than the uncertainty. In Fig. 1b, we fix the implementation parameters and show the effect of using different GRACE solutions. The GFZ solution provides the best overall fit (71.6 % of the SLR variance explained and same amplitude of the annual signal). Nonetheless, all three time series-GFZ, CSR and JPLare very close to each other and the amplitude of the annual signal is statistically equivalent (within 1σ ).
The phase estimates are not significantly affected by any of the above-mentioned factors. The differences of phase estimates compared to those based on SLR data are all within ten days. Table 1 also lists linear trend estimates when using GIA model Model-C. Note that those trends are still based on the GSM-like solutions, but we have verified that long-term trends in atmospheric pressure over land and in OBP are negligible. The table shows that both buffer and SAL effects have a large impact on the trend due to the present-day mass transport (PDMT). The estimated trend is zero without buffer and SAL, but eventually becomes 40 % larger than the SLR trend for the model that provides the best fit of the seasonal signal. The largest effect originates from the buffer, but also SAL effects are sizeable (causing a further increase of up to 14 % when the buffer width is 200 km). The trend is based on GIA realisation Model-C, where the GIA contribution to the trend has been restored. The solution SLR FULL, where the AOD1B fields have not been removed, is provided as a reference. The highlighted solution (in bold) is recommended and is available online at http://www.citg.tudelft.nl/c20 In Table 2, we list the effect of using different GIA models for the results based on DDK4-filtered CSR solutions in combination with a 150-km buffer and taking SAL effects into account (i.e., CSRDDK4 + BUF150 + SAL). Similar conclusions hold for other setups. To allow an easier comparison with previous studies, we show the obtained trends in terms ofJ 2 .

Trend estimates and GIA
The use of GIA models allows us to separate the contribution of GIA from that of PDMT. The GIA contribution is uniquely defined for each model, while the PDMT value depends on the full GIA spectrum and is, therefore, affected by implementation details.
The smallest (in absolute value)J 2 of GIA comes from the model by A et al. (2012) (Model-A) which at the same time produces a relatively large estimate for the contribu-tion of PDMT, leading to a largerJ 2 value than the model Model-B based on an incompressible earth. Substituting the Antarctic contribution of ICE-5G with results based on IJ05 (Model-C) has no impact onJ 2 of GIA alone, likely due trade offs between the different ice history and the different viscosity structure used for the Antarctic model. However, the use of IJ05 does reduce the mass loss estimate from Antarctica, leading to a smaller PDMT contribution and to the smallest totalJ 2 . A higher viscosity in the lower mantle (Model-D) leads to larger contributions from both GIA and PDMT, which compensate each other and result in the second smallest totalJ 2 .
None of the GIA models tested here provides a very good fit to theJ 2 value determined from SLR. However, our results show positive sign ofJ 2 , confirming the findings from earlier studies on the inversion ofJ 2 observed since 1998 (Cox and Chao 2002), which has been attributed to an increased contribution from PDMT (Dickey et al. 2002;Cheng and Tapley 2004;Nerem and Wahr 2011;Cheng et al. 2013).

Eustatic sea-level variability and geocentre motion
Finally, it is worth having a brief look at two byproducts of our study: the solutions for eustatic sea-level variability (see Table 3) and for geocentre motion simultaneously obtained with C 20 (see Table 4).
The eustatic sea-level variability estimated using the approach described in Sect. 3.1 has been compared with recent results based on alternative methods and technologies (Chambers et al. 2004;Rietbroek et al. 2009;Wouters et al. 2011;Siegismund et al. 2011;Hughes et al. 2012;Bergmann-Wolf et al. 2014). Our results are in line with those estimates in terms of annual amplitude and phase.
The co-estimated geocentre motion is significantly different from the one derived from the degree-1 coefficients published on the Tellus website (ftp://podaac.jpl.nasa.gov/ allData/tellus/L2/degree_1), both in terms of a trend and annual amplitudes, especially for the Z-component. However, the obtained results are statistically equivalent to those published in the Tellus website when we use the same setup as Swenson et al. (2008), where a 300 km buffer zone is used to reduce the signal leakage when estimating the total ocean mass variation, but no buffer zone is considered when define the ocean function. This leaves the question of the optimal estimation of geocenter motion somewhat open. A more thorough analysis of this issue will be the subject of a separate study. The recommended solution is shown in bold

Discussion and conclusions
Our results (available at http://www.citg.tudelft.nl/c20) show that GRACE data at higher spherical harmonic degrees are capable of estimating seasonal changes in C 20 to a level comparable with SLR solutions. In fact, the uncertainty (computed as formal error from an analysis of time series) in the amplitude of the annual cycle is smaller for the GRACEbased solutions. This is an indication that our solutions may be less noisy than the SLR one, though it may also imply an underestimation of the signal not described by the fitted curve.
The main factor controlling the amplitude of the seasonal signal is the way how the problem of signal leakage in coastal areas is dealt with. Our simple approach of extending the land mask to include the first few hundreds of kilometres of coastal waters is already capable of producing a solution in close agreement with SLR, though more advanced techniques (e.g., based on mascons) could provide a better way to improve the spatial resolution of GRACE monthly fields and avoid the use of a buffer zone.
Accounting for self-attraction and loading effects driven by the redistribution of continental water masses has the effect of significantly increasing the amplitude of both annual signal and trend.
So far, we have discussed only estimates without the contribution of atmospheric and oceanic processes, assuming that the AOD1B products are correct. In the bottom line of Table 1, we list the full values determined from the SLR time series prior to the subtraction of the AOD1B signal. Compared to the GSM-like solution in the top line, the amplitude of the annual signal is twice as large and its phase is shifted by a month. This suggests that only about half of the seasonal total C 20 signal is determined by land hydrological processes, including the cryosphere. Therefore, if the proposed methodology is used in estimating the total C 20 signal, the accuracy of the obtained estimates will be dependent on the accuracy of the atmosphere-ocean model.
The determination of a long-term trend requires the use of a model of GIA, which still carries large uncertainties of an unknown magnitude. Further investigations are warranted in the future to mitigate the uncertainties introduced by a GIA model.
One need to bear in mind that the SLR solution is not free of systematic errors and noise. The processing parameters tuned to achieve a time series that best fit the SLR solution may, therefore, be biased. Further study for validation using accurate geophysical models may enable us to claim an even better solution than that from SLR.