Improved Estimates of Geocenter Variability from Time-Variable Gravity and Ocean Model Outputs

Abstract

Geocenter variations relate the motion of the Earth’s center of mass with respect to its center of figure, and represent global-scale redistributions of the Earth’s mass. We investigate different techniques for estimating of geocenter motion from combinations of time-variable gravity measurements from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On missions, and bottom pressure outputs from ocean models. Here, we provide self-consistent estimates of geocenter variability incorporating the effects of self-attraction and loading, and investigate the effect of uncertainties in atmospheric and oceanic variation. The effects of self-attraction and loading from changes in land water storage and ice mass change affect both the seasonality and long-term trend in geocenter position. Omitting the redistribution of sea level affects the average annual amplitudes of the x, y, and z components by 0.2, 0.1, and 0.3 mm, respectively, and affects geocenter trend estimates by 0.02, 0.04 and 0.05 mm/yr for the the x, y, and z components, respectively. Geocenter estimates from the GRACE Follow-On mission are consistent with estimates from the original GRACE mission.

1. Introduction

Variations in the Earth’s geocenter reflect the largest scale variability of mass within the Earth system, and are essential inclusions for the complete recovery of surface mass change from time-variable gravity [1,2]. The Earth’s geocenter is the difference between the Earth’s center of mass and center of figure, which is represented by the degree one spherical harmonic terms [3,4]. Estimates of geocenter position have important applications in the determination of terrestrial reference frame variations, satellite altimeter orbit fluctuations, and mass transport from time-variable gravity [1,5]. Geocenter positions are not stationary. Variations in geocenter have been attributed to changes in terrestrial water storage, glacier and ice sheet mass, atmospheric and oceanic circulation, geodynamic processes and other mass transport processes, Figure 1 [1,6,7,8].

Measurements of time-variable gravity from the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) missions are set in a center of mass (CM) reference frame, in which the total degree one variations are inherently zero [9,10,11]. However, the individual contributions to degree one variations in the CM reference frame, such as from oceanic processes or terrestrial water storage change, are not necessarily zero [4]. Applications set in a center of figure (CF) reference frame, such as the recovery of mass variations of the oceans, hydrosphere and cryosphere, also require the inclusion of degree one terms to be fully accurate [2,4]. The exclusion of degree one terms can have a significant impact on estimates of ocean mass [12], ice sheet mass change [13] and terrestrial hydrology [14] due to far-field signals leaking into each regional estimate.

There are presently two methods that regularly provide estimates of geocenter variability: (1) measurements from satellite laser ranging (SLR) [5,15] and (2) calculations from time-variable gravity and ocean model outputs [2,16]. Global inversions of GRACE data, GPS data, and modeled ocean bottom pressure can also provide long-term or detrended estimates of geocenter variability [17,18]. Trends in SLR-derived solutions can be contaminated by network effects, such as station drift [19,20]. They are not necessarily related to long-term changes in geocenter position. Geocenter estimates from inversions of time-variable gravity and ocean model outputs can provide trend estimates. However, they require information about the oceanic redistribution of mass from changes in terrestrial water storage and ice mass change [2]. Here, we improve geocenter estimates from time-variable gravity and ocean model output combinations by including the effects of self-attraction and loading. We test the overall sensitivity of the estimates to uncertainties in atmospheric pressure, ocean bottom pressure and glacial isostatic adjustment. In the following sections, we discuss (1) the data utilized for this research, (2) how we implement and test our method, (3) the results from our method and (4) the overall implications of the research.

2. Data

2.1. Time-Variable Gravity

We use monthly GRACE/GRACE-FO Release-6 (RL06) gravity solutions provided by the University of Texas Center for Space Research (UTCSR), the German Research Centre for Geosciences (GeoForschungsZentrum, GFZ) and the Jet Propulsion Laboratory (JPL) for the period April 2002–May 2019 [9,10,11,21]. Each Level-2 gravity field solution consists of fully-normalized spherical harmonic coefficients (${\tilde{C}}_{lm}$, ${\tilde{S}}_{lm}$) of degree, l, and order, m. We substitute the ${\tilde{C}}_{20}$ coefficients derived from GRACE/GRACE-FO data with estimates from satellite laser ranging (SLR) from NASA Goddard Space Flight Center due to anomalous variability in the GRACE-derived coefficients [22]. For periods when the satellite pairs had a single fully-operating accelerometer, we replace the ${\tilde{C}}_{30}$ coefficients derived from GRACE/GRACE-FO data with SLR-derived estimates [22]. The monthly GRACE/GRACE-FO RL06 products have non-tidal atmospheric mass variability removed using outputs from the European Centre for Medium-Range Weather Forecasts (ECMWF) and ERA-Interim [23]. Non-tidal oceanic mass variability is removed from the Release-6 GRACE/GRACE-FO data using outputs from the Max Planck Institute ocean model (MPIOM) [24]. We correct for the effects of Glacial Isostatic Adjustment (GIA), which is the viscoelastic response to changes in ice mass since the last glacial maximum, using outputs from the A et al. [25] compressible Earth model using ICE-6G ice history. The impact of different ice histories and mantle rheologies are tested using expected GIA rates from Caron et al. [26]. We account for the elastic deformation of the solid Earth induced by variations in surface mass loading using load Love numbers of gravitational potential, $k_l$, calculated by Wahr et al. [27]. The degree one load Love number of gravitational potential, $k_1$, is updated to reflect the use of a CF reference frame [28]. We smooth the GRACE/GRACE-FO coefficients using a 300 km radius Gaussian averaging function to reduce the impact of random spherical harmonic errors [4,29]. We filter the coefficients to reduce the impact of correlated north/south “striping” errors [30].

