Decorrelation of GRACE Time Variable Gravity Field Solutions Using Full Covariance Information

: In this study the feasibility and performance of time variable decorrelation (VADER) ﬁlters derived from covariance information on decadal Gravity Recovery and Climate Experiment (GRACE) time series are investigated. The VADER ﬁlter is based on publicly available data that are provided by several GRACE processing centers, and does not need its own Level-2 processing chain. Numerical closed loop simulations, incorporating stochastic and deterministic error budgets, serve as basis for the design of the ﬁlter setup, and the resulting ﬁlters are subsequently applied for real data processing. The closed loop experiments demonstrate the impact of temporally varying error and signal covariance matrices that are used for the design of decorrelation ﬁlters. The results indicate an average reduction of cumulative geoid height errors of 15% using time-variable instead of static decorrelation. Based on the simulation experience, a real data ﬁltering procedure is designed and set up. It is applied to the ITSG-Grace2014 time variable gravity ﬁeld time series with its associated full monthly covariance matrices. To assess the validity of the approach, linear mass trend estimates for the Antarctic Peninsula are computed using VADER ﬁlters and compared to previous estimates from both, GRACE and other mass balance estimation techniques. The mass change results obtained show very good agreement with other estimates and are robust against variations of the ﬁlter strength. The DDK decorrelation ﬁlter serves as main benchmark for the assessment of the VADER ﬁlter. For comparable ﬁlter strengths the VADER ﬁlters achieve a better de-striping and deliver smaller formal errors than static ﬁlters like the DDK.


Introduction
After its mission ended in 2017, more than a decade of GRACE (Gravity Recovery and Climate Experiment; [1]) monthly gravity field solutions are available. This unique time series allows estimation of geophysical processes like ice mass change [2,3] or hydrological features [4,5]. The most important goal is the optimum separation of signal and noise to receive high spatial resolution of target signals. This study focuses on post-processing strategies, in the form of decorrelation or regularization filters, computed from signal and error covariance information, of GRACE time variable gravity field solutions. The strategy is based on previously published work on DDK decorrelation filters by [6,7]. In [8] the processing concept proposed by [6] is continued and amended by using non-stationary error covariance information from their in-house GRACE processing, which is not publicly available and thus cannot be applied to gravity time series of other GRACE processing centers. Contrary to the approach presented within this study, the signal variance they use is computed in an iterative process from the GRACE observations themselves. As of today, their published series of filtered spherical harmonic coefficients is not maintained anymore, but is available only for the period from February 2003 until November 2010.
In contrast to [8], the strategy proposed by [6] applies simplifications (stationary filter characteristics) within the process of computing the filter kernels assuming a constant error covariance and signal variance information throughout the GRACE mission lifetime. Furthermore, the error covariance is approximated by a synthetic error covariance model derived from GRACE orbit information, which might lead to a significant decrease in recoverability of geophysical signals throughout the GRACE mission lifetime.
In this study we investigate the feasibility and potential of creating non-stationary decorrelation filters based entirely on publicly available data, which can be applied to any GRACE temporal gravity time series that is also complemented by full covariance information. Other key differences between the approach by [8] and our method are (1) the different set of error covariance information, and (2) a different process to determine the month-to-month signal variance. Additionally, we also investigate the impact of signal attenuation and induced phase shift of target signals by the filtering applied, putting focus on the different design options for the decorrelation filter and comparing with the performance of empirical decorrelation [9] and standard Gaussian filtering. This comparison intends to allow selecting a suitable post-processing strategy for the respective application.
Recently more and more processing centers make their month-to month error covariance matrices available. As of now (July 2018), to our knowledge the ITSG-Grace2014 [10] and ITSG-Grace2016 [11] time series from TU Graz, the time series from the Center of Space Research at the University of Texas [12], and the GFZ time series [13] are provided with full month-to-month error covariance matrices. The increasing number of processing centers distributing their error covariance matrices allows for applying decorrelation filters like the VADER filter to a growing selection of GRACE time series and to better exploit the signal content of the GRACE data.
The concept of the decorrelation process is introduced in Section 2. Section 3 comprises the introduction of the closed-loop simulation environment (Section 3.1) and the determination of a favorable filter design (Section 3.2). The impact of filter induced phase shifts of retrieved basin scale mass change signals in the closed-loop environment are evaluated in Section 4. Section 5 presents the filter design that is applied to real GRACE data and compared to other post-processing methods. The last section provides the main conclusions and a short outlook.

