Optimizing estimates of annual variations and trends in geocentermotion and J2 from a
combination of GRACE data and geophysical models
Sun, Yu; Riva, Riccardo; Ditmar, Pavel
DOI
10.1002/2016JB013073
Publication date
2016
Document Version
Final published version
Published in
Journal of Geophysical Research
Citation (APA)
Sun, Y., Riva, R., & Ditmar, P. (2016). Optimizing estimates of annual variations and trends in
geocentermotion and J2 from a combination of GRACE data and geophysical models. Journal of
Geophysical Research, 121(11), 8352-8370. https://doi.org/10.1002/2016JB013073
Important note
To cite this publication, please use the final published version (if applicable).
Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
First, we conduct an end-to-end numerical simulation study. We generate input time-variable gravity field observations by perturbing a synthetic Earth model with realistically simulated errors. We show that it is important to avoid large errors at short wavelengths and signal leakage from land to ocean, as well as to account for self-attraction and loading effects. Second, the optimal implementation strategy is applied to real GRACE data. We show that the estimates of annual amplitude in geocenter motion are in line with estimates from other techniques, such as satellite laser ranging (SLR) and global GPS inversion. At the same time, annual amplitudes of C10and C11are increased by about 50% and 20%, respectively, compared to estimates based on Swenson et al. (2008). Estimates of J2variations are by about 15% larger than SLR results in terms of annual amplitude. Linear trend estimates are dependent on the adopted GIA model but still comparable to some SLR results.
1. Introduction
The Gravity Recovery and Climate Experiment (GRACE) satellite mission [Tapley et al., 2004] has been mon-itoring the Earth system for more than a decade since launched in 2002. Monthly gravity field solutions produced on its basis in the form of Stokes coefficients are being released by a number of data analysis cen-ters [Bettadpur, 2012; Dahle et al., 2013; Watkins, 2012; Liu et al., 2010]. Such solutions lack the three degree 1 Stokes coefficients, which are needed for a complete representation of the mass redistribution in the Earth system. The three degree 1 coefficients are proportional to geocenter motion, here defined as the motion of center of mass (CM) of the Earth system with respect to center of figure (CF) of the solid Earth surface [Ray, 1999].
In principle, the degree 1 coefficients can be observed directly by tracking satellites, which are orbiting the CM, with ground stations, whose positions are fixed to the solid Earth’s surface. Accurate satellite tracking can be achieved by making use of geodetic techniques, such as Doppler orbitography and radiopositioning integrated by satellite (DORIS), Global Positioning System (GPS), and satellite laser ranging (SLR). However, DORIS- and GPS-based results are still contaminated by technique-specific errors and degraded by parameter collinearity [Altamimi et al., 2011; Meindl et al., 2013; Rebischung et al., 2014]. SLR solutions [e.g., Cheng et al., 2013b; So´snica et al., 2013] are among the most accurate estimates available but still have large uncertainties, as the sparse SLR tracking network is sensitive to the solid Earth deformation caused by surface mass loading, which makes it difficult to realize an accurate CF reference frame [Wu et al., 2012].
As shown by Blewitt et al. [2001], the translation of an elastic Earth surface caused by the degree 1 mass loading is accompanied by a specific deformation, which is detectable from globally distributed GPS measurements and can then be used to estimate degree 1 coefficients. This method has been further developed later by combining GPS data with GRACE data and an Ocean Bottom Pressure (OBP) model in a joint inversion scheme, which helps to reduce the aliases caused by uneven distribution of GPS sites as well as higher-degree loading [e.g., Kusche and Schrama, 2005; Wu et al., 2006; Rietbroek et al., 2009]. Nevertheless, the GPS measurements are still biased by draconitic errors [Griffiths and Ray, 2013] and solid Earth deformation unrelated to loading, such as thermal expansion of the bedrock [Dong et al., 2002] and tectonic movements.
elimination of them to a large extent
Correspondence to: Y. Sun,
y.sun-2@tudelft.nl
Citation:
Sun, Y., R. Riva, and P. Ditmar (2016), Optimizing estimates of annual variations and trends in geocenter motion and J2from a combination of GRACE data and geophysical models, J. Geophys. Res. Solid Earth, 121, doi:10.1002/2016JB013073.
Received 8 APR 2016 Accepted 3 NOV 2016
Accepted article online 8 NOV 2016
©2016. American Geophysical Union. All Rights Reserved.
Swenson et al. [2008] proposed a new methodology to estimate the degree 1 coefficients, which is based on GRACE data and an OBP model (from here on called the GRACE-OBP method). With this method, regional mass variations predicted from an OBP model are used as an additional constraint to transform the GRACE-based global mass anomalies from the CM to the CF frame. The estimates by Swenson et al. [2008], however, yielded a much smaller annual amplitude of C10variations than other approaches. Since the C10coefficient variations
correspond to the Z component of geocenter motion, inaccuracies in their determination can lead to large errors in estimated mass variations at high latitudes. For example, a 1 mm change in the Z component of geo-center position is equivalent to a 70 Gt mass change over Antarctica and 11 Gt mass change over Greenland [Wu et al., 2012]. In previous studies, simple validations of the GRACE-OBP method were performed by using simulated time-variable gravity fields based on geophysical models [Swenson et al., 2008; Bergmann-Wolf et al., 2014]. However, the purpose of those tests was limited to verifying the correctness of the theory. No effort has been made to evaluate the method performance in the context of real GRACE data.
Next to geocenter motion, the GRACE-OBP method is able to determine variations in the Earth’s dynamic oblateness (J2) [Sun et al., 2016]. Those variations are directly related to variations of the C20Stokes coefficient
(J2= −
√
5C20). The GRACE-based C20coefficient is subject to large uncertainties [Chen et al., 2016],
presum-ably due to tide-like aliases, and it is a common practice to replace it with estimates from other techniques, such as SLR. The GRACE-OBP method solves the problem by replacing the GRACE-based C20estimate with an alternative one that does not require an explicit observation of that coefficient. Unlike the estimates based on GRACE data alone, the annual variations of J2obtained with the GRACE-OBP method are comparable to independent results based on SLR observations.
In this study, for the first time, we evaluate the performance of the GRACE-OBP method with respect to the determination of both geocenter and J2variations by means of end-to-end simulations. Apart from simulating the time-variable gravity fields using synthetic models, we also take into account the errors that are present in the real GRACE data and OBP data. As a result, we determine optimal implementation parameters of the GRACE-OBP method. The primary aim of the study is twofold: (i) to optimize the methodology of Swenson et al. [2008] in order to improve estimates of annual variations and long-term trend in degree 1 and C20coefficients
and (ii) to demonstrate the impact of this optimization on estimates based on real GRACE data.
In the following, we first introduce the GRACE-OBP method and its input data (section 2). Realistic input data are then simulated with a synthetic Earth model plus realistically defined errors (section 3). Subsequently, we discuss the key implementation details that may significantly affect the results (section 4). Afterward, degree 1 and C20solutions based on both synthetic data (section 5) and real GRACE data (section 6) are computed and presented. Finally, we discuss our results in section 7.
2. Methodology and Input Data
2.1. Methodology
In this section, we intended not to reproduce the derivation of all the equations used in the GRACE-OBP method (for that, a reader is referred to Swenson et al. [2008]) but rather to bring attention to several important aspects.
Time variations of Stokes coefficients of spherical harmonic degree l, order m from GRACE data, ΔC′
lmand ΔS
′
lm, are related to the time variations of mass coefficients, ΔClmand ΔSlm, depicting the surface mass redistribution at the Earth’s surface. Let us drop, for simplicity, the symbol Δ, even though all coefficients used in this study should be understood as time variations. Then, this relation is [Wahr et al., 1998]
{ Clm Slm } = a𝜌Earth(2l + 1) 3(1 + kl) { C′ lm S′ lm } , (1)
where a is the average Earth radius,𝜌Earthis the average density of the Earth, and klis the gravitational elastic load Love number of degree l [Farrell, 1972]. As a result, Stokes coefficients can be obtained from known mass coefficients. Geocenter motion (CM with respect to CF) in Cartesian coordinates (X, Y, Z) is then related to the three Stokes coefficients through [Blewitt, 2003]
⎧ ⎪ ⎨ ⎪ ⎩ X Y Z ⎫ ⎪ ⎬ ⎪ ⎭ = √ 3a 1 + k1 ⎧ ⎪ ⎨ ⎪ ⎩ C′ 11 S′ 11 C′ 10 ⎫ ⎪ ⎬ ⎪ ⎭ , (2)
where Cocean 10 , C ocean 11 , S ocean 11 , and C ocean
20 represent the degree 1 and C20coefficients of oceanic mass changes,
respectively. The I matrix and G vector are defined as follows: I11S20C= 1
4𝜋 ∫ dΩ ̄P11(cos𝜃) sin(1 × 𝜙)𝜗(𝜃, 𝜙)̄P20(cos𝜃) cos(0 × 𝜙), (4)
G11C= 1 4𝜋 ∫ dΩ ̄P11(cos𝜃) cos(1 × 𝜙)𝜗(𝜃, 𝜙) ∞ ∑ l=2 l ∑ m=m0
̄Plm(cos𝜃){Clmcos m𝜙 + Slmsin m𝜙}
(other elements of the I matrix and G vector are similarly defined),
(5)
where integrals are defined globally with dΩ = sin𝜃d𝜃d𝜙; ̄Plmare normalized associated Legendre functions; 𝜃 is colatitude in spherical coordinates; 𝜙 is longitude; 𝜗(𝜃, 𝜙) is the ocean function, which equal to 1 over ocean and 0 over land; and m0is equal to 1 if l = 2 and 0 if l> 2. The elements in the G vector are essentially estimates of the oceanic degree 1 and C20coefficients based on remaining GRACE-based mass coefficients (note that the “oceanic coefficients” are defined as mass coefficients describing only the mass redistribution in the ocean area). Surface spherical harmonics are not orthogonal to each other when the study area is limited to a part of the sphere, and this allows some coefficients to be estimated on the basis of the mass redistribution in the study area synthesized from the remaining coefficients. Cocean
10 , C ocean 11 , S ocean 11 , and C ocean 20 are calculated
from an independent (ocean) model. Hence, the difference between the two should only be attributed to the lack of degree 1 and C20coefficients in the GRACE data, which results in an imperfect synthesis of the mass redistribution over the ocean area.
The procedure is somewhat similar to the restoration of the integration constant when the value of the primi-tive function at a certain point is known. Since the number of unknown coefficients is four, knowing the ocean mass variations at just four properly chosen data points is ,in principle, sufficient to find the unknown coeffi-cients. In practice, more data points are required to make the estimation accurate and stable. Nevertheless, the ocean function𝜗(𝜃, 𝜙) does not have to include the whole ocean: particularly, noisy regions can be excluded.
2.2. Input Data
Here we discuss the input for the GRACE-OBP method in general terms. Specific input data sets will be described later.
According to equations (3)–(5), two data sets are needed to estimate the four unknown coefficients: (i) a set of mass coefficients representing the global mass redistribution at the Earth’s surface and (ii) the oceanic com-ponent of the four unknown coefficients predicted in the CF frame. The mass coefficients are obtained from the GRACE level-2 data (labeled as GSM). The GSM coefficients are free from signals due to mass redistribution in the atmosphere and ocean for which purpose the Atmosphere and Ocean De-aliasing level-1B (AOD1B) product [Flechtner and Dobslaw, 2013] is used to clean raw GRACE data from the corresponding signals. In order to use these coefficients in the GRACE-OBP method, one should also clean the GSM coefficients from signals due to the solid Earth (e.g., glacial isostatic adjustment (GIA), megathrust earthquakes) by applying equation (1). The signals from the atmosphere and ocean (described by GAD product) should also be removed
from the oceanic degree 1 coefficients used in equation (3) to keep them compatible with vector G [Swenson et al., 2008; Sun et al., 2016].
Ideally, the residual oceanic degree 1 coefficients (Cocean
1m − G1mC, Socean11 − G11S) should only reflect the signals
not modeled by the OBP models, e.g., the mass exchange between ocean and land. The accurate estimation of ocean-land mass exchange is a delicate issue due to the fact that most ocean models, including the Ocean Model for Circulation and Tides (OMCT) [Thomas, 2002], which is used to produce the AOD1B, conserve the total volume of the ocean (Boussinesq approximation). This assumption results in artificial changes in the total mass of the ocean. The mean ocean mass is then removed (by adding/removing a uniform layer of water to the ocean) to conserve the total ocean mass. As a result, the ocean model can only predict internal mass redistribution.
Variations in the total ocean mass can be taken into account by integrating the GSM-based (i.e., GRACE-observed) mass anomalies over the oceans. Then, the corresponding oceanic coefficients are esti-mated by assuming a certain spatial distribution of this mass. Note that variations derived from the GSM product are in the CM frame and, in our case, lack the C20coefficient. To restore the missing four coefficients, it
is possible to make use of an iterative approach. The four coefficients computed by means of the GRACE-OBP method are used to complement the GSM product at the next iteration. Starting from an initial guess where those coefficients are null, the estimation of the total ocean mass variation usually converges within a few iterations.
3. Simulation of GRACE Solutions
3.1. Error-Free GSM Coefficients
The error-free GSM coefficients are constructed using the updated European Space Agency Earth System Model (ESM) [Dobslaw et al., 2015]. The updated ESM employs state-of-the-art geophysical models to simulate gravity field variations related to five components of the Earth system: atmosphere (A), ocean (O), terrestrial water (H), continental ice sheets (I), and the solid Earth (S). In our simulation, we only deal with the first four components, which are related to atmosphere and water mass redistribution at the Earth surface. In other words, we assume that the solid Earth signal can be completely removed from GRACE data by means of independent models, though it is not exactly the case when dealing with real data.
The atmosphere component is based on the latest reanalysis from ECMWF, ERA-Interim [Dee et al., 2011]. Terrestrial water storage variations are based on Land Surface Discharge Model [Dill, 2008]. The cryospheric component is constructed using the surface mass balance predictions from RACMO [Ettema et al., 2009] and a simple linear ice discharge model developed by Gruber et al. [2011] when modeling the original ESM. The oceanic part is generated by summing up contributions from OMCT, mesoscale variability from MPIOM ocean model [Storch et al., 2012], and a uniform layer of water used to conserve the mass in the Earth system. The ESM model covers a 12 year period from 1995 to 2006 with a temporal resolution of 6 h and a spatial resolution of about 100 km in terms of half wavelengths (complete to spherical harmonic degree 180). In addition, the unperturbed dealiasing model named DEAL is provided [Dobslaw et al., 2016]. DEAL is different from the sum of A and O, as it ignores the barystatic sea level variability and the small-scale ocean variability, which are also omitted in current AOD1B product. Along with the DEAL product, AOerr files are also delivered. These files contain a realization of realistic errors in A and O components. These errors need to be added to DEAL to arrive at a realistically perturbed model equivalent to the AOD1B [Dobslaw et al., 2016].
The above mentioned synthetic Earth is our study object. The “error-free” GSM coefficients are then simu-lated by subtracting the DEAL product from the result. (Notice that these coefficients are strictly speaking not error free because the DEAL product lacks some signals that are present in the A and O components.) To match the temporal resolution of the real GRACE data, monthly averages of the simulated GSM coefficients are calculated.
It is worth noting that mass conservation is enforced in ESM by adding/removing a uniform layer of water over the oceans to balance variations in the total mass of all the components. In reality, most of the water that represents mass exchange between ocean and land does not spread uniformly over oceans because the distribution of the ocean water is subjected to self-attraction and loading (SAL) effects [Gordeev et al., 1977]. In order to make our synthetic model closer to the real Earth, we modify the updated ESM model such that
Figure 1. Synthesized GRACE random errors of June 2006 in terms of equivalent water heights.
the exchanging water distributes over the oceans according to the sea level equation [Farrell and Clark, 1976; Tamisiea et al., 2010]. We expect SAL effects to have a large impact on water redistribution [Clarke et al., 2005]. 3.2. Simulation of GSM Errors
To generate realistic GSM coefficients, one needs to simulate the errors that are present in real GRACE data. In this study, we consider two types of errors: (i) errors in the GAC product and (ii) errors in the level-2 GRACE data (random errors). The GAC errors are the errors in monthly averages of the AOD1B product. They are mimicked by the errors provided in the AOerr file after we compute their monthly averages. The GRACE random errors, however, have multiple contributing sources, such as uncertainties of onboard sensors, deficiencies in orbit determination, and errors in background geophysical models (including the AOD1B product) at short time scales. An extended discussion of sources of noise in GRACE data can be found in, e.g., Ditmar et al. [2012]. Random errors in GRACE level-2 data are different from month to month, which can be explained, first of all, by variations in the orbit ground track pattern [Wagner et al., 2006]. We simulate those errors r on the basis of the error variance-covariance matrices C of GRACE Stokes coefficients as follows.
r = Lx, (6)
where the vector x contains normally distributed uncorrelated random numbers with zero mean and unit variance (X ∼ B(0, 1)). L is the lower triangle matrix obtained by the Cholesky decomposition of the error variance-covariance matrix C (C = LLT). In our study, we use the matrices C for the CSR RL05 GRACE solutions, which have been released since recently (ftp://ftp.csr.utexas.edu/outgoing/grace/). The simulated random errors for a particular month form a set of Stokes coefficients complete to the maximum degree of the error variance-covariance matrix, which is 96. For the sake of consistency, we simulate error-free GSM coefficients also to degree 96. We simulate GRACE errors for a period of approximately 12 years using the error covariance matrices from 2003 to 2014. Note that missing months in the time series of real GRACE data are not included. Then the time tags of error realizations are shifted by 8 years backward to match the updated ESM time interval (1995–2006). In Figure 1 we show the simulated random errors in level-2 GRACE data for June 2006 in terms of equivalent water heights. These errors form north-south oriented stripes in the spatial domain, which are similar to those in real level-2 GRACE data. They can be largely suppressed by appropriate filtering [e.g., Klees et al., 2008; Kusche et al., 2009]. However, in our simulation, we use these errors as they are because they largely cancel out as we integrate mass anomalies over a sufficiently large area. Using a filtered solution as input carries the risk of introducing nonnegligible biases into low-frequency gravity signals.
Finally, we add the simulated GAC and random errors to the error-free GSM coefficients from the ESM to obtain the realistically perturbed coefficients, which are used as input in the further simulations. We generated 100 realizations of random errors and thus have 100 sets of simulated noisy GSM coefficients.
According to equation (3), the degree 1 and C20estimates are affected by errors in both GSM coefficients and oceanic degree 1 and C20coefficients. Here we assume that the four oceanic coefficients are error free and try
to quantify only the impact of errors in GSM coefficients. We will address the impact of errors in the four input oceanic coefficients in section 6.3.
4. Implementation Details
Apart from the input GSM coefficients and geophysical models, the solution of equation (3) is largely depen-dent on three implementation details: (i) the truncation degree of the GSM coefficients, (ii) the ocean function that partitions the Earth’s surface into ocean and land, and (iii) whether or not considering SAL effects when distributing ocean water.
We will have to find the truncation degree that is low enough to exclude the strong noise present at higher degrees but still high enough to depict the large-scale mass redistribution, which contributes significantly to the solutions.
The second implementation parameter is the ocean function, which is responsible for selecting ocean grid points that are used to constrain the solution. Due to limited resolution of GRACE data, strong continental sig-nals leak into ocean areas, so that an ocean function using the exact ocean-land boundaries will erroneously capture them as ocean signals. In order to reduce leakage from land into the oceans, we establish a buffer zone by excluding points within a few hundred kilometers (up to 400 km) from the coast. A reasonable buffer width choice will allow us to use an ocean area not polluted by leakage from land and still large enough to ensure a stable solution.
In order to use the GRACE-OBP method, we need to define the four input oceanic coefficients (Cocean 10 , … ,C
ocean 20 ;
see equation (3)). As explained in section 2.2, these coefficients should only reflect the mass variations over ocean areas resulting from the water exchange between the continents and the oceans. After fixing the trun-cation degree of GSM coefficients and the width of the buffer zone, the estimated total water exchange is also fixed. However, the solution is still affected by how water distributes over the oceans. Here we consider two scenarios as follows: (i) ignoring the SAL effects, so that water redistributes uniformly over oceans as in Swenson et al. [2008]; the resulting ocean height change is known as barystatic (also called eustatic) sea level variation, and (ii) taking into account SAL effects by computing the resulting fingerprints [Mitrovica et al., 2001] using the sea level equation.
In the following section, we analyze the effect of different choices of the implementation details. We discuss degree 1 and C20coefficients separately.
5. Results of Numerical Experiments
In the numerical experiments, we estimate degree 1 and C20time series based on different combinations of the implementation parameters. The truncation degrees tested for the input “GRACE” solutions are from 10 to 96 with a 1∘ step. We show results for degrees 10 to 70 because solutions based on higher truncation degrees are in all cases too noisy due to large errors in high-degree coefficients. For each truncation degree, we test buffer widths from 0 to 400 km with a 100 km increment. For a particular combination of truncation degree and buffer zone width, we consider water distributions with and without taking SAL effects into account. Different degree 1 and C20time series (TC10, TC11, TS11, and TC20) corresponding to specific parameter combinations are
obtained. Amplitude (Tamp), phase (Tpha), and linear trend (Ttrnd) of the annual signal are estimated. Note that
Tampand Tphaare defined by approximating the annual cycle with the expression Tampcos(𝜔(t − t
0) − Tpha),
where𝜔 is equal to 2𝜋 rad/yr and t0refers to 1 January of a particular year.
The resulting degree 1 and C20time series vary significantly for different synthetic GSM solutions reflecting
different error realizations, which means that results based on just one set of GSM coefficients are not rep-resentative. Therefore, we show results based on all 100 sets of simulated GSM solutions. That is, for each combination of implementation parameters, we obtain 100 estimates of annual amplitude, phase, and linear
Figure 2. Results of simulated data processing: The mean annual amplitudes of the estimated GSM degree 1 time series
(in millimeter of geocenter motion) using different implementation parameters based on 100 sets of simulated GSM solutions. Their standard deviations are indicated by light colored bands. The true amplitudes are marked in all panels as black horizontal lines. (a–c) The amplitude estimates are based on the uniform exchanging water distribution (NoSAL). (d–f ) Results after distributing the exchanged water according to its gravitational fingerprints (SAL).
trend. Then, we show the mean value of the 100 estimates and use the corresponding standard deviation as the uncertainty range.
5.1. Annual Variations of Geocenter Motion
In Figure 2, we show mean annual amplitude estimates of degree 1 time series and their standard deviations (indicated by light colored bands). The annual amplitude estimates are sensitive to all three implementation parameters. With a truncation at around degree 45, we notice that applying a 200 km buffer zone increases the amplitude estimates of C10and C11by about 40% and 10%, respectively, compared to the solutions using no buffer zone (Figures 2a and 2b). A further 25% increase is found for both coefficients after using the realistic exchanging water distribution (Figures 2d and 2e). In contrast, accounting for SAL effects and using a buffer zone have contradicting effects on the estimated amplitude of S11time series. Increasing the buffer tends
to reduce further the already underestimated annual amplitude (Figure 2c). Considering SAL effects, on the other hand, increases the estimates, so that the resulting amplitude is close to the synthetic truth (Figure 2f ). Generally speaking, the figure shows that without taking into account the SAL effects in the data processing procedure, one cannot retrieve the true amplitude, no matter how the truncation degree and buffer width are chosen (Figures 2a–2c). However, by considering only the SAL effects without properly choosing the buffer zone width (e.g., using zero buffer width), one still cannot match the true amplitude (Figures 2d–2f ). Note that the results based on the error-free GSM solutions are not shown because they are very close to the mean amplitudes. This is because the monthly averages of the AOD1B errors are negligible compared to GRACE ran-dom errors, as will be shown below. Therefore, taking the mean of 100 noisy realizations reduces the noise level by approximately a factor of 10 (√100). Best amplitude estimates are obtained after accounting for SAL effects and applying a buffer zone wider than 100 km. With these two implementation parameters chosen, the truncation degree becomes a less decisive factor, and a broad range of values (from 30 to 60) allows for
Figure 3. Results of simulated data processing: (a–c) The mean annual phases of the estimated GSM degree 1 time
series including their standard deviation (indicated by light colored bands) using different implementation parameters based on 100 sets of simulated GSM solutions. The true phases are marked in all panels as black horizontal lines.
estimating the amplitude accurately. Truncation degree higher than 60 clearly results in an increased uncer-tainty because of large errors in the high-degree coefficients.
Similar conclusions can be drawn for the C20coefficient, as well as for the phases of annual variations and linear trends. In the following figures, we limit ourselves to presenting solutions using 0 km buffer width and ignoring SAL effects to show how the results look like if the original methodology of Swenson et al. [2008] is reproduced. In addition, we show solutions based on a 200 km buffer zone with and without considering SAL effects to illustrate the impact of a properly chosen buffer zone and considering SAL effects.
In Figure 3, we show the annual phase estimates for degree 1 coefficients. For all three coefficients, a good estimation of the phase can be obtained by accounting for the SAL effects and using the buffer width of 200 km. The buffer zone width does not play a critical role here as the use of a 200 km buffer zone barely alters the phase estimates. Considering SAL effects has a much larger effect (up to 10 days), especially for the cases of C11and S11.
5.2. Linear Trend in Geocenter Motion
In Figure 4, we show the effect of implementation details on estimates of a long-term linear trend. An accurate estimation of C10trend also requires the use of a buffer zone and the inclusion of SAL effects. In contrast, the
impact of implementation details onto the estimation of C11is minor. Finally, we see a better estimation of the S11trend by using a buffer zone and taking into account SAL effects, but the uncertainties are relatively
large, which implies that the trend can still be overestimated or underestimated by up to 25%.
Figure 5. Similar to Figure 3 but for annual amplitude, annual phase, and linear trend of theC20coefficient.
From this test, it is tempting to conclude that the GRACE-OBP method is able to retrieve the correct linear trend, especially in the C10and C11coefficients. However, the synthetic model is not affected by the solid
Earth contributions. Uncertainties in modeling those contributions might dominate in trend estimates when working with real data. We will discuss this issue further in section 6.
5.3. C20Variations
As mentioned earlier, the GRACE-OBP method has been already used to determine J2variations (or
equiva-lently C20variations) by Sun et al. [2016], who validated the results against estimates from SLR over
approxi-mately the same interval as in our study. The optimal strategy found in that study consisted of using GRACE monthly solutions complete to spherical harmonic degree 60 and filtered by means of the DDK4 filter [Kusche et al., 2009], in combination with a buffer width of 150 km, and accounting for SAL effects. However, a specific SLR solution [Cheng et al., 2013a] may not represent the truth and that strategy may offer a biased parameter setting. Here we make use of the synthetic model for an independent validation of the GRACE-OBP method as well as for choosing its optimal implementation. Results are shown in Figure 5.
As far as the amplitude of the annual signal is concerned, the use of a buffer zone of 200 km has a rela-tively small effect of increasing the estimated value, which remains 50% smaller than the truth (Figure 5a). Accounting for SAL effects (Figure 5b) doubles the amplitude and allows us to recover the truth almost exactly (in combination with a 200 km buffer and a truncation degree between 30 and 50). The phase of the annual signal, on the other hand, is fairly independent from the buffer size and SAL effects and varies within not more than 1 week for truncation degree between 30 and 50 (Figure 5b).
The linear trend estimate turns out to be highly dependent on all three implementation details. Including a 200 km buffer zone and accounting for SAL effects have again a large impact. At the same time, it is also important to choose a truncation at around degree 45 to obtain the true values.
5.4. Optimal Implementation Parameter Setup
In this section, we apply a formal criterion to identify the optimal parameter setup for estimating degree 1 and C20coefficients simultaneously with the GRACE-OBP method. We define the quality indicator for annual
variations (Qann) and linear trend (Qtrnd) separately as follows:
Qann= C20 ∑ cf =C10 n ∑ i=1 {(
Tcf,iampsin(Tcf,ipha)− Tcf(t)ampsin(Tcf(t)pha)
2)
+ (
Tcf,iampcos(Tcf,ipha)− Tcf(t)ampcos(Tcf(t)pha)
2)} (7) and Qtrnd= C20 ∑ cf =C10 n ∑ i=1 {( Ttrnd cf,i − Tcf(t)trnd )2} , (8)
where Tcfrepresents the time series of the coefficient indicated by index cf and the index of cf runs over the four estimated coefficients, namely, C10, C11, S11, and C20. Note that these coefficients are scaled to describe
mass changes in terms of equivalent water height (equations (1) and (2)). Subscript i indicates that the result is based on the ith realization of noisy GSM coefficients; n is the number of simulated GSM solutions,
Figure 6. Results of simulated data processing: The quality indicator for different combinations of implementation
parameters. We show the quality of estimated (a) annual variations and (b) linear trends.
which is equal to 100. Superscripts amp, pha, and trnd represent the annual amplitude, phase, and trend of the corresponding time series, respectively; superscript (t) refers to the synthetic truth. We are looking for the parameter settings that yield the smallest Qannand Qtrndvalue.
In Figure 6, we show the quality of selected parameter setups. The quality of the estimated annual variations and linear trends of the four coefficients is much lower if SAL effects are not taken into account. Therefore, that option is not addressed in the figure. Clearly, the best results are obtained from the buffer width from 150 to 250 km in combination with a proper truncation. Ultimately, we select the option “T45 BUF200 SAL,” which uses a truncation at degree 45, buffer width of 200 km, and accounting for SAL effects, as the recommended parameter setup for estimating both annual variations and linear trend.
5.5. The Impact of GSM Errors
It is clear that both random errors in GRACE data and errors in the monthly averages of atmosphere and OBP (GAC) contribute to the uncertainties in the obtained GSM coefficients. It is interesting to quantify the rel-ative contributions of these two types of errors to the error budget of the obtained estimates of degree 1 and C20coefficients. Here after setting the implementation parameters to the identified optimal setup, we
consider the following: (i) GSM coefficients only subject to GAC errors and (ii) those only subject to GRACE random errors. Both are compared with results based on the error-free GSM coefficients. As shown in Table 1, the uncertainty due to GAC errors are less than 10% of that due to GRACE random errors for all four coeffi-cients. Compared to the impact of GRACE random errors, the impact of GAC errors is negligible. Therefore, it is enough to just quantify the impact of GRACE random errors (with 100 realizations) on the estimates of annual variations (the impact of GAC errors cannot be shown realistically since we have only one realization for that type of error). These findings are likely applicable to real data as well, since the provided AOerr represent errors in AOD1B product used in the GRACE data processing [Dobslaw et al., 2016].
Table 1. Results of Simulated Data Processing: RMS Errors (RMS Values of the Difference Between Obtained Time Series
and the Synthetic Truth) of the Degree 1 andC20Solutions Due to Atmosphere and OBP Errors and GRACE Random Errorsa
C10(mm) C11(mm) S11(mm) C20(10−11)
RMS error GAC 0.04 0.02 0.02 0.01
GRACE 0.49 0.50 0.56 0.13
Ann Var error Amp 0.1 0.1 0.1 0.2
Pha (day) 2 2 3 2
aWe also show the impact of RMS errors on the estimation of annual variations. RMS errors as well as errors in annual amplitude of degree 1 solutions are reported in terms of geocenter motion.
Figure 7. Results of real data processing: The annual amplitudes of the estimated GSM degree 1 time series
(in millimeter of geocenter motion) using different implementation parameters based on GRACE data (CSR RL05). The uncertainty ranges are estimated from least squares fitting error. (d–f ) The values obtained with the recommended data processing setup are marked with blue dots.
6. Results Based On Real GRACE Data
In this section, we present degree 1 and C20solutions based on the GSM coefficients provided by CSR RL05 (complete to degree 96) covering a 12 year period (from August 2002 to June 2014). The pole tide (mainly affecting the C21and S21coefficients) has been corrected for according to Wahr et al. [2015]. Unlike the syn-thetic input, real GSM coefficients are not free of time-variable signals from the interior of the Earth. Here we only attempt to remove the prominent signals caused by GIA. Since the use of a GIA model does not affect annual variations, we postpone the GIA discussion until we show the linear trend estimates.
6.1. Geocenter Motion
Even though we have already identified the optimal parameter setup, we find it more informative to pro-duce degree 1 solutions for different combinations of implementation details as it was done in the numerical experiments. Since the results are now based on just one set of real GSM solutions, the uncertainty ranges are computed based on least squares fitting error only.
When looking at the estimated annual amplitudes (Figure 7), we see that the effects of changing the buffer size and including SAL effects are extremely similar to the synthetic case (Figure 2). Consequently, we believe that the recommended parameter setup based on the numerical experiments can also be applied to real data, even though the synthetic Earth employed in the numerical studies may not model the Earth system perfectly. The impact of the implementation parameters on the annual phase (Figure 8) is similar to the synthetic case (Figure 3) only for C11, while C10and S11show a larger dependence on the buffer size and a smaller dependence on taking SAL effects into account. Nonetheless, the qualitative agreement between the results obtained with synthetic and real data still supports the use of a buffer and accounting for SAL effects.
So far, we have discussed only degree 1 time series without the contribution of atmospheric and oceanic processes (GSM-like coefficients). In Table 2, we show the results in terms of full geocenter motion, i.e., after
Figure 8. Results of real data processing: (a–c) The annual phases of the estimated GSM degree 1 time series using
different implementation parameters. The uncertainty ranges are calculated based on least squares fitting errors.
restoring the degree 1 coefficients of the GAC products. In order to quantify the GAC-modeled atmosphere and ocean contribution, we additionally list the results for GSM-like coefficients as a reference. The ampli-tude and phase of C10are changed after restoring the GAC product insignificantly: by less than 1% and by
about half a month, respectively. In contrast, C11and S11are largely affected. The GAC products change their
phase estimates by about a month and account for about 30% and 40% of the total amplitude of C11and S11,
respectively. Table 2 also shows the solutions from the GRACE TELLUS website that are based on the methodology by Swenson et al. [2008] (ftp://podaac.jpl.nasa.gov/allData/tellus/L2/degree_1/, downloaded in December 2015; GAC product is restored), as well as other solutions from recent literature. Compared to the results based on Swenson et al. [2008], our estimates have a considerably larger annual amplitude for the C10
and C11coefficients, respectively, by about 50% and 20%. Since the same GAC product is restored to both
solu-tions, the increase can only be caused by a different setup of the implementation parameters. According to Figure 7, 54% of the increase in C10amplitude is due to the use of a 200 km ocean buffer, and the rest is due to
Table 2. Results of Simulated Data Processing: Estimated Amplitudes and Phases of the Annual Geocenter Motiona
C10 C11 S11
Amp (mm) Pha (day) Amp (mm) Pha (day) Amp (mm) Pha (day) Time Span (year)
GSM T45 BUF200 SAL OMCT 2.8 ± 0.1 85 ± 2 1.6 ± 0.1 85 ± 3 1.7 ± 0.1 302 ± 3 2002.6–2014.5
GSM T45 BUF200 SAL ECCO_Polar_zero 2.6 ± 0.1 79 ± 2 1.6 ± 0.1 75 ± 3 1.8 ± 0.1 309 ± 2 2002.6–2014.5 GSM T45 BUF200 SAL ECCO_Polar_OMCT 2.7 ± 0.1 82 ± 2 1.6 ± 0.1 82 ± 3 1.7 ± 0.1 306 ± 2 2002.6–2014.5
CSR T45 BUF200 SAL OMCT 2.9 ± 0.2 69 ± 4 2.3 ± 0.1 52 ± 3 2.8 ± 0.1 327 ± 2 2002.6–2014.5
CSR T45 BUF200 SAL ECCO_Polar_zero 2.5 ± 0.2 58 ± 4 2.0 ± 0.1 53 ± 3 3.0 ± 0.1 333 ± 2 2002.6–2014.5 CSR T45 BUF200 SAL ECCO_Polar_OMCT 2.7 ± 0.2 59 ± 4 2.0 ± 0.1 54 ± 3 3.0 ± 0.1 334 ± 2 2002.6–2014.5
GFZ T45 BUF200 SAL OMCT 3.0 ± 0.2 69 ± 4 2.3 ± 0.1 54 ± 3 2.8 ± 0.1 326 ± 2 2002.6–2014.5
JPL T45 BUF200 SAL OMCT 2.9 ± 0.2 72 ± 4 2.2 ± 0.1 55 ± 4 2.7 ± 0.1 328 ± 2 2002.6–2014.5
Swenson et al. [2008] 1.9 ± 0.1 65 ± 4 1.9 ± 0.1 53 ± 3 2.5 ± 0.1 319 ± 2 2002.6–2014.5
SLR [So´snica et al., 2013] 3.0 ± 0.3 - 2.75 ± 0.15 - 2.2 ± 0.1 - 2002.6–2014.5
SLR [Cheng et al., 2013b] 4.2 ± 0.3 33 ± 2 2.9 ± 0.4 49 ± 4 2.7 ± 0.1 324 ± 2 2002.6–2014.5
INV [Rietbroek et al., 2012a] 2.6 25 2.1 60 3.4 326 2003.0–2009.0
KAL-1 [Wu et al., 2015] 3.9 ± 0.1 21 ± 1 2.1 ± 0.1 45 ± 1 2.7 ± 0.1 321 ± 1 2002.3–2009.3
KAL-2 [Wu et al., 2015] 3.3 ± 0.1 22 ± 3 1.9 ± 0.1 54 ± 2 2.6 ± 0.1 322 ± 1 2002.3–2009.3
KAL-3 [Wu et al., 2015] 3.5 ± 0.1 19 ± 1 1.9 ± 0.1 52 ± 1 3.0 ± 0.1 337 ± 1 2002.3–2009.3
aThe contribution of atmosphere and ocean (GAC) is restored. As a reference, the row with leading ’GSM’ represents the optimal GSM-like geocenter motion based on CSR RL05.The top part of the table presentes results of this study, while lower part show results from other studies. T45 means a truncation at degree-45, BUFxxx denotes the buffer width of xxx km; BUF0 indicates the absence of a buffer zone. SAL indicates that the SAL effects have been taken into consideration. The same nomenclature is also used in other tables in this study. Solutions labeld with ’GFZ’ and ’JPL’ are based on GSM coefficients from GFZ RL05a and JPL RL05. OMCT and ECCO indicate the OBP model addopted in the solution.
Figure 9. Results of real data processing: The degree 1 time series based on the proposed implementation parameter
setup after the restoration of the atmosphere and ocean signal using the GAC product, together with the solution based on Swenson et al. [2008] and an SLR solution [Cheng et al., 2013b]. The uncertainties of the SLR solution are shown as gray bands.
the consideration of SAL effects. The increase in C11amplitude, on the other hand, should only be attributed
to SAL effects. The new estimates compare much better with the annual signals detected by other techniques (Table 2). Full degree 1 time series based on alternative GRACE monthly solutions, such as JPL RL05 [Watkins, 2012] and GFZ RL05a [Dahle et al., 2013], give very close estimates.
Table 3. Results of Real Data Processing: Geocenter Motion Trends Due to Surface Mass Transport Estimated Using
Different GIA Models, Based on Solution CSR T45 BUF200 SAL OMCTa
C10(mm yr−1) C
11(mm yr−1) S11(mm yr−1) Time Span (year)
ICE-5G_VM2 −0.21 ± 0.04 −0.03 ± 0.03 0.11 ± 0.02 2002.6–2014.5
ICE-6G_VM5a −0.33 ± 0.04 −0.06 ± 0.03 0.07 ± 0.02 2002.6–2014.5
Swenson et al. [2008] −0.11 ± 0.03 −0.08 ± 0.02 −0.02 ± 0.02 2002.6–2014.5
Wu et al. [2010] −0.16 ± 0.07 −0.10 ± 0.01 0.29 ± 0.05 2002.3–2009.0
Rietbroek et al. [2012b] −0.37 −0.14 0.12 2003.0–2009.0
Figure 10. Same as Figure 8 but for linear trend estimates.
In Figure 9, we plot the proposed degree 1 solution together with the solution by Swenson et al. [2008], as well as a selected SLR solution [Cheng et al., 2013b], linear trends being removed. The annual cycle is the most prominent feature in all solutions. Although we have increased the annual amplitude significantly for C10and
moderately for C11, as compared to Swenson et al. [2008], large discrepancies still exist with respect to the SLR
solution from Cheng et al. [2013b].
Since the GRACE-OBP methodology implies that mass transport takes place only at the Earth’s surface, we subtracted GIA signals from the input GSM coefficients. In Table 3, we list linear trend estimates obtained with two GIA models and compare them with trend estimates from literature. GIA model ICE-5G_VM2 [A et al., 2012] is based on the ICE-5G ice history and a simplified version of mantle viscosity model VM2 [Peltier, 2004] and computed for a compressible Earth model. GIA model ICE-6G_VM5a is based on the ICE-6G ice history and the VM5a viscosity profile [Peltier et al., 2015]. Differences between the surface mass trend estimates based on two GIA models suggest that uncertainties in GIA models play an important role in geocenter trend estimates. However, the optimal estimation of geocenter motion trend caused by surface mass transport, as well as of the full trend (which can be obtained by adding back the GIA-induced trend), is still under investigation and not discussed in this study.
In Figure 10, we demonstrate that the linear trend estimates based on various combinations of implementa-tion parameters also show a similar behavior, as compared to that observed in our numerical experiments. For example, the linear trend estimates in C10are smaller for higher truncation degrees and more narrow buffer
zones. Estimates in S11increase for wider buffer zones while are relatively insensitive to truncation degrees
larger than 30.
6.2. C20Variations
Results for C20variations based on real data are shown in Figure 11. As for the case of degree 1 coeffi-cients, both the annual signal and the linear trend show a similar behavior to the synthetic test, in terms of
Figure 12. Results of real data processing: TheC20time series based on the proposed implementation parameter setup together with the solution from Sun et al. [2016] and an SLR solution (uncertainties are indicated by light gray band) [Cheng et al., 2013a]. The contribution of atmosphere and ocean is not restored in the GRACE-OBP solutions, as the same contribution is removed from the SLR solution. A linear trend has been removed.
dependence on the truncation degree, as well as the impact of buffer width and SAL effects. Hence, we expect the best estimation of the annual amplitude to be obtained when using a 200 km wide buffer and accounting for SAL effects, where the latter is particularly important.
In Figure 12 we show the detrended GSM-like C20time series, supplemented by the SLR solution (GSM like) by Cheng et al. [2013a] and the GRACE-OBP solution by Sun et al. [2016], which was optimized to match the SLR results. Note that the optimal GRACE-OBP solution from this study generally falls within the SLR uncertainty range but is less volatile. However, the amplitude estimate over the whole time interval is considerably larger than the estimate from the SLR solution (Table 4). Comparisons with annual amplitude estimates of other SLR solutions [Lemoine et al., 2013; So´snica et al., 2014] are even worse. Further investigation is needed to understand these discrepancies.
As far as the linear trend is concerned, the results shown in Figure 11c again support the use of the same imple-mentation parameters as the synthetic test prompts. Statistics for selected solutions are shown in Table 5, where we have adopted the same GIA models as used for the degree 1 solutions. Note that to facilitate
Table 4. Results of Real Data Processing: Estimated Annual Amplitudes and Phases of theC20Time Seriesa
GSM Like Full
Amp (10−11) Pha (day) Amp (10−11) Pha (day) Time Span (year)
CSR T45 BUF200 SAL OMCT 8.6 ± 0.3 73 ± 2 16.1 ± 0.7 52 ± 3 2002.6–2014.5
CSR T45 BUF200 SAL ECCO_Polar_zero 8.0 ± 0.3 71 ± 2 12.8 ± 0.7 55 ± 3 2002.6–2014.5 CSR T45 BUF200 SAL ECCO_Polar_OMCT 8.2 ± 0.3 74 ± 2 13.8 ± 0.7 55 ± 3 2002.6–2014.5
GFZ T45 BUF200 SAL OMCT 8.5 ± 0.3 71 ± 2 16.1 ± 0.7 50 ± 3 2002.6–2014.5
JPL T45 BUF200 SAL OMCT 8.7 ± 0.3 74 ± 2 16.1 ± 0.7 52 ± 3 2002.6–2014.5
GRACE-OBP [Sun et al., 2016] 6.8 ± 0.3 77 ± 2 - - 2003.0–2013.4
SLR [Cheng et al., 2013b] 7.0 ± 0.4 81 ± 3 14.1 ± 0.7 53 ± 3 2002.6–2014.5 SLR [Lemoine et al., 2013] 6.0 ± 0.4 62 ± 3 14.0 ± 0.7 43 ± 3 2002.6–2014.5 SLR [So´snica et al., 2014] 4.5 ± 0.5 88 ± 7 11.6 ± 0.8 49 ± 4 2002.6–2014.5 aAmplitude and phase of the annual variations in all listed SLR solutions are estimated by ourselves. Results in the first two columns are free of atmosphere-ocean contributions (GSM like). These contributions are restored in the results shown in the next two columns (full).
Table 5. Results of Real Data Processing: Estimated Linear Trends FromJ2Time Series Based on CSR T45 BUF200 SAL OMCT Solution and Two Different GIA Models, as Well as From the Literaturea
Surface (10−11yr−1) GIA (10−11yr−1) Total (10−11yr−1) Time Span (year)
ICE5G_VM2 7.7 ± 0.4 −3.2 4.5 ± 0.4 2002.6–2014.5
ICE6G_VM5a 7.3 ± 0.4 −3.5 3.8 ± 0.4 2002.6–2014.5
GRACE-OBP [Sun et al., 2016] 7.4 −3.3 4.1 ± 0.2 2003.0–2013.4
SLR [Cheng et al., 2013b] - - 2.0 ± 0.3 2002.6–2014.5
SLR [Lemoine et al., 2013] - - 3.1 ± 0.4 2002.6–2014.5
SLR [So´snica et al., 2014] - - 4.0 ± 0.4 2002.6–2014.5
aSLR trends are estimated for the same time interval as considered in our study. The contribution of atmosphere and ocean has been restored in the GRACE surface mass estimates (surface). The total GRACE-based trends are obtained as the sum of surface- and GIA-related trends.
comparison with other studies, we list ̇J2values. Also, since J2varies nonlinearly, it is only meaningful to compare linear trends for approximately the same time span. Therefore, our results are only compared with several SLR solutions with approximately the same time interval. The full linear trend estimates from the GRACE-OBP approach supported with either GIA model are close to the estimate from the SLR solution by So´snica et al. [2014] but about twice as large as that from the SLR solution by Cheng et al. [2013a]. The full linear trend of the SLR solution by Lemoine et al. [2013] is relatively close to the GRACE-OBP solution based on the ICE-6G_VM5a model.
6.3. The Impact of Errors in the Oceanic Degree 1 and C20Coefficients
Until now, we have analyzed the GRACE-OBP approach under the assumption that the oceanic degree 1 and C20coefficients are error free. In reality, the four oceanic coefficients will likely contain systematic errors that
could be magnified by the GRACE-OBP approach. Through analytical error propagation based on equation (3), a 1 mm error in Cocean 10 , C ocean 11 , S ocean 11 , and C ocean
20 propagates to about 1.8, 1.4, 1.6, and 1.9 mm error in the
resulting C10, C11, S11, and C20coefficients, respectively. It is difficult to synthesize realistic error realizations for low-degree oceanic coefficients as they may be contaminated by time-correlated errors. Therefore, we limit ourselves to comparing annual variations in degree 1 and C20coefficients based on two alternative OBP
models in order to show the potential impact of such errors.
That is, we repeat the computations presented above, having replaced the OMCT model with the Estimating the Circulation and Climate of the Ocean (ECCO) [Fukumori, 2002; Kim et al., 2007]. The ECCO OBP field is from ECCO’s Near Real-Time Kalman filter estimate version kf080, which assimilates altimetry data as well as in situ temperature profiles. The ECCO OBP fields are available in monthly averages, on 1∘ × 1∘ grids covering oceans between 80∘S and 80∘N (http://grace.jpl.nasa.gov/data/get-data/ocean-bottom-pressure/). We fill the polar gaps in the ECCO OBP estimates in two ways: (i) with zeros (ECCO_Polar_zero) and (ii) with data points from OMCT OBP (ECCO_Polar_OMCT).
In Tables 2 and 4, we show the amplitude and phase estimates of the annual variations for degree 1 and C20coefficients based on different ECCO OBP models. The annual variations of all four coefficients are not significantly changed. For GSM-like degree 1 and C20coefficients, the largest differences are seen when
ECCO_Polar_zero model is used instead of OMCT. However, the annual amplitudes are only different within 10% and annual phases are different within 10 days.
When considering full coefficients, the impact of using a different OBP model becomes larger. Again, the largest difference is seen between OMCT and ECCO_Polar_zero models. Nevertheless, even in that case, the annual amplitudes of degree 1 coefficients are different within 15%. For annual phase estimate, the largest difference is 11 days. Annual amplitude of full C20coefficient is about 20% smaller when ECCO_Polar_zero is
used. The annual phases of full C20coefficients, on the other hand, agree within 1 sigma (3 days).
It is also important to realize that the optimal choice for implementation parameter based on the numerical study ignores the errors in the oceanic degree 1 and C20coefficients. Nevertheless, they are still valid in case
of estimating annual variations, judging from the relatively small differences in those estimates when using different OBP models.
cients analyzed and for different quantities of interest: the amplitude/phase of annual variations or the linear trend. Nevertheless, we find it essential to identify a single optimal setup in order to ensure a consistency of the obtained estimates. The recommended setup consists of truncation at degree 45, a buffer width of 200 km, and accounting for SAL effects.
We have also shown that in real GRACE data processing, the dependence of the solution on the implementa-tion details is similar to the synthetic case. Therefore, the optimal implementaimplementa-tion setup discussed above is likely also suitable to deal with real data. The resulting time series of geocenter motion and J2variations are expected to be significantly improved with respect to previous results based on the original methodology of Swenson et al. [2008]. Notable differences are in the amplitude of the annual signal of GSM-like C10, which is
50% larger than in Swenson et al. [2008], and in the annual amplitude of full J2, which is 15% larger than in Cheng et al. [2013a].
When we restore the atmosphere-ocean contribution to arrive at the full coefficients, the annual amplitude of C10is barely changed and the phase is shifted by only 2 weeks, which suggests that continental
hydro-logical processes, changes in the cryosphere and total ocean mass variations are responsible for most of the seasonal variations in the Z component of geocenter motion. In contrast, the annual amplitudes of C11and
S11are increased by about 30% and 70%, respectively, which implies that atmosphere-ocean variations are
largely driving seasonal variations in geocenter motion at the X and Y components. Also, J2variations are largely affected by the atmosphere-ocean contribution, which almost doubles the annual amplitude, while introducing a phase shift of 3 weeks.
Errors in both GRACE coefficients and the oceanic degree 1 and C20coefficients contribute to the uncertainties
in the obtained degree 1 and C20time series. However, GRACE errors only account for annual amplitude errors in geocenter motion at the 0.1 mm level, which is below 10% of the total signal. For C20coefficients, GRACE
errors also contribute marginally (with less than 10%). The effect of this error on to the annual phase estimates is also minor (less than 3 days). In contrast, errors in oceanic degree 1 and C20coefficients seem to play a somewhat larger role. By comparing estimates based on different OBP models, we conclude that the errors in oceanic degree 1 coefficients account for up to 15% of the annual amplitude estimates. While errors in oceanic C20coefficients could cause a difference of about 20%. Again, annual phase estimates are not significantly
affected (for all cases, the differences are well within 2 weeks).
Regarding trends in geocenter motion, we only calculate those driven by surface mass. GIA-induced trends are not reliably predicted by GIA models and are not restored. Moreover, Swenson et al. [2008] believe that geocenter motion trends caused by GIA are of negligible magnitude and do not have to be restored at all. As for the case of full ̇J2, trend estimates are highly dependent on the adopted GIA model [Sun et al., 2016]. GIA has a direct effect in terms of solid earth contribution, but it also has an important secondary effect com-ing from the use of GRACE data to constrain surface mass redistribution at higher degrees (the GRACE-OBP method establishes a link between low- and high-degree gravity variations). The GIA effect on ̇J2estimates clearly appears from Table 5, where the difference between the ̇J2 contribution of the two GIA models
(direct effect) is about 10% (0.3⋅10−11 yr−1), whereas the difference in the derived surface contribution
(indirect effect) is even larger (0.4 ⋅ 10−11 yr−1). Finally, the estimates of total ̇J
2 differ by as much as
For both geocenter and J2trends estimated with the GRACE-OBP method, a wider range of GIA models needs to be analyzed in order to draw more definitive conclusions regarding the optimal data processing scheme and its accuracy.
In this study, the discussion has concentrated on annual variations and the linear trends, which are two promi-nent features of the geocenter motion and J2time series. However, a consideration of them alone does not guarantee that the GRACE-OBP method provides the optimal estimates at other time scales (or frequencies). To quantify the overall quality of the resulting time series in a numerical study, it would be necessary to calcu-late the RMS difference between the synthetic truth and the obtained time series. However, unlike the errors related to annual cycle and linear trend, the RMS error does not reduce in case of a relatively long time series. As a result, one needs to pay more attention to mitigating the high-frequency noise in the resulting geocenter and J2variations. This will be subject of further studies.
References
A, G., J. Wahr, and S. Zhong (2012), Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: An application to Glacial Isostatic Adjustment in Antarctica and Canada, Geophys. J. Int., 192(2), 557–572, doi:10.1093/gji/ggs030.
Altamimi, Z., X. Collilieux, and L. Métivier (2011), ITRF2008: An improved solution of the international terrestrial reference frame, J. Geod., 85(8), 457–473, doi:10.1007/s00190-011-0444-4.
Bergmann-Wolf, I., L. Zhang, and H. Dobslaw (2014), Global eustatic sea-level variations for the approximation of geocenter motion from Grace, J. Geod. Sci., 4(1), 37–48, doi:10.2478/jogs-2014-0006.
Bettadpur, S. (2012), UTCSR level-2 processing standards document, Tech. Version 4, Univ. Texas, Austin.
Blewitt, G. (2003), Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth, J. Geophys. Res., 108(B2), 2103, doi:10.1029/2002JB002082.
Blewitt, G., D. Lavallée, P. Clarke, and K. Nurutdinov (2001), A new global mode of Earth deformation: Seasonal cycle detected, Science, 294(5550), 2342–2345, doi:10.1126/science.1065328.
Chen, J., C. Wilson, and J. Ries (2016), Broad-band assessment of degree-2 gravitational changes from GRACE and other estimates, 2002–2015, J. Geophys. Res. Solid Earth, 121, 2112–2128, doi:10.1002/2015JB012708.
Cheng, M., B. D. Tapley, and J. C. Ries (2013a), Deceleration in the Earth’s oblateness, J. Geophys. Res. Solid Earth, 118, 740–747, doi:10.1002/jgrb.50058.
Cheng, M. K., J. C. Ries, and B. D. Tapley (2013b), Geocenter variations from analysis of SLR data, in Reference Frames for Applications in Geosciences, Int. Assoc. of Geod. Symp., vol. 138, edited by Z. Altamimi and X. Collilieux, pp. 19–25, Springer, Berlin.
Clarke, P. J., D. A. Lavallée, G. Blewitt, T. M. van Dam, and J. M. Wahr (2005), Effect of gravitational consistency and mass conservation on seasonal surface mass loading models, Geophys. Res. Lett., 32, L08306, doi:10.1029/2005GL022441.
Dahle, C., F. Flechtner, C. Gruber, D. König, R. König, G. Michalak, K.-H. Neumayer, and D. G. GFZ (2013), GFZ GRACE Level-2 Processing Standards Document for Level-2 Product Release 0005, Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany.
Dee, D. P., et al. (2011), The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, Q. J. R. Meteorol. Soc., 137(656), 553–597, doi:10.1002/qj.828.
Dill, R. (2008), Hydrological model LSDM for operational Earth rotation and gravity eld variations, Sci. Tech. Rep. 08/09, Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany, doi:11.2312/GFZ.b103-08095.
Ditmar, P., J. T. D. Encarnação, and H. H. Farahani (2012), Understanding data noise in gravity field recovery on the basis of inter-satellite ranging measurements acquired by the satellite gravimetry mission GRACE, J. Geod., 86(6), 441–465, doi:10.1007/s00190-011-0531-6. Dobslaw, H., I. Bergmann-Wolf, R. Dill, E. Forootan, V. Klemann, J. Kusche, and I. Sasgen (2015), The updated ESA Earth System Model for
future gravity mission simulation studies, J. Geod., 89, 505–513, doi:10.1007/s00190-014-0787-8.
Dobslaw, H., I. Bergmann-Wolf, E. Forootan, C. Dahle, T. Mayer-Gürr, J. Kusche, and F. Flechtner (2016), Modeling of present-day atmosphere and ocean non-tidal de-aliasing errors for future gravity mission simulations, J. Geod., 5, 423–436, doi:10.1007/s00190-015-0884-3. Dong, D., P. Fang, Y. Bock, M. K. Cheng, and S. Miyazaki (2002), Anatomy of apparent seasonal variations from GPS-derived site position time
series, J. Geophys. Res., 107(B4), 2075, doi:10.1029/2001JB000573.
Ettema, J., M. R. van den Broeke, E. van Meijgaard, W. J. van de Berg, J. L. Bamber, J. E. Box, and R. C. Bales (2009), Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling, Geophys. Res. Lett., 36, L12501, doi:10.1029/2009GL038110. Farrell, W. E. (1972), Deformation of the Earth by surface loads, Rev. Geophys., 10(3), 761–797, doi:10.1029/RG010i003p00761. Farrell, W. E., and J. A. Clark (1976), On postglacial sea level, Geophys. J. R. Astron. Soc., 46(3), 647–667,
doi:10.1111/j.1365-246X.1976.tb01252.x.
Flechtner, F., and H. Dobslaw (2013), AOD1b Product Description Document for Product Release. 05, GFZ German Res. Cent. for Geosci., Potsdam, Germany.
Fukumori, I. (2002), A partitioned Kalman filter and smoother, Mon. Weather Rev., 130(5), 1370–1383, doi:10.1175/1520-0493(2002)130<1370:APKFAS>2.0.CO;2.
Gordeev, R. G., B. A. Kagan, and E. V. Polyakov (1977), The effects of loading and self-attraction on global ocean tides: The model and the results of a numerical experiment, J. Phys. Oceanogr., 7(2), 161–170, doi:10.1175/1520-0485(1977)007<0161:TEOLAS>2.0.CO;2. Griffiths, J., and J. R. Ray (2013), Sub-daily alias and draconitic errors in the IGS orbits, GPS Solutions, 17(3), 413–422,
doi:10.1007/s10291-012-0289-1.
Gruber, T., J. L. Bamber, M. F. P. Bierkens, H. Dobslaw, M. Murböck, M. Thomas, L. P. H. van Beek, T. van Dam, L. L. A. Vermeersen, and P. N. A. M. Visser (2011), Simulation of the time-variable gravity field by means of coupled geophysical models, Earth Syst. Sci. Data, 3(1), 19–35, doi:10.5194/essd-3-19-2011.
Kim, S.-B., T. Lee, and I. Fukumori (2007), Mechanisms controlling the interannual variation of mixed layer temperature averaged over the Niño-3 region, J. Clim., 20(15), 3822–3843, doi:10.1175/JCLI4206.1.
Klees, R., E. A. Revtova, B. C. Gunter, P. Ditmar, E. Oudman, H. C. Winsemius, and H. H. G. Savenije (2008), The design of an optimal filter for monthly GRACE gravity models, Geophys. J. Int., 175(2), 417–432, doi:10.1111/j.1365-246X.2008.03922.x.
Kusche, J., and E. J. O. Schrama (2005), Surface mass redistribution inversion from global GPS deformation and Gravity Recovery and Climate Experiment (GRACE) gravity data, J. Geophys. Res., 110, B09409, doi:10.1029/2004JB003556.
Acknowledgments
We thank CSR for making the full noise variance-covariance matrices publicly available. Y.S. would like to thank his sponsor, the Chinese Scholarship Council, and two guaran-tors, Xiaotao Chang and Jinyun Guo. Y.S. is also partly supported by the 973 program under grant 2013CB733302-2. R.R. acknowledges support by the Netherlands Organization for Scientific Research (NWO), through VIDI grant 864.12.012. All data used in this study are publicly available. SLR [Cheng et al., 2013b] is downloaded from from ftp://ftp.csr.utexas.edu/pub/slr/ geocenter/GCN_RL05.txt. SLR [Lemoine et al., 2013] is extracted from GRGS RL03 GRACE/SLR monthly GRACE gravity field solutions (http://gravitegrace.get.obs-mip.fr/ grgs.obs-mip.fr/). SLR [So´snica et al., 2014] is extracted from AIUB SLR monthly gravity field (6 by 6∘) solutions (ftp://ftp.unibe.ch/aiub/ GRAVITY/SLR/). Links for other data are provided in the main text.