As an independent assessment of geocenter variability, we use monthly estimates derived from satellite laser ranging (SLR) [15]. The traditional SLR-derived geocenter solutions reflect the best fit to the center of the SLR network (CN) [31]. However, SLR-derived solutions can be affected by network-effects when individual stations drift relative to each other due to tectonics, surface loading or other processes [19]. Here, we use both the traditional CN-CM solutions [15] and the new UTCSR CF-CM solutions that try to account for these network-effect range biases [5]. We correct both sets of SLR solutions for the gravitational effects of non-tidal atmospheric and oceanic variability using the 6-h Release-6 de-aliasing product provided by GFZ [23].

2.2. Atmospheric Reanalyses and Ocean Models

We estimate the effect of uncertainty in atmospheric pressure by comparing the geocenter outputs from the GRACE/GRACE-FO atmospheric de-aliasing products (GAA) with estimates from ERA-Interim [32], MERRA-2 [33], NCEP-DOE-2 [34] and JRA-55 [35] reanalysis outputs. ERA-Interim is computed by ECMWF and is available starting from 1979 [32]. MERRA-2 is computed by the NASA Global Modeling and Assimilation Office (GMAO) and is available starting from 1980 [33]. NCEP-DOE-2 is computed by the National Centers for Environmental Prediction (NCEP) and is available starting from 1979 [34]. JRA-55 is computed by the Japan Meteorological Agency (JMA) and is available starting from 1958 [35]. Atmospheric pressure anomalies for each reanalysis are calculated relative to the mean over 2001–2002.

We estimate the effect of uncertainty in oceanic circulation by substituting the GRACE/GRACE-FO ocean bottom pressure (OBP) product (GAD) derived from MPIOM with estimates from ECCO-JPL near real-time Kalman-filtered (kf080i) simulations [36,37] and ECCO Version 4 Release 3 (V4r3) simulations [38,39]. Each ocean model incorporates different model physics, data assimilation schemes, and atmospheric forcings. MPIOM is a global ocean model forced with atmospheric outputs from ECMWF medium-range forecasts, and is coupled with a prognostic sea ice model [24]. ECCO-JPL is a regional ocean model configured to resolve tropical ocean circulation (79.5°S–78.5°N latitudinal limits) that is forced with atmospheric outputs from NCEP Reanalysis, and the model does not include sea ice [36]. ECCO V4r3 is a global ocean model that is forced with atmospheric outputs from ERA-Interim, and is coupled with a sea ice model [38,40,41]. The OBP data from ECCO V4r3 were converted from anomalies relative to depth, ${\varphi }_{bot}$, into a time series of absolute OBP in pascals using the model standard density, ${\rho }_0$, and the average gravitational acceleration of the Earth, g. As Boussinesq-type models conserve volume rather than mass, OBP anomalies for each ECCO model were calculated relative to the global average OBP at each time step [42]. Temporal anomalies in OBP are calculated relative to the mean over 2003–2007, following the JPL Tellus OBP product documentation [36].

3. Methods

A change in surface mass density, $\sigma (\theta ,\varphi ,t)$, of a thin layer at the Earth’s surface at time, t, colatitude, $\theta$, and longitude, $\varphi$, can be decomposed into a series of fully-normalized spherical harmonic coefficients, ${\tilde{C}}_{lm}\left(t\right)$ and ${\tilde{S}}_{lm}\left(t\right)$, after allowing for the Earth’s elastic response [4].

$$\left\{\begin{array}{c}{\tilde{C}}_{lm}\left(t\right)\\ {\tilde{S}}_{lm}\left(t\right)\end{array}\right\}=\frac{3}{4\pi a{\rho }_e}\frac{1+k_l}{2l+1}\int \sigma (\theta ,\varphi ,t) P_{lm}\left(\cos \theta \right)\left\{\begin{array}{c}\cos m\varphi \\ \sin m\varphi \end{array}\right\} d\mathsf{\Omega }$$

(1)