Decorrelation Strategy and Setup of Experiments
As mentioned in Section 1, we build our experiment on the decorrelation approach by [6,7] for filtering GRACE time variable gravity field solutions. The main relation between unfiltered x and filtered x α gravity field solutions, represented by vectors of spherical harmonic coefficients, is given by with the corresponding normal equation matrix N (inverse of the error variance-covariance matrix), the inverse signal variance matrix M, and the scaling factor α. These three entities form the filter matrix W α . A detailed discussion of this approach can be found in [6,14]. As mentioned before, we intend to avoid any simplification for computing the filter matrix W α , especially with respect to the error covariance matrix. One should note that the gravity field solutions represented by x and x α are residual fields with respect to a static field, expressing temporal gravity signals averaged over the time span of one month.
With the scaling factor α the filter strength is adjusted. The larger the α, the stronger the filtering is constrained towards the a-priori signal variance, resulting in smoother filtered gravity fields. The filter strength can be approximated similarly to [7] by a mean half-weight radius of the filter kernel. It depends on location (latitude and longitude; [15]) and on the direction (cf. Figure 1). The characteristics of the GRACE error covariances typically lead to a smaller half-weight radius in a north-south than in an east-west direction (deviation between 50% and 60%). The reason can be found in the way the gravity field is sampled by the along-track ranging mission GRACE, in an orbit with an inclination of 89.5 degrees. This set-up results in higher sensitivity in a north-south than in an east-west direction. The absolute values of the parameter α are depending on the relative values of the signal and error covariance matrices. Therefore the values for α presented in this study should not be generalized, since other relations of magnitude between signal and error variances can occur. However, the relative differences of α matter, since they are steering the filter strength. found in the way the gravity field is sampled by the along-track ranging mission GRACE, in an orbit with an inclination of 89.5 degrees. This set-up results in higher sensitivity in a north-south than in an east-west direction. The absolute values of the parameter are depending on the relative values of the signal and error covariance matrices. Therefore the values for presented in this study should not be generalized, since other relations of magnitude between signal and error variances can occur. However, the relative differences of matter, since they are steering the filter strength.
(a) (b) Figure 1. (a) Normalized 2D representation of the filter Kernel in the spatial domain, exemplarily shown for a geographical longitude of λ = 0 degrees, and latitude ϕ = 30 degrees north (center of figure). (b) Transection through the normalized smoothing kernel. In blue the east-west section and in red the north-south section through the center of the kernel are shown. Figure 1a shows an example of a 2D representation in the spatial domain for a convolution kernel derived from the filter matrix . x-and y-axis indicate the distance in kilometers from the kernel center which is for the given example a geographical longitude λ = 0 (Greenwich meridian) and latitude ϕ = 30 degrees north. One can observe anisotropic behavior including negative side lobes especially in a north-south direction. This side-lobe characteristic is illustrated also in Figure 1b, showing the traverse through the center in a north-south and an east-west directions for the same kernel as presented in Figure 1a. The negative side lobes in a north-south direction have an edge preserving or sharpening effect for the along track direction. The absence of strong negative side lobes in an east-west direction induces stronger smoothing in the cross track direction. In combination this causes a reduction of the typical GRACE striping pattern [16]. The filter kernel is shown in spatial domain only for illustrative purposes. From the perspective of computational costs, the application of the filter in frequency domain is numerically more efficient than in the spatial domain by several orders of magnitude. As Equation (1) shows, the filtering process is a matrix vector multiplication in frequency domain. This is a benefit of the spherical harmonic representation of the gravity field. In spatial domain one would need to evaluate and perform the convolution on the whole grid for every single grid point separately due to the anisotropy and inhomogeneity of the filter.
To compare different decorrelation filters with each other, the half weight radius of the convolution kernel can be estimated. One option, also proposed by [7], is to estimate the half weight radius of the convolution kernel at specific latitudes in all cardinal directions. This method has some drawbacks, since such an approach is always based on the choice of samples, which are finite, and does not necessarily represent a global picture. Nevertheless we decided to accept these drawbacks and adopt the method applied by [17] for all eight DDK filter strengths, to allow comparing the decorrelation filters in a consistent manner. Accordingly we evaluate all monthly filter kernels at the same locations (ϕ = 0, 30, and 60 degrees north at a longitude of λ = 0 degrees) and estimate the half weight radius in all four cardinal directions. The resulting average radii are summarized in Table 1 Figure 1a shows an example of a 2D representation in the spatial domain for a convolution kernel derived from the filter matrix W α . x-and y-axis indicate the distance in kilometers from the kernel center which is for the given example a geographical longitude λ = 0 (Greenwich meridian) and latitude ϕ = 30 degrees north. One can observe anisotropic behavior including negative side lobes especially in a north-south direction. This side-lobe characteristic is illustrated also in Figure 1b, showing the traverse through the center in a north-south and an east-west directions for the same kernel as presented in Figure 1a. The negative side lobes in a north-south direction have an edge preserving or sharpening effect for the along track direction. The absence of strong negative side lobes in an east-west direction induces stronger smoothing in the cross track direction. In combination this causes a reduction of the typical GRACE striping pattern [16]. The filter kernel is shown in spatial domain only for illustrative purposes. From the perspective of computational costs, the application of the filter in frequency domain is numerically more efficient than in the spatial domain by several orders of magnitude. As Equation (1) shows, the filtering process is a matrix vector multiplication in frequency domain. This is a benefit of the spherical harmonic representation of the gravity field. In spatial domain one would need to evaluate W α and perform the convolution on the whole grid for every single grid point separately due to the anisotropy and inhomogeneity of the filter.
To compare different decorrelation filters with each other, the half weight radius of the convolution kernel can be estimated. One option, also proposed by [7], is to estimate the half weight radius of the convolution kernel at specific latitudes in all cardinal directions. This method has some drawbacks, since such an approach is always based on the choice of samples, which are finite, and does not necessarily represent a global picture. Nevertheless we decided to accept these drawbacks and adopt the method applied by [17] for all eight DDK filter strengths, to allow comparing the decorrelation filters in a consistent manner. Accordingly we evaluate all monthly filter kernels at the same locations (ϕ = 0, 30, and 60 degrees north at a longitude of λ = 0 degrees) and estimate the half weight radius in all four cardinal directions. The resulting average radii are summarized in Table 1 for different DDK and VADER filter strengths. The numbers for the VADER filter are derived as mean of all monthly solutions of a time series for the respective filter strength. Table 1. Mean Gaussian radii estimated from filter kernels at 0, 30, and 60 degrees latitude over 124 months (2003-02 until 2013-12) for the ITSG-Grace2014-VADER filter. Radii for the DDK results are computed as simple mean from numbers given in [17] for 0, 30, and 60 degrees latitude.

