Reducing filter effects in GRACE-derived polar motion excitations

Polar motion is caused by mass redistribution and motion within the Earth system. The GRACE (Gravity Recovery and Climate Experiment) satellite mission observed variations of the Earth’s gravity field which are caused by mass redistribution. Therefore GRACE time variable gravity field models are a valuable source to estimate individual geophysical mass-related excitations of polar motion. Since GRACE observations contain erroneous meridional stripes, filtering is essential to retrieve meaningful information about mass redistribution within the Earth system. However filtering reduces not only the noise but also smoothes the signal and induces leakage of neighboring subsystems into each other. We present a novel approach to reduce these filter effects in GRACE-derived equivalent water heights and polar motion excitation functions which is based on once- and twice-filtered gravity field solutions. The advantages of this method are that it is independent from geophysical model information, works on global grid point scale and can therefore be used for mass variation estimations of several subsystems of the Earth. A closed-loop simulation reveals that due to application of the new filter effect reduction approach the uncertainties in GRACE-derived polar motion excitations can be decreased from 12–48% to 5–29%, especially for the oceanic excitations. Comparisons of real GRACE data with model-based oceanic excitations show that the agreement can be improved by up to 15 percentage points.


Introduction
Redistribution and motion of masses within the Earth system cause polar motion and length-of-day (LOD) variations. Initially, individual geophysical excitation mechanisms of Earth rotation were estimated from geophysical fluid models for atmosphere, oceans and continental hydrosphere (e.g., Gross et al. 2004;Seitz and Schmidt 2005;Chen and Wilson 2005;Nastula et al. 2007;Brzezinski et al. 2009;Dobslaw et al. 2010). The modeling of geodynamic processes within and between the individual subsystems of the Earth is very complex and suffers from uncertainties in the process descriptions, parametrization and forcing (Güntner 2008). Thus, large discrepancies exist between different model solutions for the individual contributions to Earth rotation, especially for the continental hydrosphere.
Between 2002 and 2017, the gravity satellite mission GRACE (Gravity Recovery And Climate Experiment) observed the time variable gravity field of the Earth arising from mass redistribution within the Earth system. Hence, geophysical mass-related excitations of Earth rotation can be determined from time variable gravity field models. Since LOD variations are related to the degree-2 spherical harmonic coefficient C 20 , which still suffers from ocean tide model errors within the GRACE data processing (Chen and Wilson 2008), most GRACE Earth rotation studies are restricted to polar motion investigations like this study. Initially, GRACE time variable gravity field models were used to estimate the integral mass-related excitations of polar motion (e.g., Cheng and Tapley 2004;Nastula et al. 2007;Chen and Wilson 2008;Gross 2009;Seoane et al. 2009;Seoane and Gambis 2012), later they were also used to estimate the individual contributions Göttl et al. Earth, Planets and Space (2019) 71:117 of the continental hydrosphere, oceans and cryosphere separately (e.g., Jin et al. 2010Jin et al. , 2012Seoane et al. 2011;Chen et al. 2012Chen et al. , 2013Chen et al. , 2017Göttl et al. 2012Göttl et al. , 2015Göttl et al. , 2018Göttl 2013;Adhikari and Ivins 2016;Meyrath and van Dam 2016).
Since GRACE level-2 products contain meridional error stripes due to the GRACE sensor characteristics and mission geometry, filtering is essential to retrieve meaningful information about mass redistribution within the Earth system (Wahr et al. 1998). However, filtering reduces not only the noise but also smooths the signal and induces that neighboring subsystems leak into each other. Therefore, the separation into the contributions of the individual subsystems of the Earth still remains a challenge. There exists a large number of papers concerning the reduction of leakage effects in GRACE-derived mass variations of the continental hydrosphere (e.g., Klees et al. 2007;Longuevergne et al. 2010;Landerer and Swenson 2012;Vishwakarma et al. 2016Vishwakarma et al. , 2017Khaki et al. 2018), oceans (e.g., Chambers 2006Chambers and Bonin 2012;Peralta-Ferris et al. 2014) and cryosphere (e.g., Baur et al. 2009;King et al. 2012;Velicogna and Wahr 2013;Chen et al. 2015;Mu et al. 2017). Many of these techniques depend on geophysical model information such as the forward modelling approaches (additive, multiplicative and scaling) which are designed for specific subsystems of the Earth or regions. Chen et al. (2012) demonstrated that leakage corrections are important to quantify adequate hydrological and cryospheric polar motion excitations from GRACE data. We present a novel filter effect reduction approach which does not depend on geophysical model information and which can be applied to several subsystems of the Earth. The paper is organized as follows: the next section provides information on time variable gravity field solutions derived from real and simulated GRACE observations and an explanation of the processing steps to determine the mass-related part of the equatorial effective angular momentum functions which describe polar motion excitations of the continental hydrosphere, oceans, Antarctica and Greenland. We introduce a novel approach to reduce the filter effects (attenuation and leakage) in GRACE-derived polar motion excitation functions which is based on once and twice filtered GRACE gravity field solutions. To validate this new method, we perform a closed-loop simulation and we apply the method to real GRACE data to study the impact on GRACE-derived polar motion excitation functions. Comparisons with ocean model results are performed. Finally, the last section provides the conclusions.