where a is average radius of the Earth, ${\rho }_e$ is the average density of the Earth, $k_l$ is the gravitational load Love number of degree l, $P_{lm}$ is the associated Legendre Polynomial of degree l and order m, and $d\mathsf{\Omega }$ is the element of solid angle $\sin \theta d\theta d\varphi$. The change in the surface mass density field, $\sigma (\theta ,\varphi ,t)$, can be calculated through a summation of the fully-normalized spherical harmonics as shown in Wahr et al. [4].

$$\sigma (\theta ,\varphi ,t)=\frac{a{\rho }_e}{3}\sum _{l=0}^{\infty }\sum _{m=0}^l\frac{2l+1}{1+k_l} P_{lm}\left(\cos \theta \right)\left\{{\tilde{C}}_{lm}\left(t\right)\cos m\varphi +{\tilde{S}}_{lm}\left(t\right)\sin m\varphi \right\}$$

(2)

For this work, we use surface mass spherical harmonic coefficients, denoted here as $C_{lm}\left(t\right)$ and $S_{lm}\left(t\right)$, which are calculated from the fully-normalized spherical harmonic coefficients, ${\tilde{C}}_{lm}\left(t\right)$ and ${\tilde{S}}_{lm}\left(t\right)$.

$$\left\{\begin{array}{c}{C}_{lm}\left(t\right)\\ S_{lm}\left(t\right)\end{array}\right\}=\frac{a{\rho }_e}{3}\frac{2l+1}{1+k_l}\left\{\begin{array}{c}{\tilde{C}}_{lm}\left(t\right)\\ {\tilde{S}}_{lm}\left(t\right)\end{array}\right\}$$

(3)

Time-variable gravity measurements from GRACE and GRACE-FO are set in a reference frame of instantaneous center of mass (CM) [9,10,11]. In this reference frame, spherical harmonic coefficients of degree one, $C_{10}$, $C_{11}$, and $S_{11}$, are inherently zero [4]. However, applications set in a center of figure (CF) reference frame require estimates of degree one variability to be accurate [2,4]. Here, we estimate the degree one variations by combining time-variable gravity measurements with ocean model outputs following Swenson et al. [2]. We start by partitioning the surface mass density changes into the individual land and ocean components (Equation (4)). For this, we use an ocean function, $\vartheta (\theta ,\varphi )$, with coastlines buffered by 300 km in order to limit the leakage of mass from the land into the ocean estimate [4].

$$\begin{array}{cc}\hfill \sigma (\theta ,\varphi ,t)& ={\sigma }_{land}(\theta ,\varphi ,t)+{\sigma }_{ocean}(\theta ,\varphi ,t)\hfill \\ \hfill {\sigma }_{ocean}(\theta ,\varphi ,t)& =\vartheta (\theta ,\varphi ) \sigma (\theta ,\varphi ,t)\hfill \end{array}$$

(4)

The oceanic components of the degree one spherical harmonics ($C_{10}^{ocean}$, $C_{11}^{ocean}$, and $S_{11}^{ocean}$) can be calculated from the changes in ocean mass, ${\sigma }_{ocean}(\theta ,\varphi ,t)$ [2,4].

$$\begin{array}{cc}\hfill C_{10}^{ocean}\left(t\right)& =\frac{1}{4\pi }\int P_{10}\left(\cos \theta \right) {\sigma }_{ocean}(\theta ,\varphi ,t) d\mathsf{\Omega }\hfill \\ \hfill C_{11}^{ocean}\left(t\right)& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) {\sigma }_{ocean}(\theta ,\varphi ,t) \cos \varphi d\mathsf{\Omega }\hfill \\ \hfill S_{11}^{ocean}\left(t\right)& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) {\sigma }_{ocean}(\theta ,\varphi ,t) \sin \varphi d\mathsf{\Omega }\hfill \end{array}$$

(5)

Assuming that the changes in oceanic mass can be determined from the global surface mass density field using an ocean function (Equation (4)), the oceanic contributions to degree one can also be estimated from the global set of spherical harmonics [2,4]. Here, we separate the degree one terms in the spherical harmonic summation from the higher degree terms.