Results from Closed-Loop Simulations
This section describes the simulation environment used for the closed-loop experiments (Section 3.1) and assesses the different design options of the filter to derive a favorable filter design (Section 3.2).

The Closed-Loop Simulation Environment
In order to perform the closed-loop simulations, we employ the closed-loop mission simulator described by [18] to simulate GRACE-like observations and to calculate monthly GRACE-like time variable gravity field solutions. Acceleration differences along the line of sight between the two satellites serve as observations. GRACE-like white noise with standard deviation of 2 × 10 −9 m/s 2 is added to the observations as stochastic error on the level of accelerations. A discussion of different noise models, such as white noise and colored noise, can be found in [14]. This GRACE-like white noise can be accounted for being a good first approximation for the stochastic error budget. This is confirmed by the matching magnitude of the formal error degree root mean square (RMS) and the degree RMS of the error derived from residual ITSG-Grace2014 solutions from spherical harmonic (SH) degree 20 onwards. The input data for the observations are almost 12 years of 6 hour hydrological, ice, and solid earth signal (HIS) data from the ESA Earth system model [19]. The spherical harmonic degree of expansion used within the simulations is 120. To simulate non-tidal atmospheric and ocean signal (AO) de-aliasing errors, AO observations are added with a weighting of 10% and 20% to the HIS components at the level of normal equations. The choice of 10% and 20% is derived from the investigations presented by [20]. Another option, not pursued within this study, would be the use of AO errors described in [21]. Ocean tide (OT) de-aliasing errors are accounted for by ocean tide model differences between FES2004 [22] and EOT08a [23] in the same way as the AO de-aliasing errors. The AO and OT error budget are considered as coarse upper bounds of error budget contributions.
Other key inputs for the simulations are the orbits. For the presented study we use GRACE orbits from [24] for the period between January 2003 and November 2014. Obviously the timeframe of the orbits and the Earth model do not match. For the investigations regarding filter design and performance this does not cause any disadvantage. The simulated time series just contains geophysical signal from another decade but the conclusions drawn are nevertheless meaningful and valid as if the timeframe was matching because the orbit positions are not used as observations, but only the gravitational acceleration differences computed are at the orbit positions from the underlying time variable gravity field.
Main outputs from the simulations are the monthly sets (timeframe is defined by the respective GRACE orbits) of SH coefficients by solving the corresponding full normal equation systems (NEQs). The input N in Equation (1) is derived from the NEQs of the simulation runs.
The closed-loop simulations are run for each of the 143 months of the timeframe investigated. The different post-processing strategies are compared on the basis of the residuals with respect to the true reference, which is the mean HIS signal of each individual month.