GRACE gravity field models
Time variable gravity field models, derived from GRACE observations, can be used to estimate both the integral mass-related excitations of polar motion and the individual contributions of the continental hydrosphere, oceans and cryosphere. In this study, we use the new Release-06 (RL06) GRACE gravity field solutions from the German Research Centre for Geosciences (GFZ), Potsdam (Dahle et al. 2019), Center for Space Research (CSR), Austin (Bettadpur 2018), and NASA's Jet Propulsion Laboratory (JPL), Pasadena (Yuan 2018) as well as the new ITSG-Grace2018 gravity field solutions provided by the Institute of Theoretical Geodesy and Satellite Geodesy (ITSG) of Graz University of Technology (Mayer-Gürr et al. 2018) which all incorporate the new RL06 background models. Due to the latest release update from RL05 to RL06 the consistency of the GRACE gravity field solutions could be significantly improved as well as the signal-to-noise ratio (Göttl et al. 2018). They are provided as Level-2 products GSM, GAC and GAD. While the GSM products represent the Earth gravity field excluding tidal effects and short-term atmospheric and oceanic variations, the GAC products represent the non-tidal short-term effects of the atmosphere and oceans which are reduced within the GRACE gravity field processing to minimize aliasing effects. The GAD products are provided especially for oceanographers, they represent ocean bottom pressure variations. We use the sum of the GSM and GAD products to identify individual excitations of polar motion. The degree-1 Stokes coefficients of the GSM products are equal to zero by definition, because the GRACE data processing is performed in the Earth's center-of-mass (CM) frame. But redistribution of masses at the Earth's surface and within its interior are referred to a coordinate system attached to the Earth's crust which moves relatively to the CM. Therefore, the degree-1 Stokes coefficients are replaced by solutions from Swenson et al. (2008) derived from GRACE data and ocean model outputs. Furthermore, as recommended, we replace the inaccurate C 20 coefficient of the GSM products by an improved external satellite laser ranging (SLR) solution from Cheng et al. (2013) based on GRACE RL06 background models. To separate mass variations in the continental hydrosphere, oceans and cryosphere, the effect of glacial isostatic adjustment (GIA) must be removed from the GRACE observations. We use the global GIA model ICE-6G_D (VM5a) described by Peltier (2012) and Peltier (2015). To retrieve meaningful information about mass redistribution within the Earth system the noisy GRACE gravity field solutions are filtered by the anisotropic decorrelation filter DDK (Kusche 2007) to reduce the meridional error stripes. The DDK filter is based on a regularization of the normal equation system applying the signal covariance matrix and the error covariance matrix of one specific month. Thus the DDK filter is static instead of time varying and its kernel in the space domain is different for each location. Stripes and spurious patterns can be reduced more efficiently than by the isotropic Gaussian filter. Nevertheless the DDK filter also reduces not only the noise but also smooths the signal and induces that neighboring subsystems leak into each other. In this study we developed a novel approach to reduce these filter effects. The method is described in the next section. By applying suitable masks and the global spherical harmonic analysis hydrological, oceanic and cryospheric potential coefficients are obtained that describe mass variations of the individual subsystems of the Earth in the spectral domain. The degree-2 spherical harmonic potential coefficients C 21 and S 21 are directly related to the mass-related part of the equatorial effective angular momentum functions which describe polar motion excitations (Barnes et al. 1983). We use the conversion formulas given by Göttl (2013) to determine the polar motion excitation functions of the continental hydrosphere χ H , oceans χ O , Antarctica χ A and Greenland χ G .