$$\begin{array}{cc}\hfill C_{10}^{ocean}\left(t\right)=\frac{C_{10}\left(t\right)}{4\pi }& \int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega } +\hfill \\ \hfill \frac{C_{11}\left(t\right)}{4\pi }& \int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{S_{11}\left(t\right)}{4\pi }& \int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{1}{4\pi }& \int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill C_{11}^{ocean}\left(t\right)=\frac{C_{10}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega } +\hfill \\ \hfill \frac{C_{11}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{S_{11}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{1}{4\pi }& \int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill S_{11}^{ocean}\left(t\right)=\frac{C_{10}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega } +\hfill \\ \hfill \frac{C_{11}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{S_{11}\left(t\right)}{4\pi }& \int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega } +\hfill \\ \hfill \frac{1}{4\pi }& \int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \end{array}$$

(8)

If the oceanic contributions to degree one variability ($C_{10}^{ocean}$, $C_{11}^{ocean}$, and $S_{11}^{ocean}$) can be estimated from an ocean model, then the unknown complete degree one terms ($C_{10}$, $C_{11}$, and $S_{11}$) in Equations (6)–(8) can be calculated from the residual between the oceanic degree one terms and the measured mass change over the ocean calculated using all other degrees of the global spherical harmonics [2].

$$\left[\begin{array}{c}{C}_{10}\left(t\right)\\ C_{11}\left(t\right)\\ S_{11}\left(t\right)\end{array}\right]={\left[\begin{array}{ccc}{I}_{10C}^{10C}& I_{11C}^{10C}& I_{11S}^{10C}\\ I_{10C}^{11C}& I_{11C}^{11C}& I_{11S}^{11C}\\ I_{10C}^{11S}& I_{11C}^{11S}& I_{11S}^{11S}\end{array}\right]}^{-1}\left[\begin{array}{c}{C}_{10}^{ocean}\left(t\right)-G_{10C}\left(t\right)\\ C_{11}^{ocean}\left(t\right)-G_{11C}\left(t\right)\\ S_{11}^{ocean}\left(t\right)-G_{11S}\left(t\right)\end{array}\right]$$

(9)

The I-matrix in Equation (9) is comprised of the degree one terms in the spherical harmonic summations from Equations (6)–(8) [2,4]:

$$\begin{array}{cc}\hfill I_{10C}^{10C}& =\frac{1}{4\pi }\int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega }\hfill \\ \hfill I_{11C}^{10C}& =\frac{1}{4\pi }\int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega }\hfill \\ \hfill I_{11S}^{10C}& =\frac{1}{4\pi }\int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega }\hfill \\ \hfill I_{10C}^{11C}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega }\hfill \\ \hfill I_{11C}^{11C}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega }\hfill \\ \hfill I_{11S}^{11C}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \cos \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega }\hfill \\ \hfill I_{10C}^{11S}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{10}\left(\cos \theta \right) d\mathsf{\Omega }\hfill \\ \hfill I_{11C}^{11S}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \cos \varphi d\mathsf{\Omega }\hfill \\ \hfill I_{11S}^{11S}& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right) \sin \varphi \vartheta (\theta ,\varphi ) P_{11}\left(\cos \theta \right) \sin \varphi d\mathsf{\Omega }\hfill \end{array}$$

(10)

The measured mass changes over the ocean, $G_{10C}$, $G_{11C}$, and $G_{11S}$, are estimated from time-variable gravity measurements of degree 2 and greater as listed in Equations (6)–(8) [2].

$$\begin{array}{cc}\hfill G_{10C}\left(t\right)& =\frac{1}{4\pi }\int P_{10}\left(\cos \theta \right) \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \\ \hfill G_{11C}\left(t\right)& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right)\cos \varphi \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \\ \hfill G_{11S}\left(t\right)& =\frac{1}{4\pi }\int P_{11}\left(\cos \theta \right)\sin \varphi \vartheta (\theta ,\varphi ) \sum _{l=2}^{\infty }\sum _{m=0}^l P_{lm}\left(\cos \theta \right)\left\{C_{lm}\left(t\right)\cos m\varphi +S_{lm}\left(t\right)\sin m\varphi \right\} d\mathsf{\Omega }\hfill \end{array}$$

(11)

3.1. Eustatic Sea Level from Land Surface Fluxes

The MPIOM model used to correct GRACE/GRACE-FO data for ocean circulation does not include the effects of land-sea exchange and the corresponding sea level response [24]. As the land-sea flux signal would be present in the ocean mass variations, we estimate the monthly land-sea flux using GRACE/GRACE-FO time-variable gravity measurements when calculating the oceanic contributions to degree one ($C_{10}^{ocean}$, $C_{11}^{ocean}$, and $S_{11}^{ocean}$) [12]. We use the same 300 km buffered coastline mask to calculate our land function ($\mathcal{L}=1-\vartheta$). At each time step, we use solutions to the sea level equation to calculate the spatial pattern of sea level variation induced by the change in land mass [43,44,45]. These changes in sea level are often referred to as the effects of self-attraction and loading (SAL) or as the sea level fingerprints (SLF) of the land mass change [46]. As the SLFs can differ significantly from global ocean averages, geocenter estimates that assume a uniform redistribution of terrestrial water fluxes can be negatively impacted [16]. Here, we use a pseudo-spectral approach for solving the sea level equation in which we assume the Earth deforms elastically from the modern-day mass change without the manifestation of viscoelastic effects [47,48]. When calculating each SLF, we use a static ocean function with the same 300 km buffered coastlines to verify that mass is being conserved.

3.2. Iterated Solutions