Determination of a Favorable Filter Design
With the closed-loop simulator the effect of employing different populations of the error covariance matrices and different designs of the signal variance matrices, building the filter matrix W α , are tested. They are for N (cf. Equation (1)) in particular full, order-wise block diagonal, and diagonal population. For M the effect of using signal variance composed by the true reference for each month, the variance of the annual mean (calendar year), a calendar month mean, and an overall mean (in time) are evaluated. Monthly mean HIS signal of the Earth model is the basis for deriving the signal variance used within the simulation runs. The true reference is rejected for further application since it is not known in real world environments and since applying the same signal in a circular way is not advisable. Accordingly, more generalized representations of the signal variance are applied for further processing. In addition, the respective signal variance is not applied directly, but as a degree dependent approximation, leading to a population of the signal variance matrix only along its diagonal.
The SH degree l dependent approximation of the signal variance σ M reads Since the approximate signal variance is now only degree dependent, it is equivalent to a Kaula rule [25]. Figure 2 illustrates the variability of the a and b coefficients in Equation (2), derived from the 143 month time series of the HIS dataset. a and b coefficients are estimated from each month, and the median for each calendar month is indicated in the figure. For comparison, not only the a and b coefficients derived from the Earth model are shown, but also the ones derive from 143 months of real GRACE data, filtered with the DDK4 filter. The pre-filtered datasets are downloaded from the International Centre for Global Earth Models (ICGEM) web portal (http://icgem.gfz-potsdam.de/ ICGEM/). A similar analysis regarding the variability of signal variance can be found in [15]. The amplitude of the a coefficients from the DDK4 filtered ITSG-Grace-2014 solution is higher than the one derived from the model (about factor 2; illustrated by dashed red line). One reason is that the signal contained in the filtered GRACE solutions also contains-contrary to the coefficients derived from the Earth model-signal from residual AO and OT de-aliasing errors. A second effect is the difference in signal strength itself; analysis showed that the Earth model tends to show smaller signal amplitudes than the real world situation [26]. The difference of a factor 2 for the a coefficient can be accounted for when selecting a suitable α value. The dashed red line reveals that the characteristics of the a coefficients derived from the Earth model data and the DDK4 filtered ITSG-Grace2014 data is, apart from the scaling issue, very similar. The b coefficients, representing the exponent in Equation (2), agree very well (apart from the difference of the local maxima in June and July) with the ones derived from the DDK4 filtered solution. But one should keep in mind that the underlying datasets do not cover the same time period. Accordingly some minor deviations must be expected although both results represent decadal average values. Due to the fact that decadal average values are used for building the signal variance for the VADER filter, a bias towards the Earth models characteristics is avoided. The impact of using static, cyclo-stationary, and non-stationary signal variance is evaluated in detail later in this section. The impact of different populations of the error covariance matrices is analyzed in the following. Figure 3 shows cumulative geoid errors in mm for different VADER filters. The median from all 143 monthly solutions is shown for the respective filter strength, steered by the choice of the value, which is expressed by the Gaussian smoothing radius according to Table 1. The impact of using different setups of the inverse of the normal equation matrix in equation (1) is demonstrated. The smallest residual is found using the full information from the covariance matrix (red curve). Approximating the matrix and keeping only blocks containing the correlations among coefficients of the same order (order blocks) causes inferior recoverability of the original signal (black curve). Degrading the covariance matrix even more to a diagonal structure (blue curve), which implies not considering any correlation among coefficients, delivers the largest residual geoid error of the three options. For the diagonal setup, not only is the minimum residual achievable significantly higher than for the other two options, but also stronger filtering is necessary (larger average filter radius) to reach this minimum. However, still all VADER filter setups are able to achieve better signal recoverability than the Gaussian filter (cyan dashed line). Another way for determining the signal variance is the one embedded in the DMT-1 series [27], which estimates the signal variance from the solution itself, and is therefore on the one hand independent of external model data. On the other hand, separation of different geophysical signal constituents is not possible with such a technique. Using the spectral characteristics of specific target signals like HIS for building the filter allows tailoring of the regularization more precisely for the respective filter scenario.
The impact of different populations of the error covariance matrices is analyzed in the following. Figure 3 shows cumulative geoid errors in mm for different VADER filters. The median from all 143 monthly solutions is shown for the respective filter strength, steered by the choice of the α value, which is expressed by the Gaussian smoothing radius according to Table 1. The impact of using different setups of the inverse of the normal equation matrix N in equation (1) is demonstrated. The smallest residual is found using the full information from the covariance matrix (red curve). Approximating the matrix and keeping only blocks containing the correlations among coefficients of the same order (order blocks) causes inferior recoverability of the original signal (black curve). Degrading the covariance matrix even more to a diagonal structure (blue curve), which implies not considering any correlation among coefficients, delivers the largest residual geoid error of the three options. For the diagonal setup, not only is the minimum residual achievable significantly higher than for the other two options, but also stronger filtering is necessary (larger average filter radius) to reach this minimum. However, still all VADER filter setups are able to achieve better signal recoverability than the Gaussian filter (cyan dashed line).
The average smoothing radius should not be misinterpreted as the exact number of the spatial resolution of the final result. The average smoothing radius can, due to the way it is derived, be seen as an approximation of the spatial resolution. If the convolution kernel was Gaussian type shaped (meaning isotropic), the half weight radius would correspond exactly to the spatial resolution of the filtered solution [28]. The filter in general acts as a kind of low-pass filter omitting high frequency signals. Since the VADER filter is highly anisotropic and the indicated average smoothing radius is a global approximation, the real resulting spatial resolution cannot be described by a single number.
The results of the closed-loop simulations so far indicate receiving best performance by using fully populated error covariance matrices. The analysis regarding the approximated signal variance indicates a varying behavior of the signal variance in a cyclo-stationary manner. A stationary signal variance is expected to deliver inferior results as this state does not represent the real variability. Using true month-to-month signal variance is on the other hand most promising from a theoretical perspective, but not possible to realize in real world applications, because the true signal is simply not known. The average smoothing radius should not be misinterpreted as the exact number of the spatial resolution of the final result. The average smoothing radius can, due to the way it is derived, be seen as an approximation of the spatial resolution. If the convolution kernel was Gaussian type shaped (meaning isotropic), the half weight radius would correspond exactly to the spatial resolution of the filtered solution [28]. The filter in general acts as a kind of low-pass filter omitting high frequency signals. Since the VADER filter is highly anisotropic and the indicated average smoothing radius is a global approximation, the real resulting spatial resolution cannot be described by a single number.
The results of the closed-loop simulations so far indicate receiving best performance by using fully populated error covariance matrices. The analysis regarding the approximated signal variance indicates a varying behavior of the signal variance in a cyclo-stationary manner. A stationary signal variance is expected to deliver inferior results as this state does not represent the real variability. Using true month-to-month signal variance is on the other hand most promising from a theoretical perspective, but not possible to realize in real world applications, because the true signal is simply not known.
In the following an in-depth analysis of the impact of stationary or non-stationary error covariance and signal variance is presented. Figure 3 shows the median values of cumulative geoid error from residuals of the filtered set of solutions. Analyzing the scattering of the respective clusters of filtered solutions gives valuable insight into the properties and characteristics of the filter design. Figure 4 shows the individual monthly cumulative geoid errors in mm derived from the residuals (filtered solution minus true reference). The colored clusters indicate different filter strengths. The black x-marks indicate the median of the respective cluster in terms of cumulative geoid error (y-axis) and average smoothing radius (x-axis). Accordingly the x-marks correspond to the median curves like the ones shown in Figure 3. Figure 4a shows the results, building the VADER filter with static (stationary) error covariance and variable (cyclo-stationary) signal variance. Panel (b) shows results with both ingredients, signal and error covariance, with non-stationary properties (non-stationary N and cyclostationary M). Especially for the clusters around the local minimum (green to blue colored clusters) In the following an in-depth analysis of the impact of stationary or non-stationary error covariance and signal variance is presented. Figure 3 shows the median values of cumulative geoid error from residuals of the filtered set of solutions. Analyzing the scattering of the respective clusters of filtered solutions gives valuable insight into the properties and characteristics of the filter design. Figure 4 shows the individual monthly cumulative geoid errors in mm derived from the residuals (filtered solution minus true reference). The colored clusters indicate different filter strengths. The black x-marks indicate the median of the respective cluster in terms of cumulative geoid error (y-axis) and average smoothing radius (x-axis). Accordingly the x-marks correspond to the median curves like the ones shown in Figure 3. Figure 4a shows the results, building the VADER filter with static (stationary) error covariance and variable (cyclo-stationary) signal variance. Panel (b) shows results with both ingredients, signal and error covariance, with non-stationary properties (non-stationary N and cyclo-stationary M). Especially for the clusters around the local minimum (green to blue colored clusters) the scattering is significantly reduced. One can conclude that using non-stationary error covariance significantly improves the quality of filtered time series. Panel (c) indicates in black the performance of a stationary design of the VADER filter, meaning stationary M and N as used also by the DDK filter. For this scenario the average smoothing radius is the same for each month since the filter matrix W α is identical for each month. For comparison, the performance of a 350 km Gaussian filter in magenta and an empirical decorrelation filter [9] in green are shown (label S&W). the scattering is significantly reduced. One can conclude that using non-stationary error covariance significantly improves the quality of filtered time series. Panel (c) indicates in black the performance of a stationary design of the VADER filter, meaning stationary M and N as used also by the DDK filter. For this scenario the average smoothing radius is the same for each month since the filter matrix is identical for each month. For comparison, the performance of a 350 km Gaussian filter in magenta and an empirical decorrelation filter [9] in green are shown (label S&W).    Figure 4 by giving some performance metrics for the different post-processing methods. The average smoothing radii are indicated to allow a rough comparison of the respective filter strength. For the empirical S&W decorrelation filter a window size of 5 with a second degree polynomial was chosen. The minimum SH order for the first part of the filter cascade is 10 and for the second part of the filter cascade a 250 km Gaussian filter is applied. For details regarding the assessment of the empirical decorrelation refer to [9,14]. The lowest residuals are derived when using the VADER filter in non-stationary design. All results originate from the respective local minima of the different post-processing techniques. The impact of switching the properties of the error covariance matrix are significantly larger than for the signal variance (as long as it is chosen within reasonable bounds). The scattering is drastically reduced using time variable decorrelation, comparing the max values of the outliers scaling down from 6.84 and 7.14 to 0.44 and 0.43. The empirical decorrelation (S&W) does show better performance than the standard Gaussian filter, but is clearly outperformed by the VADER filter. Months showing large deviations from the median mostly show problems in data coverage and orbital configuration (short repeat cycles). The improvement of the median of cumulative geoid errors from residuals is about 15% switching from stationary to non-stationary filter design. For individual months this improvement is significantly higher as one can see from the single scattered results (Figure 4). The maxima are reduced by more than one order of magnitude. Figure 5 shows the performance of post-processing (VADER compared with Gaussian) within the closed-loop experiments for one individual month. Panel (a) shows SH degree RMS in mm geoid height of the residuals of the raw and the filtered solutions. Residual means in this case the difference between the true HIS signal, represented by the mean HIS input signal, as an approximation, for the simulation run of the respective month (black curve). To analyze the effect of deterministic error sources like ocean tide (OT) signal, an OT difference signal is added to the observations, leading to the residual SH degree RMS drawn in the dashed lines. As described in [18], characteristic RMS peaks, induced by incomplete OT de-aliasing affecting specific SH order bands, can be seen in the dashed gray line around SH degrees 46, 61, and 107. As one can conclude from the dashed red line, the OT signal is almost completely removed from the observations and the residuals are almost as small as for the case without deterministic errors. For the Gaussian filter, this is not the case. Here one can still see significantly increased residuals after filtering (dashed blue line). Another notable effect is the intersection of the solid blue and solid red line with the solid black line, the latter representing the average HIS signal for the specific month. The VADER filter allows signal recovery (positive signal to noise ratio) up to SH degree 55 instead 30, using the Gaussian filter. representing the average HIS signal for the specific month. The VADER filter allows signal recovery (positive signal to noise ratio) up to SH degree 55 instead 30, using the Gaussian filter.  Figure 5b shows cumulative geoid errors in mm for different filter radii computed from the residuals of all 143 monthly solutions (the median is shown for each filter strength) with respect to the corresponding HIS signal filtered with VADER (solid lines) and Gaussian (dashed lines). Complete de-aliasing (blue) is compared with effects of ocean tide model differences (red) and of 20% of the non-tidal atmospheric and oceanic (AO) signals (black). The minima of the solid curves are always lower (for same colors) than the dashed lines and correspond to a shorter average smoothing radius. The increase of the residuals is stronger in the case where deterministic OT error contributions are added instead of the 20% AO error. This circumstance is of course related to the choice of error contributions and will change using different assumptions regarding the error budgets. The minimal residuals for both filter techniques are shifted to higher average smoothing radii, implying the need for stronger filtering, by inducing additional errors. This effect is stronger for the Gaussian filter than for the VADER filter which shows hardly any noticeable shifts in the smoothing radius. The ability to better separate between target signals and noise is one of the key advantages of the VADER filter compared to the Gaussian filter.  Figure 5b shows cumulative geoid errors in mm for different filter radii computed from the residuals of all 143 monthly solutions (the median is shown for each filter strength) with respect to the corresponding HIS signal filtered with VADER (solid lines) and Gaussian (dashed lines). Complete de-aliasing (blue) is compared with effects of ocean tide model differences (red) and of 20% of the non-tidal atmospheric and oceanic (AO) signals (black). The minima of the solid curves are always lower (for same colors) than the dashed lines and correspond to a shorter average smoothing radius. The increase of the residuals is stronger in the case where deterministic OT error contributions are added instead of the 20% AO error. This circumstance is of course related to the choice of error contributions and will change using different assumptions regarding the error budgets. The minimal residuals for both filter techniques are shifted to higher average smoothing radii, implying the need for stronger filtering, by inducing additional errors. This effect is stronger for the Gaussian filter than for the VADER filter which shows hardly any noticeable shifts in the smoothing radius. The ability to better separate between target signals and noise is one of the key advantages of the VADER filter compared to the Gaussian filter.

Impact of Post-Processing Methods on the Phase of Seasonal Signals with the Closed-Loop Environment
Another aspect of importance, is the impact of the applied post-processing method not only on the amplitude, but also on the phase of geophysical signals. For this purpose seasonal signals (linear, annual, semi-annual, and 161-day) are estimated from a 143 month closed-loop time series in a least squares procedure as described in [14]. An analysis method as described by [29], using the Hilbert transform of the time series of seasonal signals, is applied to estimate the phase difference between the true solution and the respective filtered solutions. In this study we analyze selected river basins, taken from the STN-30p dataset [30]. For the regional integration a simple zero/one kernel is applied (zero outside the region of interest and one inside). The performance of four post-processing strategies or filter designs is investigated. Table 3 gives a summary of the VADER filter with non-stationary error covariance and cyclo-stationary signal variance (VADER), the VADER filter with non-stationary error covariance and stationary signal variance (VADER med.), the empirical S&W decorrelation (S&W), and the Gaussian filter approach (Gauss). The label VADER med. is chosen since the signal variance is computed by taking the median of the cyclo-stationary signal variance, thus leaving a stationary signal variance. The main reason for phase shifts of mass signals derived on basin scale is leakage. If two neighboring basins embed signal, the signal is leaking over the basin boundary, thus influencing the estimated mass signals of each other. If the signals show, for example, a seasonal variation and the phase of this variation is different for the two neighboring basins, the estimated amplitude and phase of the two basins is correlated with each other and possibly altered from the true state. Investigations on such effects and how to deal with them are discussed, for example, in [31]. The best situation is of course using a filter method, which minimizes the effect of leakage and potential phase shifts already in the first place. The VADER filter causes, as Table 3 shows, the smallest phase shift for all analyzed river basins. The VADER med filter with stationary instead of cyclo-stationary signal variance performs almost equally well. For Ganges and Parana, the cyclo-stationary signal variance delivers the smallest phase shift RMS compared to the stationary signal variance, whereas for Amazon and Mississippi the results are almost equal. The main impact on the VADER filter performance can be attributed to the use of the correct month-to-month error covariance matrices. For the Mississippi river, the phase shift in general is relatively small since it is fairly isolated from other sources showing high seasonal amplitudes. The opposite holds for the Parana River, neighboring the Amazon basin. The Amazon basin itself show relatively small phase shifts since it is very large in size and has strong signal amplitudes. The larger the basins and the stronger the signals are, the smaller the influence by the surrounding river basins becomes. As an intermediate conclusion to be drawn from the closed-loop experiments, the VADER filter in non-stationary design seems to have to the best target signal (amplitude and phase) retrieval capability compared to all other post-processing methods evaluated in this study. The next section intends to substantiate these findings by testing the favorable VADER filter design within a real data environment.

Application to Real Data
Different post-processing methods were applied to monthly gravity field solutions of ITSG-Grace2014 [10], for which both the monthly coefficient solutions and the corresponding full variance-covariance matrices are available. As an example, monthly solutions of April and September 2004 were picked for further detailed analysis. April 2004 is considered to represent a nominal state of the GRACE mission, whereas September suffers from a short repeat cycle, resulting in non-ideal data coverage. The DDK filtered series used for comparison were obtained as a ready to use pre-filtered product from the ICGEM web portal (http://icgem.gfz-potsdam.de/ICGEM/). Figure 6 shows in the left column post-filter results of a residual GRACE solution. The static field used for computing a residual field is GOCO02S [32]. This conclusion is also underlined by investigating the impact of different post-processing techniques, especially the use of static (DDK) versus time-variable (VADER) covariance information, on the stability of long-term mass trend estimates. Based again on ITGS-Grace2014 data [10] for the period 2003 to 2008, mass trends have been estimated for the Antarctic Peninsula. This region is defined by basin number 24, 25, 26, and 27, following [33]. Different DDK and VADER filter strength are evaluated. Apart from these different filters, all processing steps are the same to ensure comparability. Figure 7 shows the mass trend results in dependence of the filter strength, expressed by the average Gaussian smoothing radius, for the different DDK (black) and VADER (blue) filters. Additionally the error bars are shown, demonstrating, as expected, smaller error bars for the stronger filters. As a reference, in purple color the 2012 ice sheet mass balance inter-comparison exercise (IMBIE) mass balance estimate based on the gravity method, in green color the average of all methods applied and compared in [2] are plotted. Applying different VADER filter strengths, very similar All four filter strategies deliver results that allow identifying geophysical features. However, comparing, for example, the signal in Greenland or the Gulf of Alaska reveals higher amplitudes and less leakage into the ocean using the VADER filter instead of Gaussian or empirical decorrelation. For the solution in September 2004 the situation is very different. With the exception of the VADER filtered solution, all other solutions show significant residual striping, or noise in low and mid latitudes. The VADER filter is able to account for the particular circumstances present in such a short repeat cycle orbit since it is data driven, presuming that these circumstances are correctly reflected in the error covariances N of the respective month. The empirical S&W decorrelation filter is data driven as well, but relies on certain assumptions regarding the correlation within SH coefficients, which are not met in this particular month. The Gaussian filter delivers a comparably smooth result since it is not data driven. But compared to the VADER filtered result it shows significant signal attenuation, visible, for example, by comparing the signals in Greenland. For the DDK4 filtered solution comparably strong residual striping is visible in the solution for April 2004. Since the pre-filtered products, provided via the ICGEM web portal are filtered all with the same DDK filter kernels, this circumstance is not surprising. Analysis of the different available time series with month-to-month error covariance matrices revealed significant differences in the correlation structure between different time series for same monthly solutions. Applying one generalized or uniform error correlation structure for different time series, using different processing approaches, bears the risk of introducing significant deviations and artefacts into the filtered solution. This issue should be kept in mind when using pre-filtered DDK solutions from the ICGEM portal. Upon the tests we conducted, we strongly recommend to use the respective error covariance matrices for the time series to be filtered.
This conclusion is also underlined by investigating the impact of different post-processing techniques, especially the use of static (DDK) versus time-variable (VADER) covariance information, on the stability of long-term mass trend estimates. Based again on ITGS-Grace2014 data [10] for the period 2003 to 2008, mass trends have been estimated for the Antarctic Peninsula. This region is defined by basin number 24, 25, 26, and 27, following [33]. Different DDK and VADER filter strength are evaluated. Apart from these different filters, all processing steps are the same to ensure comparability. Figure 7 shows the mass trend results in dependence of the filter strength, expressed by the average Gaussian smoothing radius, for the different DDK (black) and VADER (blue) filters. Additionally the error bars are shown, demonstrating, as expected, smaller error bars for the stronger filters. As a reference, in purple color the 2012 ice sheet mass balance inter-comparison exercise (IMBIE) mass balance estimate based on the gravity method, in green color the average of all methods applied and compared in [2]

Conclusions
Based on the results of this study, we conclude that is feasible to construct a powerful postprocessing filter from full GRACE covariance information. We regard the VADER filter approach to be a valid and performant post-processing strategy based on publicly available information only. This makes it easy to be applied to different GRACE time series produced by various GRACE processing centers. The statement by [8] regarding the suitability of such a filter approach can be

Conclusions
Based on the results of this study, we conclude that is feasible to construct a powerful post-processing filter from full GRACE covariance information. We regard the VADER filter approach to be a valid and performant post-processing strategy based on publicly available information only. This makes it easy to be applied to different GRACE time series produced by various GRACE processing centers. The statement by [8] regarding the suitability of such a filter approach can be confirmed by the performance with respect to short repeat cycles and induced phase shifts of hydrological signals. Evaluating the closed-loop simulations as well as real data filtering one can clearly observe better performance in terms of signal recoverability than with all other filters tested in this study. Since the VADER and DDK filter series follow deviant philosophies, stationary versus non-stationary design, a difference in the performance of the two filters must be expected. The VADER filter makes use of as much information available including level 2 error covariance information, whereas the DDK filter is based on more generalized assumptions using synthetic and static error covariance information.
The real data experiments comparing VADER and DDK show better performance for the VADER filter. The introduction of explicitly time variable filters does take into account the varying quality of the monthly solutions. Static filters in contrast do not take into account different orbital heights and therefore changing resonances, varying ground track patters and related changing spatial coverage of observations, as well as instrument degradation. The results show clearly that the assumption of a constant correlation structure is only partially valid and not optimal.
Special attention should be drawn when estimating trends and seasonal variations (e.g., annual or semiannual signals) with time series containing different signal and noise contents from month to month. The VADER filter accounts for varying noise content and thus performs significantly better than other filter strategies especially in critical months with decreased data coverage and/or quality. According to the assumed signal variances it also accounts for seasonally varying geophysical signal content. We assumed the 12 calendar month cyclo-stationary signal variances to be suitable for the VADER decorrelation approach. Data from models of course does not perfectly represent the truth, and wrong assumptions on the signal variances might affect trend, amplitude, and phase estimates. A static signal variance on the other side also does not reflect the true signal variance. Results in the closed-loop experiments show significantly larger residuals using static signal variance (although the error covariance is the main driver for the filter performance).
The main advantage of the presented VADER filter is its robustness and applicability on any time series providing full covariance information. Since the VADER filter is data driven, we expect advantages also with respect to the results of the GRACE-Follow on (GRACE-FO) gravity field mission, which was successfully launched in May 2018. Compared to GRACE, whose K-band microwave system had a ranging accuracy of about 1 µm, GRACE-FO has an improved inter-satellite ranging accuracy in the order of 40-50 nm due to laser interferometry. The VADER filter concept with its improved de-striping capabilities might help to fully exploit this improved ranging accuracy. Additionally, one of the main goals of GRACE-FO is to extend the 15 years of time series obtained by GRACE. Linking these two mission together, their different error characteristics, as they are reflected in the VADER filtered temporal gravity time series, will have to be taken into account, especially if consistent long-term trend estimates are to be derived.
With the GRACE-FO mission in place now, we encourage all processing centers to deliver time series with full covariance information to boost the exploitation of time variable gravity field solutions.

Conflicts of Interest:
The authors declare no conflict of interest.