Simulation data
To validate the new filter effects reduction approach we use 4 years of the full-scale closed-loop simulation data described in Flechtner et al. (2016). The observations there are based on a realistic orbit scenario and error assumptions for instruments and background models. In this study we use the monthly GRACE-like gravity field solutions in terms of spherical harmonic coefficients up to degree and order 60. Figure 1 shows the degree amplitudes in terms of geoid heights for the true input signals and the estimated GRACE GFZ RL06 coefficients. The simulated GRACE errors are retrieved by comparing the true and estimated coefficients. It can be clearly seen that the simulated errors of the GRACE gravity field solutions are significantly larger than the formal errors. Therefore in the next section empirical errors of the GRACE gravity field solutions will be estimated.

Reduction of filter effects
This section presents a novel filter effect reduction approach on global grid point basis which has been developed following the methods published by Landerer and Swenson (2012) as well as by Vishwakarma et al. (2016Vishwakarma et al. ( , 2017. Landerer and Swenson (2012) determined gain factors that reduce the differences between the original unfiltered and once filtered terrestrial water storage data from the GLDAS (Global Land-Data Assimilation System, Rodell (2004)) model in a least-squares sense. They show that due to the multiplication of the GRACE filtered gravity field models with model-derived gain factors the filter effects, attenuation and leakage, can be reduced adequately in regions, where the quality of the model is high, but not at locations where specific processes are absent in the model such as melting of ice sheets and glaciers. Therefore this method works only if an accurate geophysical model exists for the area of interest and if the leakage-in effect of the surrounding subsystem is small (e.g., continental hydrosphere). Vishwakarma et al. (2016Vishwakarma et al. ( , 2017) developed a modelindependent approach on catchment scale. They showed that the ratio between true leakage and leakage from once filtered fields is approximately the same than the ratio between leakage from once-and twice-filtered fields in case the catchment is surrounded by catchments with similar signal strength.
Based on the ideas of these two methods we developed a filter effect reduction approach that is independent from geophysical model information and that works for several subsystems of the Earth on global grid point scale. The procedure is presented as a flowchart in Fig. 2.
In a first step once filtered spherical harmonic potential coefficients ( C 1DDK nm , S 1DDK nm ) and equivalent water heights ( ewh 1DDK ) are determined from noisy GRACE gravity field models by applying the DDK filter. Due to the fact that the formal errors of the GRACE gravity field solutions are too optimistic we have to estimate empirical errors in the next step. Therefore the empirical standard deviations of the noisy GRACE gravity field solutions are calculated via Fig. 1 Degree amplitudes in terms of geoid height error for 4 years: simulated GRACE errors (red), formal GRACE GFZ RL06 errors (orange), empirical GRACE errors (green). The grey lines represent the true input signals whereas the blue lines represent the estimated GRACE GFZ RL06 coefficients. The medians of the monthly solutions are shown as bold lines from the unfiltered fully normalized global spherical harmonic potential coefficients C p nm and S p nm provided by specific GRACE processing centers p ∈ 1, ..., P and the average of the spherical harmonic potential coefficients of all GRACE processing centers Here we use the new RL06 gravity field solutions from the GRACE processing centers GFZ, CSR, JPL and ITSG. The time series of the potential coefficients are determined at discrete times t = t k with k = 1, ..., K where K is the total number of months. To obtain the monthly empirical errors of the GRACE gravity field solutions the empirical standard deviations are multiplied by standard normally distributed random values. Figure 1 shows the degree amplitudes of the simulated, formal and empirical GRACE errors. It is clearly visible that the empirical GRACE errors are more realistic than the formal GRACE errors. The empirical errors seem to be an adequate approximation of the real GRACE gravity field errors.
In the next step we add the empirical errors to the oncefiltered potential coefficients and apply the DDK filter a second time to obtain twice-filtered spherical harmonic potential coefficients ( C 2DDK nm , S 2DDK nm ) and equivalent water heights ( ewh 2DDK ). To reduce the discrepancies between the original equivalent water heights (only known in simulation studies) and the once-filtered equivalent water heights we determine global grid point gain factors k org by minimizing the target function The gain factors are mainly larger than 1 over continental grid points especially along the coastlines while they are mainly smaller than 1 over ocean grid points, see Fig. 3. This arises from the large signal strength differences, the larger continental signals leak into the weaker oceanic signals. In reality the original equivalent water heights are not known therefore we determine the global grid point gain factors k from once-and twice-filtered equivalent water heights by minimizing the target function These gain factors k are damped according to the original gain factors k org and they are significantly smoother as shown in Fig. 3. But the basic structure is adequately approximated. By multiplying the once-filtered equivalent water heights with these global grid point gain factors improved equivalent water heights are obtained. Figure 4 shows the relative standard deviations (RSD) of the equivalent water height residuals derived from the simulated true and noisy GRACE gravity field solutions by applying on the on hand only the DDK filter and on the other hand the DDK filter and the global grid point gain factors. One can see that due to the application of the global grid point gain factors the discrepancies with respect to the true data can (7) ewh(θ , , t) = k(θ, ) · (ewh 1DDK (θ, , t)) be significantly reduced, especially in the oceans and in Antarctica. Naturally the gain factors derived from original and once-filtered equivalent water heights (only available in simulation studies) yield better results than the gain factors derived from the once-and twice-filtered equivalent water heights. But the latter appears to be an adequate approximation.

Results and validation
In this section we compare GRACE-derived polar motion excitation functions for the continental hydrosphere, oceans, Antarctica and Greenland with and without consideration of the filter effect reduction approach described in the previous section. The comparisons are performed both for simulation data and real GRACE GFZ RL06 gravity field solutions. Figure 5 shows differences of polar motion excitation functions for the continental hydrosphere, oceans, Antarctica and Greenland derived from simulated noisy GRACE gravity field solutions using only the DDK filter or by applying the DDK filter in combination with the filter effect reduction approach developed within this study with respect to the true polar motion excitation functions. For the oceans χ 1 is the prominent component, whereas χ 2 is dominant for continental hydrosphere, Antarctica and Greenland. The corresponding correlation coefficients and relative standard deviations (RSD) are listed in Table 1. By applying the global grid point gain factors the uncertainties can be reduced from 12-48% to 4-28% for χ i (k org ) and 8-38% for χ i (k) . The gain factors derived from once-and twicefiltered equivalent water heights are a good approximation but they are damped due to the filtering process. To counteract this phenomenon we determine scaling factors f i for each polar motion excitation function χ i (k) from the ratio While the scaling factors for Antarctica and Greenland are larger than 1 ( f A 1 = 1.1 , f A 2 = 1.2 , f G 1 = 1.2 , f G 2 = 1.2 ), they are smaller than 1 for the oceans ( f O 1 = 0.7 , f O 2 = 0.9 ) and equal to 1 for the continental hydrosphere. This can be explained by the fact that the gain factors k are mainly smaller than the original gain factors k org for Antarctica and Greenland whereas for the oceans the opposite is true. Figure 5 shows that due to the  (1) only the DDK filter (green), (2) the DDK filter and the gain factors k org (blue), (3) the DDK filter and the gain factors k (black), (4) the DDK filter, the gain factors k and the scaling factors f i (red) with respect to the effective angular momentum functions derived from the true GRACE gravity field solutions for the continental hydrosphere χ H , oceans χ O , Antarctica χ A and Greenland χ G application of these scaling factors the agreement with the true polar motion excitation functions can be further improved. The uncertainties are reduced to 5-29% which is close to the results achieved with the original gain factors. The largest improvement can be seen for the GRACE-derived polar motion excitation time series of the oceans (19 percentage points), followed by Greenland (7 percentage points), Antarctica (5 percentage points) and the continental hydrosphere (3 percentage points).

GRACE gravity field models
After the successfully performed closed-loop simulation the novel filter effect reduction approach is applied to the GRACE GFZ RL06 time variable gravity field models. Figure 6 shows the global grid point gain factors derived from once-and twice-filtered GRACE GFZ RL06 equivalent water heights. Compared to the simulated gain factors these gain factors are often slightly smaller. Figure 7 shows the standard deviations of equivalent water heights over 13 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) derived from GRACE GFZ RL06 gravity field solutions using the DDK filter and by applying the DDK filter and the global grid point gain factors. The signal over the continents is significantly larger than the signal over the oceans (the colorbar of the figure is selected to display the differences over the oceans). One can see that the signal over the oceans is significantly decreased (up to 20% ) due to reducing the large leakage-in effect from the continental hydrosphere. Furthermore, significant changes over Antarctica are seen because the leakage effect of the larger mass changes in the bounder areas are reduced. The same phenomenon can be notice in Greenland. Looking at the GRACEderived polar motion excitation functions for the continental hydrosphere, oceans, Antarctica and Greenland in Fig. 8 one can see that due to the filtering the signal of the continental hydrosphere, Antarctica and Greenland is damped because the leakage-out effects are significantly larger than the leakage-in effects from the smaller signal of the oceans. In the oceans it is exactly the opposite the leakage-out effects are small and the leakage-in effects from the continental hydrosphere, Antarctica and Greenland are large. Therefore the oceanic mass effect is strengthened due to the filtering process. By applying our filter effect reduction approach one can counteract these phenomena, especially if the scaling factors are considered. The impact on polar motion excitations of the oceans ( RSD = 16% ), Greenland ( RSD = 16% ) and Antarctica ( RSD = 13% ) is large, while the impact on the hydrological excitations ( RSD = 3% ) is significantly smal ler.
Comparison with external solutions is performed exemplarily for the oceanic mass-related excitations because geophysical model results for excitations of the continental hydrosphere and the two continental ice sheets are characterized by higher uncertainties. We use the mass-related parts of the oceanic excitation series from the ocean model ECCO (Estimating the Circulation and Climate of the Ocean; ftp://euler  Table 1 Correlation coefficients and relative standard deviations (RSD) in percentage of the polar motion excitation functions for the continental hydrosphere χ H , oceans χ O , Antarctica χ A and Greenland χ G derived from simulated noisy GRACE gravity field solutions using the DDK filter or by applying the DDK filter in combination with the global grid point gain factors k org or k (inclusive scaling factors f i ), respectively, with respect to polar motion excitation functions derived from true GRACE gravity field solutions  .jpl.nasa.gov/sbo/oam_globa l/ECCO_kf079 .chi) and from the Earth System Modelling group at German Research Centre for Geosciences (ESMGFZ) which are based on data from the ocean model MPIOM [Max-Planck-Institute Ocean Model, Dobslaw and Dill (2018)]. The ESMGFZ time series allow to combine the effective angular momentum functions for the dynamic ocean (OAM) and for the barystatic sea-level (SLAM) to take into account the inflow of terrestrial water into the oceans. The trends of the ocean model results are not realistic due to the fact that both ocean models are Boussinesq model that conserve volume rather than mass. Mass conservation is achieved by adding a homogeneous shell of mass at every time step. Figure 9 shows the differences between the GRACE-based oceanic polar motion excitation functions with and without considering the filter effect reduction approach and the ocean model results. The trends in all time series are reduced. The corresponding correlation coefficients and relative standard deviations (RSD) are given in Table 2. By applying the global grid point gain factors the agreement of GRACE-derived and modelbased oceanic polar motion excitation functions can be improved by up to 8 percentage points. By taking into account the scaling factors the agreement can be improved even more by up to 15 percentage points. In general, the accordance with the ESMGFZ data is higher than with the ocean model results from ECCO.

Summary and conclusions
In this study we present a novel approach to reduce filter effects (attenuation and leakage) in GRACE-derived equivalent water heights and polar motion excitation functions. The advantages of this method are that it is independent from geophysical model information, works on a global grid point scale and can therefore be used for mass variation estimations of several subsystems of the Earth. The gain factors are derived from once-and twice-filtered GRACE gravity field solutions. Over the continents the gain factors are mostly positive because the larger land signals leak into   (1) only the DDK filter (green), (2) the DDK filter and the gain factors k (black) and (3) the DDK filter, the gain factors k and the scaling factors f i (red) with respect to the polar motion excitation functions derived from the ocean model ECCO (left) and from ESMGFZ data (right) the weaker ocean signals. Whereas over the oceans the gain factors are mostly negative to counteract the large leakagein effects of the continental hydrosphere.
To demonstrate the benefit of this method a closedloop simulation is performed. It can be shown that due to the application of the global grid point gain factors the discrepancies with respect to the true data can be significantly reduced, especially for the oceans, Greenland and Antarctica. The uncertainties in GRACE-derived polar motion excitations can be reduced from 12-48% to 5-29% by taking into account not only the gain factors but also the scaling factors which counteract the damping due to the filtering. The largest improvement is observed for the oceanic excitations (19 percentage points), followed by the polar motion excitations of Greenland (7 percentage points), Antarctica (5 percentage points) and continental hydrosphere (3 percentage points).
We applied the filter effect reduction approach to real GRACE gravity field solutions from GFZ. The agreement with ocean model results for the oceanic mass-related excitation functions could be improved by up to 15 percentage points. Thus the time series taking into account the filter effect reduction approach are closer to the truth.
In this study we demonstrated that our method allows to determine improved polar motion excitation functions for the continental hydrosphere, oceans and the continental ice sheets Antarctica and Greenland. This filter effect reduction may also be applied to estimate improved mass changes for specific regions, such as continents, oceans, ice sheets, or river basins. The achieved improvement will, however, be different for each region according to its size, its location and the performance of the DDK filter in the respective area.