In order to estimate the full component of land-sea mass transport, an initial estimate of the geocenter variability is needed to calculate the total land mass [2]. Initially, we use annual and semi-annual geocenter estimates from Chen et al. [7], which are calculated from land-surface model estimates of soil moisture and snow water equivalent. We then calculate the eustatic sea level change induced from the land-sea flux and use it to generate a full geocenter estimate with Equation (9). The initial geocenter estimate from Chen et al. [7] is then replaced with the newly calculated estimate, and the sequence is repeated until the difference between successive iterations falls below a threshold value (see flowchart in Figure 2).

3.3. Spherical Harmonics of Atmospheric and Oceanic Variability

Changes in atmospheric pressure, p, and ocean bottom pressure, $p_{bot}$, at colatitude, $\theta$, and longitude, $\varphi$, will impact the Earth’s gravitational field. These changes can be decomposed into a series of fully-normalized spherical harmonic coefficients representing the induced gravitational change after allowing for the Earth’s elastic response [49,50]:

$$\left\{\begin{array}{c}{\tilde{C}}_{lm}\left(t\right)\\ {\tilde{S}}_{lm}\left(t\right)\end{array}\right\}=\frac{3}{4\pi {\rho }_e}\frac{1+k_l}{(2l+1)}\int {\xi }_l(\theta ,\varphi ,t) P_{lm}\left(\cos \theta \right)\left\{\begin{array}{c}\cos m\varphi \\ \sin m\varphi \end{array}\right\} d\mathsf{\Omega }$$

(12)

The vertical integral, ${\xi }_l(\theta ,\varphi ,t)$, in Equation (12) is determined based on assumptions of the Earth’s geometry and the vertical structure of the atmosphere or ocean [49].

Here, spherical harmonics from NCEP-DOE-2 and JRA-55 reanalyses are calculated assuming a thin layer two-dimensional atmosphere with a realistic Earth geometry incorporating the model orography, $h(\theta ,\varphi )$, and estimates of geoid height, $N(\theta ,\varphi )$ (Equations (12) and (13)), while harmonics from ERA-Interim and MERRA-2 reanalyses are calculated assuming a three-dimensional atmospheric geometry integrating over the model layers (Equations (12) and (14)) [49].

$${\xi }_l(\theta ,\varphi ,t)={\left(\frac{a+h(\theta ,\varphi )+N(\theta ,\varphi )}{a}\right)}^{l+2}\frac{p_0(\theta ,\varphi ,t)}{g(\theta ,\varphi )}$$

(13)

$${\xi }_l(\theta ,\varphi ,t)=-{\int }_{p_0}^0{\left(\frac{a+z(\theta ,\varphi )+N(\theta ,\varphi )}{a}\right)}^{l+2}\frac{dp}{g(\theta ,\varphi ,z)}$$

(14)

Similarly, spherical harmonics from ECCO kf080i and ECCO V4r3 ocean bottom pressure are calculated assuming a thin layer two-dimensional ocean with a realistic Earth geometry incorporating estimates of geoid height and ocean bathymetry, $d(\theta ,\varphi )$ (Equations (12) and (15)).

$${\xi }_l(\theta ,\varphi ,t)={\left(\frac{a+N(\theta ,\varphi )-d(\theta ,\varphi )}{a}\right)}^{l+2}\frac{p_{bot}(\theta ,\varphi ,t)}{g(\theta ,\varphi )}$$

(15)

3.4. Time Series Analysis

We calculate the average geocenter change by simultaneously fitting a least-squares model with constant, linear and quadratic terms with annual, semiannual and 161-day oscillating terms to the estimate geocenter time series [51]. The 161-day oscillating terms account for aliasing of the S2 tidal constituent in the monthly GRACE/GRACE-FO time-variable gravity fields [52].

4. Results

4.1. Simulated Geocenter Estimates

We test the efficacy of our methodology for estimating geocenter variations by running experiments using fields of simulated global mass variability [2]. Sets of coefficients are constructed using estimates of terrestrial water storage (TWS) calculated from the Global Land Data Assimilation Systems (GLDAS) NOAH land surface model [53], estimates of surface mass balance (SMB) change for glaciers and ice caps calculated from the Regional Atmospheric and Climate Model (RACMO2.3) [54,55] and estimates of Greenland and Antarctic ice sheet mass balance from the mass budget method (MBM) [56]. Monthly GLDAS TWS estimates were calculated by combining the snow water equivalent, canopy water storage and soil moisture variables for non-glaciated regions [57]. We include the mass changes from glaciers and ice sheets to incorporate more processes that affect inter-annual geocenter variability [8]. The SMB of a glacier represents the sum of mass accumulation from snow and rain minus the surface ablation from meltwater runoff, sublimation, and snow drift erosion [58,59]. Cumulative anomalies in SMB were calculated for glaciers and ice caps in reference to a 1961–1990 baseline [60]. MBM estimates of ice sheet mass balance were calculated combining SMB outputs from RACMO2 with estimates of total ice discharge [56,61]. The sea level fingerprints of the synthetic data were calculated as an estimate of oceanic variability [48]. The ocean function used when calculating the sea level fingerprints in the synthetic is an update of Hall et al. [62] that incorporates Antarctic grounded ice delineations [63] and more accurate Greenland coastlines [64].

We test three different scenarios: (1) a uniform ocean redistribution of the land mass change with a static seasonal geocenter estimate similar to Swenson et al. [2], (2) a uniform ocean distribution of the land mass change that iteratively solves for the geocenter and (3) an oceanic redistribution taking into account self-attraction and loading effects that iteratively solves for the geocenter. Each model run uses spherical harmonic coefficients of degree two and greater that are calculated from the synthetic data. The harmonics are truncated to degree and order 60 and processed using the same 300 km Gaussian averaging function and decorrelation filter used with the GRACE/GRACE-FO data [4,29,30]. In order to calculate the full component of the land mass change, we include seasonal estimates of degree one variation from Chen et al. [7]. In the two iteration scenarios, we replace the degree one coefficients for each run with the calculated results of the previous run. The process is repeated until the difference between estimates on successive iterations falls below a threshold value.

We compare the results of each scenario against the exact time series of degree one variations calculated from the synthetic dataset (Figure 3). The RMS differences between the scenarios and the original time series are 0.32 mm, 0.21 mm, and 0.11 mm for the static, iterated, and iterated SLF scenarios, respectively. This represents differences of 21%, 14%, and 7%, respectively, for each of the three scenarios compared to the standard deviation of the original time series. Calculating the land mass change using a static seasonal estimate of geocenter variability underestimates both the short-term and long-term variability in geocenter motion. Self-attraction and loading effects limit the recovery of geocenter variations in the two scenarios that assume a uniform redistribution of land-sea fluxes. The recovery of the seasonal variability of each coefficient is significantly improved when self-attraction and loading effects are incorporated (Figure 4).

4.2. Recovered Geocenter Estimates

We estimate sets of degree one coefficients that are corrected for the effects of non-tidal atmospheric and oceanic variation [2]. These coefficients are directly applicable to time-variable gravity applications for estimating surface mass change [4]. We calculate geocenter variability estimates for each of the processing centers (CSR, GFZ, JPL) of the GRACE/GRACE-FO Science Data System (SDS) using the iterated self-consistent geocenter method that incorporates self-attraction and loading effects (Figure 5). Geocenter results from all three processing centers largely agree in terms of the amplitude and phase of the seasonal signal, along with the long-term trend in each coefficient. Self-attraction and loading effects impact the trend of the estimated geocenter solution for all three processing centers (Figure 6). Using a static seasonal geocenter to estimate the total land-sea flux results in weaker trends for all three processing centers. The trend and annual amplitudes of each coefficient calculated using time-variable gravity fields from each processing center are listed in

Table 1. Geocenter motion annual amplitudes, annual phase, and trends for 2002–2017 derived using satellite laser ranging (SLR) and time-variable gravity fields from the Center for Space Research (CSR), the German Research Centre for Geosciences (GFZ) and the Jet Propulsion Laboratory (JPL) corrected for the effects of non-tidal atmospheric and oceanic variation. Errors denote the 95% confidence level.

	x	y	z
Annual Amplitude [mm]
CSR	1.34 ± 0.11	1.54 ± 0.12	2.25 ± 0.16
GFZ	1.38 ± 0.14	1.56 ± 0.13	2.30 ± 0.16
JPL	1.31 ± 0.11	1.52 ± 0.12	2.20 ± 0.17
SLR CN-CM	1.93 ± 0.38	1.17 ± 0.38	4.25 ± 0.57
SLR CF-CM	1.29 ± 0.29	1.48 ± 0.23	2.97 ± 0.46
Annual Phase [day]
CSR	356.5 ± 5.0	151.9 ± 4.6	9.6 ± 4.2
GFZ	352.7 ± 5.7	150.7 ± 4.9	4.5 ± 4.1
JPL	355.0 ± 4.8	151.4 ± 4.7	7.9 ± 4.4
SLR CN-CM	5.5 ± 11.7	194.7 ± 18.9	52.4 ± 7.7
SLR CF-CM	347.9 ± 13.3	169.6 ± 9.0	46.3 ± 9.1
Trend [mm/yr]
CSR	−0.15 ± 0.02	0.10 ± 0.02	−0.62 ± 0.03
GFZ	−0.19 ± 0.03	0.21 ± 0.03	−0.66 ± 0.03
JPL	−0.15 ± 0.02	0.11 ± 0.03	−0.63 ± 0.03

. Annual amplitudes for solutions derived from satellite laser ranging (SLR) after correcting for non-tidal atmospheric and oceanic variability are listed for comparison [5,15]. The listed uncertainties for each coefficient are 95% confidence intervals, but do not take into account spherical harmonic errors or uncertainties in the geophysical corrections.

4.3. Uncertainty Estimates

The total uncertainty in the estimated geocenter time series represents a combination of GRACE/GRACE-FO measurement error, signal leakage between ocean and land, and uncertainty in the geophysical corrections. We assess the impact of these uncertainties, which we assume are uncorrelated, on the trend, annual amplitude and annual phase to the 95% confidence level. Monthly errors in the GRACE.GRACE-FO Level-2 spherical harmonics are estimated by first smoothing each coefficient with a 13-month Loess-type algorithm [13] and then calculating the residuals between the original and smoothed coefficients [65]. Spherical harmonic uncertainties contribute 0.11, 0.09 and 0.19 mm to the x, y, and z components for coefficients derived from CSR time-variable gravity fields.

We estimate the contribution from uncertainty in ocean circulation using ocean bottom pressure outputs from ECCO-JPL near real-time Kalman-filtered simulations and ECCO V4r3 simulations [36,38]. The time series of recovered degree one coefficients derived from CSR time-variable gravity fields using ocean bottom pressure outputs from ECCO-JPL kf080i [36], ECCO V4r3 [38] and MPIOM [23] are shown in Figure 7. The trend and annual amplitudes of each coefficient derived using the different ocean bottom pressure estimates are listed in

Table 2. Geocenter motion annual amplitudes, annual phase, and trends for 2002–2015 derived using time-variable gravity fields from the Center for Space Research (CSR) using ocean bottom pressure outputs from ECCO-JPL real-time Kalman-filtered simulations (kf080i) [36] and ECCO Version 4 Release 3 simulations (V4r3) [38]. Errors denote the 95% confidence level.

	x	y	z
Annual Amplitude [mm]
ECCO-JPL kf080i	1.46 ± 0.20	1.28 ± 0.17	1.80 ± 0.31
ECCO V4r3	1.63 ± 0.18	1.21 ± 0.16	2.31 ± 0.27
MPIOM	1.34 ± 0.11	1.55 ± 0.13	2.24 ± 0.16
Annual Phase [day]
ECCO-JPL kf080i	307.3 ± 7.9	165.5 ± 7.9	327.8 ± 9.9
ECCO V4r3	323.0 ± 6.5	150.8 ± 7.7	326.2 ± 6.9
MPIOM	358.3 ± 5.0	150.3 ± 4.8	9.6 ± 4.3
Trend [mm/yr]
ECCO-JPL kf080i	−0.32 ± 0.04	0.16 ± 0.03	−0.48 ± 0.06
ECCO V4r3	−0.10 ± 0.04	0.12 ± 0.04	−0.44 ± 0.06
MPIOM	−0.12 ± 0.02	0.08 ± 0.03	−0.62 ± 0.03

. Solutions from ECCO-JPL agree with ECCO V4r3 and MPIOM derived solutions for the y and z components, but differs in terms of trend in the x-component [2]. Solutions from ECCO V4r3 and MPIOM agree for all coefficients between 2002 and 2014.

We estimate the uncertainty in atmospheric circulation by using an ensemble of reanalysis outputs [32,33,34,35]. Differences between the GRACE/GRACE-FO atmospheric de-aliasing product (GAA) and outputs from ERA-Interim [32], MERRA-2 [33], NCEP-DOE-2 [34], and JRA-55 [35] reanalyses will affect the estimates of geocenter motion. The average monthly uncertainty in geocenter variation due to uncertainties in atmospheric circulation is 0.07, 0.10, and 0.24 mm for the x, y, and z components, respectively.

Glacial Isostatic Adjustment (GIA) affects both ocean mass and land-sea flux calculations used to reconstruct the geocenter variation. The contribution of GIA uncertainty is calculated by comparing our best case estimate against the expected GIA rate from Caron et al. [26]. Using coefficients from Caron et al. [26] affects trend estimates by 0.04, 0.003, and 0.20 mm/yr for the x, y, and z components, respectively. This represents trend differences in the x, y, and z components of 21%, 3%, and 28%, respectively, compared with the values derived using GIA outputs from A et al. [25] with ICE-6G ice history. Uncertainty in GIA is a limiting factor for determining rates of geocenter change from reconstructions with time-variable gravity and ocean models.

5. Discussion

There is better agreement between the estimates derived from GRACE/GRACE-FO and the new CF-CM SLR-derived solutions compared with SLR-derived CN-CM estimates Figure 8, Table 1 [5]. Self-attraction and loading effects improve the correspondence with SLR-derived solutions, particularly in terms of average annual amplitude (Figure 9). However, there are still differences in the seasonal amplitudes of the z-component of geocenter motion ($C_{10}$). Using a different ocean model to derive the oceanic component of geocenter variability cannot fully explain the difference between the solutions (Figure 8). Uncertainties in the GRACE spherical harmonics and uncertainties in atmospheric variation can partially explain the differences in the solutions. In addition, uncertainties in the ability to correct for network-effects in the SLR-derived CF-CM solutions can also help explain some discrepancies with the solutions derived from GRACE/GRACE-FO.

There is strong agreement between the three processing centers for each of the degree one coefficients during both the GRACE and GRACE-FO periods. However, between the end of 2015 and the end of the GRACE mission, there is disagreement between centers in estimates of the x and y components of geocenter motion ($C_{11}$ and $S_{11}$). Operational procedures enacted to maintain the battery life of the GRACE satellites during the latter stages of the mission affected the quality of the time-variable gravity fields. The accelerometer onboard GRACE-B was turned off in September 2016 to reduce the battery load and maintain the operation of the microwave ranging instrument. Independent methods have been developed by the GRACE processing centers to spatiotemporally transplant the accelerometer data retrieved from GRACE-A for GRACE-B [66]. The increased uncertainty in the gravitational fields likely affects our ability to estimate the geocenter variability, particularly for the months with a single operating accelerometer.

Swenson et al. [2] provides a method for deriving geocenter motion from time-variable gravity measurements and estimates of ocean bottom pressure. The full component of oceanic mass is estimated in Swenson et al. [2] by calculating the land-sea fluxes using time-variable gravity measurements. We expand upon this method by using an iterated, self-consistent geocenter when calculating the full component of the land mass change, and by including self-attraction and loading effects when redistributing the mass over the ocean. The necessity of these inclusions only became evident with a longer record of time-variable gravity measurements. Using a static, seasonal geocenter omits the contribution of inter-annual fluctuations in degree one to the total land mass change. Self-attraction and loading effects impact both the seasonal amplitudes and the long-term trend of the recovered geocenter estimate (Figure 6 and Figure 9). This is due to the spatial patterns of sea level fingerprints that can deviate strongly from the uniform average, particularly over the long-term from changes in glaciers and ice sheets [48,67].

Sun et al. [16] expanded on Swenson et al. [2] by simultaneously estimating oblateness, $C_{20}$, variations along with geocenter variations. They test the sensitivity of geocenter estimation methods to geophysical processes, such as glacial isostatic adjustment and ocean redistributions due to self-attraction and loading. In Sun et al. [16] the GRACE spherical harmonics are truncated to degree and order 45 and not processed with a decorrelation filter. Here, we use the full expansion of spherical harmonics for each processing center and filter for correlated errors in the GRACE/GRACE-FO harmonics [30]. The test the efficacy of the two methods using our synthetic reconstruction of global mass change. While statistically similar between the two techniques, we find that expanding to higher degree and orders and filtering produces more consistent estimates compared with synthetic estimate. The Sun et al. [16] method is used when calculating the JPL Tellus geocenter product. The differences between the geocenter estimates computed here and the estimates provided by JPL Tellus are largely during the final months of the GRACE mission and initial months of the GRACE-FO mission. Estimates of Antarctic ice sheet mass balance from time-variable gravity are sensitive to uncertainties in geocenter variation [68]. We find that differences between the two geocenter estimates affect Antarctic ice sheet mass balance estimates by 8–10 Gt/yr between 2002 and 2019.

Here, we use a buffered 300 km ocean function to calculate the sea level fingerprints in order to conserve mass and to reduce the leakage between land and ocean [4]. A full-resolution ocean function with realistic coastlines could be used if sets of scaling factors could be derived for both the land and ocean mass change. Typically, the use of scaling factors is applicable only for deriving seasonal fluctuations in land mass [69]. A more-complete scaling factor, such as from Hsu and Velicogna [67] for the land mass change, could possibly improve estimates of geocenter variation. The effects of self-attraction and loading would likely become more evident with a full-resolution ocean function as the effects can be more pronounced in coastal areas [48].

6. Conclusions

Geocenter variations are an important representation of global mass change, with applications in satellite gravimetry and in satellite orbit determination. Here, we investigate the effects of different calculations of the eustatic sea level caused by changes in land mass for estimating geocenter variations from combinations of ocean model outputs and time-variable gravity measurements from the GRACE and GRACE Follow-On missions. We find that using an iterated, self-consistent geocenter incorporating self-attraction and loading effects provides the best estimate of geocenter variation. Uncertainties from glacial isostatic adjustment, ocean circulation and atmospheric circulation limit the determination of geocenter position from this technique. The annual amplitudes of the z-component of geocenter variation differs between estimates from this technique and estimates derived from satellite laser ranging (SLR). The effects of self-attraction and loading improve the correspondence with SLR-derived coefficients but there are still discrepancies worth further investigation. Estimates of geocenter variations using data from the GRACE Follow-On mission are consistent with estimates from the GRACE mission, which enables the extension of the geocenter time series going forward.