Satellite Gravity Field Recovery Using Variance‐Covariance Information From Ocean Tide Models

Monthly gravity field recovery using data from the GRACE and GRACE Follow‐On missions includes errors limiting the spatial and temporal resolution of the estimated gravity fields. The major error contributions, besides the noise of the accelerometer instruments, arise from temporal aliasing errors due to imperfections in the non‐tidal atmospheric and oceanic de‐aliasing models and ocean tide models. We derive uncertainty information for the eight major tidal constituents from five different ocean tide models and introduce it into the gravity field recovery process in terms of a constrained normal equation system while expanding the parameter space by additional tidal parameters to be adjusted. We prove the effectiveness of the ocean tide variance‐covariance information through realistic simulations and we assess its potential based on microwave and laser interferometer observations from the GRACE Follow‐On mission. We show that errors are reduced by more than 20% ocean wRMS for a Gaussian filter radius of 300 km if uncertainty information for ocean tides is considered and stochastic modeling of instrument errors is applied, compared to the latest GFZ release 6.1. Our results also show the limited visibility of the effectiveness of the ocean tide variance‐covariance information due to the dominating error contribution of non‐tidal atmospheric and oceanic mass variations. Additionally, we investigate the option of estimating ocean tide parameters over a 1‐year period while including ocean tide uncertainty information in order to improve ocean tide background modeling.


Earth and Space Science
HAUK ET AL.

10.1029/2023EA003098
2 of 26 low-low-SST (ll-SST) systems.ll-SST observations on GRACE were based on K/Ka-band ranging (KBR) between two identical satellites in identical near-polar orbits at a separation of about 200 km operating at an initial altitude of 500 km.On GRACE-FO, which has an identical orbit configuration to the GRACE mission, ll-SST observations are operated by two independent systems in parallel.In addition to KBR, laser ranging interferometry (LRI, Kornfeld et al., 2019;Ghobadi-Far et al., 2020) is used as a technology demonstrator with increased nanometer precision compared to micrometer precision KBR (V.Müller et al., 2022).However, the potential of the LRI cannot be fully exploited since errors in sub-optimal modeling of high-frequency geophysical phenomena, such as non-tidal atmospheric and oceanic mass variations as well as ocean tide signals, dominate the total error budget of the gravity field recovery process in terms of temporal aliasing errors (Flechtner et al., 2015).
GRACE and GRACE-FO gravity field solutions are typically provided at a monthly time scale.High-frequency, (i.e., periods of a few days or even a few hours) geophysical mass variations can contaminate the monthly mean solutions due to undersampling and the so-called temporal aliasing error.Hence, processing strategies depend on background models to subtract the gravimetric high-frequency variability.Together with the anisotropic error behavior resulting from only along-track inter-satellite ranging, the typical spatial error structure in the GRACEand GRACE-FO gravity field maps arises in terms of meridional stripes.Besides temporal aliasing errors, the ACC instrument's error greatly contributes to the total error budget, especially in the longer wavelength spectrum of the time-varying gravity field (Flechtner et al., 2015).In GRACE-FO, this error source is even more prominent due to the deficiency of the ACC sensor on the second GRACE-FO satellite (denoted as GF2 in the following), resulting in a significant degradation of the ACC observations.For this reason, the GF2 ACC data need to be replaced by synthetic data, the so-called transplant data (Harvey et al., 2022), officially generated by the GRACE-FO Science Data System (SDS).The SDS transplant data is derived from the first GRACE-FO satellite (also denoted as GF1) ACC measurements by applying time shift and attitude corrections as well as model-based residual accelerations due to thruster firings on GF2.The transplant method helps to keep the level of uncertainty of the ACC observations at a reasonable level.However, the noise of the ACC system on GRACE-FO is still assumed to be larger than on GRACE before the thermal control was switched off in 2014 (Bandikova et al., 2019).
One standard procedure to mitigate striping errors in the recovered gravity fields is the a posteriori application of de-striping and filtering techniques (e.g., Horvath, 2017;Kusche, 2007;Swenson & Wahr, 2006).As a result, the adjusted gravity field models are smoothed and reduced in amplitude, implicating a loss of potential gravimetric information.Therefore, the constant improvement of background model information, such as ocean tide models, that are used to predict tidal mass variability, is indispensable for improving gravity field solutions.
Ocean tide models typically show larger uncertainties in polar regions, due to the coverage of satellite altimetry missions, such as TOPEX/POSEIDON, ERS-1, ERS-2 or Jason-1, which do not provide adequate data in high latitudes.Therefore, ocean tide models show largest error signals around Antarctica and the Arctic ocean (Andersen et al., 2023;Stammer et al., 2014).Moreover, regions covered with sea ice provide limited and reduced-quality tide gauge data, and regions with shallow waters suffer from land-contamination of along-track altimetry influencing the retrieval of sea level in the coastal region (Hart-Davis et al., 2021).These facts lead to increased errors in coastal regions and over continental shelves in the ocean tide models.Figure 1 shows the error characteristics in terms of complex ocean tide amplitude ζ, derived by the root sum square (RSS) of the standard deviations δ(ζ s ) (cf.Equation 1 in Section 2) of the eight major tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , and Q 1 ).The complex ocean tide amplitude ζ s,j for ocean tide model j and partial tide s combines the in-phase component as the real part and the quadrature component as the imaginary part of the tidal oscillation into a single complex quantity.Therefore, it considers amplitude differences, but also phase differences between individual ensemble members.These modeling errors cause a large part of the ocean tide aliasing effect in the satellite gravity field retrieval process (Stammer et al., 2014).
One way to improve ocean tide models for satellite gravity field recovery is analyzing GRACE and GRACE-FO data with co-estimation of the tidal gravity field parameters together with the nominal gravity field coefficients.This method has been used by the Institute of Geodesy at Graz University of Technology (ITSG) for the calculation of monthly GRACE and GRACE-FO products since their ITSG-Grace2018 release (Kvas et al., 2019).The method makes use of the entire length of GRACE data (more than 15 years), which is necessary to separate the different tidal constituents from their alias periods having much longer durations, for example, ∼7.5 years for K 1 and ∼3.7 years for K 2 , compared to the daily GRACE orbits (Mayer-Gürr et al., 2011).The estimated tidal HAUK ET AL.In this study, we reduce tidal aliasing errors in an alternative manner using uncertainty information from ocean tide models.These uncertainties pertain to differences derived from an ensemble of five ocean tide models (cf.Section 2).The idea of accounting for background model errors to find a more appropriate stochastic model for GRACE/GRACE-FO observations has already been discussed in Kvas and Mayer-Gürr (2019).The authors derived uncertainty estimates for the non-tidal atmospheric and oceanic de-aliasing product instead of ocean tides and investigated the consideration of these model uncertainties by co-estimating constrained model corrections, or by augmenting the covariance matrix of observations.The latter is also discussed in Abrykosov et al. (2021) using the identical ocean tide variance-covariance information presented in the scope of our study.The authors demonstrated a potential benefit in reducing tidal aliasing errors in GRACE/GRACE-FO-like gravity field recovery.They described an approach where the signal due to ocean tides is introduced into the system in terms of observations while propagating covariance information from ocean tide model errors directly to the level of observations.In this approach, no specific ocean tide parameters are setup, that is, the functional model is not extended by corresponding ocean tide parameters to be estimated.However, their analysis was restricted to the simulation environment with some major simplifications.For example, they used one-dimensional observations in terms of linear accelerations along the line-of-sight based on circular repeat orbits.Further, the GPS observation component as well as non-gravitational forces, such as atmospheric drag and solar radiation pressure, are missing.Therefore, investigating the application of ocean tide variance-covariance information in a more realistic simulation environment, especially in real GRACE/GRACE-FO Level-2 gravity field processing, is recommended and addressed in our study.
Next to the application of the variance-covariance information from ocean tides for monthly gravity field recovery, we again take up the idea of co-estimating ocean tides over a longer time period to improve ocean tide models for satellite gravity field recovery purposes (Kvas et al., 2019).This method is now applied in combination with the variance-covariance information for ocean tides.
Section 2 discusses the derivation of the ocean tide variance-covariance matrix (OTVCM).Section 3 illustrates the implementation of the tidal parameterization into the utilized software, followed by the description of the simulation and the real data environment with focus on differences to the current GFZ RL06.Section 4 presents the results as simulated and recovered gravity fields from real GRACE-FO data.Next to typical monthly gravity field retrievals, the potential of improving ocean tide models for satellite gravity field recovery by co-estimating ocean tide parameters when applying ocean tide uncertainty information, is addressed in Section 5.The conclusions and outlook are provided in Section 6.

Derivation of a Variance-Covariance Matrix for Ocean Tide Models
Ocean tide models have reached impressive accuracies by assimilating empirically estimated tidal constituents which achieve a near-global coverage due to satellite altimetry data.The root sum square error (RSSE) of the eight major ocean tides for the most capable models is below 1 cm in the open ocean, as measured by the complex ocean tide amplitude ζ with respect to a set of deep ocean bottom pressure recorders compiled by Ray (2013).This error level has to be benchmarked with the root-sum-square signal being 36.6 cm (Stammer et al., 2014).The models can therefore capture more than 97.5% of the open ocean tide signal.Nonetheless, the residual modeling errors induce considerable aliasing effects when processing satellite gravimetric solutions.These errors are most significant in polar regions, due to the inadequate coverage by satellite altimetry missions as well as the influence of sea ice, and in the shallow ocean, where the RSSE between data-constrained models is 4 cm for ζ, on average (Stammer et al., 2014).
We consider inter-model differences between five ocean tide models to derive OTVCMs for the eight major tides.We use FES2014 (Lyard et al., 2021), EOT11a (Savcenko & Bosch, 2012), TPX09 (Egbert & Erofeeva, 2002, updated), GOT4.10c (Ray, 2013), and DTU10 (Cheng & Andersen, 2011) excluding models that do not provide ocean tide information for the Antarctic shelf areas, as this would bias the error estimation.The necessary harmonization of the tidal models is described by Abrykosov et al. (2021).The ensemble of N = 5 ocean tide models is assumed to represent independent measurements of the respective partial tides, that is, varying in a Gaussian distribution around a specific mean value (cf. Figure 2b).The assumption of ocean tide models behaving like a Gaussian-distributed ensemble is debatable, as data sources for assimilation overlap between models.However, we argue that the created OTVCMs can nonetheless track the spatial error structure, as the variance between models in polar and coastal areas is mostly a product of processing strategies and hydrodynamical modeling.These areas connect to the better-constrained open ocean areas, which might be dominated by systematic errors due to satellite altimetry.However, external validation shows that these open ocean errors are minute.
A measure of uncertainty is given by the standard deviation with respect to the mean complex ocean tide amplitude  s , according to (1) Figure 2c presents the spatial error characteristics for the S 2 , highlighting the concentration of uncertainty in polar areas, coastal regions, and western boundary currents, which are all well-known regions are difficulty in data-constrained tide models (Hart-Davis et al., 2021).
As the native description of gravity functionals is by spherical harmonic functions, the five tidal models are converted to Stokes coefficients (see Abrykosov et al., 2021;Sulzbach et al., 2023) up to spherical harmonic (SH) degree and order (d/o) 30, which implies p = 961 base functions.Here, we employ real-valued SH functions Y {n,m} , of degree n and order m, where m < 0 denotes functions proportional to sin(mλ), m ≥ 0 identifies the respective functions proportional to cos(mλ) and λ denotes the geocentric longitude.This procedure is performed separately for the in-phase and quadrature components of the tidal constituents, which can be transformed to the prograde/retrograde notation (cf.Sulzbach et al., 2023, and references therein).The resolution was chosen in the SH domain as a compromise between the precision of approximating spatial errors and the computational demand and can be increased in future studies.
The OTVCMs describe the covariance of all 961 base function coefficients with each other, that is, they can be represented by a 961x961 matrix per partial tide.Since N < p, the empirical OTVCMs are not positive semidefinite and cannot be inverted.They would thus be useless for further applications (e.g., Hannart & Naveau, 2014;Ledoit & Wolf, 2003;Rajaratnam et al., 2008;Won et al., 2013).Even more, to obtain unbiased empirical OTVCMs the more stricter requirement N >> p would hold.To resolve this problem, it is thus necessary to regularize the estimated OTVCMs, for which we apply the Graphical Lasso (GLASSO) approach (Friedman et al., 2008), which additionally constrains the regularized OTVCM controlled by a regularization parameter Λ.With Λ increasing, the sparsity of the precision matrix (the inverse of the OTVCM) increases, that is, the number of connections in the underlying graph decreases.The regularization parameter was kept close to zero (Λ = 0.1), leading to rela tively small modifications with respect to the non-regularized, empirical VCM.While the choice of lambda is derived from qualitative arguments, tests with alternative parameters around 0.1 resulted in comparable results.The resulting correlation matrix r c (n 1 , m 1 , n 2 , m 2 ) is presented for S 2 in Figure 2a.
Hence, we allow for relatively large off-diagonal elements, considerably constraining the ocean tide solutions.It is crucial to keep in mind, that the obtained OTVCMs are constructed from a very small ensemble size of only five models, which will induce a bias that is hard to quantify.We refer to Abrykosov et al. (2021) for additional discussion.However, the chosen approach is based on reasonable assumptions and approximations, which will be supported by the results obtained when utilizing the OTVCMs (see following Sections).
The created matrices are symmetric and contain 462,241 independent elements, separately for in-phase and quadrature components of each tide.They are freely available under https://doi.org/10.5880/nerograv.2023.003and are employed in the following.

Methodology and Software
GFZ routinely processes monthly gravity field models from GRACE and GRACE-FO mission data as part of the GRACE/GRACE-FO SDS.The processing is performed by the Earth Parameter and Orbit System (EPOS) software package (https://www.gfz-potsdam.de/en/section/global-geomonitoring-and-gravity-field/topics/earth-system-parameters-and-orbit-dynamics/earth-parameter-and-orbit-system-software-epos/)applying the dynamical approach for satellite gravity field retrieval.A detailed description of the software and the latest processing strategy (RL06) can be found in Dahle et al. (2019).Simulated monthly gravity field solutions are derived using the same software and a processing scheme similar to GFZ RL06.Any differences between the simulation environment and real data processing are described in Sections 3.2 and 3.3.

Methodology
The method applied in EPOS for satellite gravity field recovery is based on the adjustment of SH coefficients.Therefore, the setup of the normal equation (NEQ) system is done with SH base functions of the first-order derivative of the Earth's disturbing potential ∂T/∂r (here given with respect to the satellite's radius vector r), also denoted as gravity disturbance, which can be expressed by a series expansion, according to Hofmann-Wellenhof and Moritz (2005).
In Equation 2, the gravitational constant is denoted by G, and M denotes the Earth's mass.The semi-major axis of the Earth is given by R and  Pnm is the normalized Legendre polynomial of SH degree n and order m.The fully normalized SH geopotential coefficients in terms of differences to a reference potential field are denoted by

𝐴𝐴
ΔCnm and  ΔSnm , and r is the mean Earth radius plus satellite altitude.Furthermore, the geocentric co-latitude and longitude are denoted by θ and λ.Equation 2 represents the disturbing accelerations acting on the satellite due to the disturbing potential T which is measured by GPS code and carrier phase observations, KBR observations (GRACE/GRACE-FO) and LRI observations (GRACE-FO).From Equation 2, we solve for  ΔCnm and  ΔSnm coefficients describing geophysical mass variations of Earth.The full gravity field potential coefficients are obtained by adding back the static reference gravity field as part of the linearization process of the observation equations.Apart from the gravity field coefficients, empirical accelerations, satellite state vectors and accelerometer calibration parameters are co-estimated.More information about these parameters and which parameters are set up during real data processing and simulations can be found in the corresponding following subsections.
Introducing uncertainty information from ocean tide models in terms of constrained NEQ systems for the eight major tidal constituents into the processing scheme of satellite gravity field recovery requires an extension of the parameter space for ocean tide parameters.The adaption to disturbing accelerations due to ocean tides can be accomplished by inserting the following equations into Equation 2 (Rieser et al., 2012)  ) . (4) Satellite gravity field recovery in EPOS applies the least squares adjustment method (Koch, 1987).According to this method, the functional model is usually written as and in terms of matrix notation where l denotes the observations from the satellites, A is the design matrix containing the partial derivatives ∂f(x)/∂x with respect to the unknown parameters  x to be estimated, and  v denotes the a posteriori residuals.Extending the parameter space to include ocean tide parameters requires the setup of partial derivatives ∂f(x)/∂x with respect to the unknown ocean tides.Thus, for every tidal constituent partial derivatives are formulated in terms of SH prograde and retrograde sine and cosine functions so that the vector of unknown parameters  x is extended by these ocean tide parameters.Analogously, the vector of observations l is extended by pseudo observations having zero values since the system is constrained with respect to the a priori ocean tide background model.In the sense of least squares, the a posteriori residuals  v should fulfill the minimum requirement  v  v → min in the L 2 -norm (Gauss-Markov model), with P as the weighting matrix for the observations.P can also be written as the inverse of the variance-covariance matrix (VCM) of the observations  ll . This matrix contains information regarding the uncertainty of observations.The variance-covariance information for a dedicated tidal constituent Q OT,s in terms of SH prograde and retrograde sine and cosine functions is given by We use a fully populated VCM for each of the eight major tidal constituents representing the uncertainty up to a maximum resolution of SH d/o = 30.Therefore, it cannot be excluded that correlations exist among parameters of specific tidal constituents.The uncertainty information for ocean tides determines basic conditions regarding the values of the tidal parameters to be adjusted.This enables the estimation of ocean tides at a reasonable level within short periods (e.g., monthly), despite the fact of very long-periodic aliases for K 1 and K 2 .
Following Equation 6, a best linear unbiased estimate of Before solving for the unknown parameters, we subdivide the total NEQ system regarding its individual observation components, hl-SST, ll-SST and ocean tides.For the latter, we introduce the variance-covariance information, according to Sulzbach et al. (2023), in terms of a constrained NEQ system.We adapted the original stand-alone OTVCM such that it includes already both, the left-hand side K = A T PA and the right-hand side k = A T Pl of the NEQ system.It contains the design matrix A, including the partial derivatives regarding the unknown ocean tide parameters, as well as the VCM, according to Equation 7. Implicitly, the right-hand side of the constrained ocean tide NEQ system is zero since the system is constrained with respect to the a priori ocean tide background model.In general, the total NEQ system can be split into two groups: the first one represents the left hand side K SST and the right hand side k SST for hl-SST and ll-SST components, as denoted in Equation 9, the second one defines the corresponding parts for the constrained ocean tide NEQ system K OT and k OT , as denoted in Equation 10.

𝟐𝟐
include the extended parameterization for ocean tide constituents, the remaining parts of K SST and k SST contain all other parameters, for example, the nominal gravity field parameters.K OT and k OT include only parameters regarding ocean tides.
The three NEQ system components (hl-SST, ll-SST and ocean tides) are summed up and the total system is solved, according to Equation 8.Among the estimated gravity field coefficients and other parameters, such as accelerometer calibration parameters (biases, scales) and satellite state vectors,  x contains estimates of the eight major tidal constituents  C ± nm,s and  S ± nm,s as well.Similar to the gravity field parameters, the ocean tide parameters are estimated in terms of a linearization process of the functional model.To obtain the full ocean tide signal, the a priori ocean tide model used for linearization purposes, is added back.

Simulation Environment
Due to the fact of knowing the true gravity field, the simulation world gives the opportunity to investigate the impact of dedicated error sources in the gravity field solution.In order to demonstrate the effectiveness of the OTVCM, we simulated monthly GRACE-FO-like gravity field solutions, including mass signals due to ocean tides, exclusively.The processing of the simulations within this study is based on the operational GRACE-FO RL06 method (Dahle et al., 2019), using some smaller adaptions, which will be addressed in this section.The simulated scenario is related to a typical GRACE/GRACE-FO near-polar orbit at 89° inclination, operating at an altitude of about 490 km.The initial orbit conditions for the trailing satellite are based on a real GRACE-FO orbit in January 2021 since the GRACE-FO mission had the desired orbit altitude at this time.The orbit of the leading satellite was generated by integrating the orbit of the trailing satellite up to an epoch where the desired inter-satellite separation of 200 km was reached.In the following, the simulation strategy including perturbation models and instrument noise assumptions, will be introduced in a short summary.
Forward modeling of the satellite state vectors in the simulation environment is performed by a precise numerical orbit integration (5-s integration step size) in terms of a multistep method using a two-times summed up Cowell approach of eighth order.The error-free hl-SST and ll-SST observations were computed from the satellite state vectors obtained from integration considering all force model-related quantities, summarized in Table 1 (forward modeling).The sum of the non-gravitational forces (air drag, solar radiation pressure and Earth albedo) was converted into error-free pseudo-ACC data in three directions.Drag compensation by simulating maneuvers was not considered during orbit integration.The Earth's gravitational attraction has been considered using background gravity field models.It includes a static gravity field model, a non-tidal time-varying gravity field model providing a six-hourly time series of potential changes from atmosphere (A), ocean (O), hydrosphere (H), Hz with L as the distance between two satellites in meter Amplitude spectral density for a frequency (f) band of 10 −4 to 10 −1 Hz for ATT noise (all directions)

Table 1
Force and Noise Models for Simulations cryosphere (I) and solid Earth (S) as well as an ocean tide model including the eight major tides.The simulation period was chosen as 1 year ( 2002), after which the satellite orbit decayed to 484.8 km (489.7 km after 1 month) when performing the orbit integration.The inter-satellite separation changes from 200 to about 176 km after one year (0.7 km after 1 month).This is comparable to changes observed for GRACE-FO, where the distance between the two satellites is usually kept within 200 ± 50 km by performing maintenance maneuvers.
The target signals of a satellite gravity field mission are geophysical signals caused by variations in TWS, mass loss in polar ice sheets and inland glacier systems.Consequently, gravity variations that do not represent these signals, such as ocean tides and non-tidal atmospheric and oceanic mass variations, need to be considered during the gravity field recovery process, that is, backward modeling.Therefore, during the recovery step, we substituted the nominal background models (time-varying tidal and non-tidal gravity field models) to simulate a realistic error level in the reference observations determining the time-varying gravity field (cf.Table 1, backward modeling).The static background model remains unchanged since errors due to the static gravity field do not contribute to the total error budget of the time-varying gravity field.To generate realistic hl-SST observations, we implemented the real GPS constellation of the year 2002, providing error-free code and phase observations with 40 cm white noise RMS for code and 0.3 cm white noise RMS for phase.ll-SST observations are computed geometrically by projecting orbit error-free positions and velocities onto the line-of-sight.The ll-SST component measures changes in distance or range-rates.In reality, it measures phase center to phase center for KBR, with a phase center offset correction applied to the measurements in order to refer them to the center-of-mass (CoM) of the corresponding satellite.LRI measurements have to be corrected in a similar way since the triple mirror assembly focal point does not coincide with the CoM exactly.These corrections are also subject to uncertainty but are neglected in the simulations, that is, we assumed already corrected ll-SST measurements, implicating that they refer to the CoM of each satellite.We generated noise time series for the key instruments, such as ACC, ATT and ll-SST, based on the noise models shown in Table 1 (backward modeling).The ACC noise time series were added to the pseudo-ACC observations in all three axes to obtain realistic non-gravitational measurements.For the ll-SST component, we assumed a GRACE-FO-like LRI, describing the in-orbit performance of the LRI sensor on GRACE-FO (Spero, 2021;Wegener et al., 2020) and added the LRI noise time series onto the line-of-sight observations.Uncertainty information for the satellite's attitude sensors was based on an approximation to the difference between SCA1B and LSM1B Level-1B star camera (SCA) and LRI steering mirror (LSM) products from the GRACE-FO mission (Goswami et al., 2021).For the difference between SCA1B and LSM1B, we formulated a noise model, which was assumed to be valid for all three directions (cf.Table 1, backward modeling).
The stochastic modeling of hl-SST GPS observations is based on unit matrices for code and phase observations, including variance information from pre-fit residuals estimated as a result of the computation of the linearized functions, including differences between the error-free and the reference observations.For the ll-SST component, we introduced a fully populated VCM derived from the covariance function of the sum of LRI noise and ACC noise in the time domain.The final weighting matrix for ll-SST is obtained by computing the inverse of the VCM using Cholesky decomposition routines.The stochastic model for ocean tides is introduced in terms of a constrained NEQ system, as described in Section 3.1.Since key instrument sensors are realistically modeled, we do not estimate empirical acceleration coefficients.
The NEQ system is setup arc-wise, where one arc is referred to as the length of 24 hr.Finally, the accumulated NEQ system is solved via least-squares adjustment, according to Equation 8.Besides SH coefficients of the Earth's gravity field (solved up to d/o = 120), orbital state vectors (position and velocity per arc and satellite) and ACC biases (three biases in radial and transverse direction and nine biases in normal direction per arc and satellite) have been adjusted.No ACC scale factors have been adjusted since the simulated pseudo-acceleration measurements do not contain scaled values.The omission of the ACC scale factors also implicates the argument of omitting ACC biases.Test simulations have shown that not estimating ACC biases does not change the quality of the gravity field solution.GPS receiver clock offsets and phase ambiguities have been pre-eliminated.If variance-covariance information of ocean tide model errors was applied during processing, ocean tide parameters for the eight major tidal constituents were resolved up to SH d/o = 30.

Real Data Environment
The processing of real GRACE-FO data in this study was performed using the EPOS environment and the latest release scheme, GFZ RL06.1, with one major exception: the stochastic modeling of ll-SST observations and the stochastic modeling of ocean tides by means of fully populated VCMs.GFZ RL06.1 differs from GFZ RL06 only by the inclusion of an updated SDS ACC transplant product for GF2.The updated transplant product defines a hybrid data version, where data from both satellites is used (McCullough et al., 2022).It has been shown that the application of the updated version of the GRACE-FO ACC transplant data during Level-2 gravity field processing leads to an improved gravity field solution for months being strongly affected by the satellite's beta prime angle, especially if the angle is close to zero degree (Dahle, 2022;Pie et al., 2022;Wiese et al., 2022).The beta prime angle defines the angle of the satellite's orbital plane with respect to the sun.An angle near zero degree is related to a time span in which the direction to the sun is in the orbital plane, that is, every half circle the satellites are in the Earth's shadow leading to a larger impact due to solar radiation pressure.Combined with the ACC transplant issue, this yields a degraded quality of gravity fields for particular months.
The stochastic modeling of observations in GFZ RL06.1 is realized through the derivation of root-mean-square (RMS) values from a priori or pre-fit residuals after completion of the a priori orbit generation (no gravity field parameters are setup), which is done to assure that convergence of a precise orbit determination (POD) using edited GPS, KBR range-rate (KRR) observations and iterated orbit parameter estimates is reached after one iteration (Dahle et al., 2019).The weighting of the GPS code-and carrier-phase observations as well as of the KRR observations, is achieved by unit matrices, including pre-fit RMS values, which are assumed to be representative of the uncertainties of GRACE/GRACE-FO key instruments.Further, empirical accelerations are setup in terms of sine and cosine parameters of once per revolution periodical models in along-track and cross-track directions.These parameters predominantly aim to reduce long wavelength gravity field retrieval errors, such as thermally-induced tone errors from the ACCs occurring at orbital frequencies (Kornfeld et al., 2019).However, the technique of empirical accelerations contains the danger of absorbing parts of the target gravity field signal.
Moreover, single pre-fit RMS values can usually not represent the realistic observation noise with correlations or frequency dependence.Therefore, the original GFZ RL06.1 processing scheme was modified.Similar to what was done in the simulation environment, a priori stochastic models have been developed for GRACE-FO KRR and laser interferometer range-rate (LRR) observations, following closely the approach of Murböck et al. (2023).
The technique described in this publication requires ACC data from instruments of both satellites.Therefore, the stochastic noise behavior of the ACC instrument was derived only for the GRACE mission due to the erroneous characteristics of the ACC on GF2 of the GRACE-FO mission.For our study, we used the GRACE ACC uncertainty information of the year 2014 (cf.Murböck et al., 2023, Figure 6).In this year, the ACC on GRACE showed degraded performance compared to previous years due to switched-off satellite thermal control.The noise behavior is assumed to represent the increased ACC noise behavior on GRACE-FO, caused by technical misconduct of the ACC sensor on GF2 shortly after the ACC sensor was put into science mode (Landerer et al., 2020).This assumption does not reflect the real ACC noise on GRACE-FO but at least provides a sufficiently good approximation to the error behavior of the ACC.The characterization and deduction of the stochastic model for the GRACE-FO KBR instrument noise has been accomplished equivalently to the method described in Murböck et al. (2023).The resulting stationary noise ASD model for GRACE-FO KBR observations can be approximated, according to for a frequency (f) band of 10 −4 to 10 −1 Hz.The combination of the ACC GRACE 2014 noise model with the KBR GRACE-FO noise model in range-rate domain, using a fully populated VCM, represents the stochastic model for Level-2 gravity field processing using GRACE-FO KRR 0.2 Hz data, provided by SDS.For LRI, we have used version 54 of GRACE-FO LRI Level1B 0.5 Hz datasets, provided by Albert-Einstein-Institute (AEI), Leibniz University of Hannover, Germany, publicly available at https://www.aei.mpg.de/grace-fo-ranging-datasets(L.Müller et al., 2022).The data include the correction due to the tilt-to-length (TTL) coupling (Wegener et al., 2020), which dominates the LRI error spectrum for frequencies between 1 and 35 mHz (GRACE/GRACE-FO gravity fields are typically solved for bandwidths between 0.2 and 20 mHz).Considering that the LRI noise model (cf.Table 1, backward modeling) is expected to be well below our ACC GRACE2014 noise model, it is justified to use only the latter.Thus, the weighting of GRACE-FO LRR observations during Level-2 processing has been realized by the ACC GRACE 2014 noise model, exclusively.The weighting scheme for GPS code and carrier-phase observations remains identical to GFZ RL06.
Finally, sets of monthly accumulated NEQ systems were generated for GPS, KRR and LRR observation components, each including additional parameterization for the ocean tides.In contrast to the ll-SST component, we do not apply realistic VCMs for the GPS component (as not yet been available for this study), leading to a certain imbalance between the NEQ systems of ll-SST and GPS.Therefore, on the NEQ level, an empirically derived relative weighting factor of 300 was applied to the KRR component in order to appropriately weight the KRR system relative to the GPS system.In the case of LRR, the relative weighting factor is about three orders of magnitude higher than for KRR due to higher sensitivity of the LRI and due to a sampling rate of 0.5 Hz instead of 0.2 Hz for KRR.Additionally, the OTVCM was scaled by a factor of 10 and completes the empirical relative weighting of the different NEQ system components.As it is done operationally, the total NEQ system of 1 month was resolved up to SH d/o = 96 for the target gravity field signal.The ocean tide parameters were resolved up to SH d/o = 30.

Results
The recovered gravity field solutions from simulations, as well as from real GRACE-FO data, have been assessed by analysis in the time-and frequency domain.Gravity signals and their errors are expressed in terms of equivalent water height (EWH).It denotes the height of a mass-equivalent column of water per unit area, based on the assumption that mass variations occur in a thin layer of the Earth's surface, approximated by the surface mass density.The relation between surface mass density coefficients and the gravitational potential coefficients is given by (Wahr et al., 1998) with  ave being the average density of the Earth (5,517 kg/m 3 ).In the frequency domain, we used amplitudes as quantity, computed per SH degree and per SH order Analogously, error information can be derived by the a posteriori formal standard deviations of the estimated gravity field coefficients  {  C nm ,  S nm } , also denoted as formal errors.These measures can be expressed in terms of EWH by applying Equation 12.In the time domain, we analyzed the gravity fields in terms of global spatial potential grids of 1° × 1° resolution by performing SH synthesis of the estimated potential coefficients with respect to the Earth's semi-major axis R. Mathematically, these grids describe the integration from the gravity disturbance into the disturbing potential T(R, θ, λ) by removing the factors -(n + 1)/r and (R/r) n+1 from Equation 2. As a further common metric to quantify recovered gravity fields over the ocean, we compute global latitude-dependent weighted RMS (wRMS) values over oceanic regions based on the global spatial grids by applying a continental mask with a 500 km buffer zone, according to We assessed the quality of the resolved gravity fields by considering all estimated coefficients of d/o ≥ 2, except for C 20 .The latter is badly determined from GRACE and GRACE-FO.Thus, it is recommended to replace this coefficient with values based on satellite laser ranging (SLR) measurements to higher orbiting geodetic satellites (Loomis et al., 2020).Additionally, we computed a posteriori residuals in order to assess the level of the fitted orbit as part of the linearization process with respect to the true orbit.The evaluation of the co-parameterized ocean tide parameters for every tidal constituent was done by computing the disturbing potential due to ocean tides T s (t i , R, θ, λ) (cf.Equations 2 and 3) for dedicated hourly epochs t i covering a period of 24 hr to consider the characteristics of semi-diurnal as well as diurnal tides.The global 1° × 1° grids were then transformed into half peak-to-peak amplitudes T peak, s , according to Equation 17, in order to incorporate the maximum and minimum characteristics of the ocean tides within 1 day.
The half peak-to-peak amplitudes were computed for the first of June in 2019.Since the characteristics of ocean tides are largely time-invariant, the computational epoch is irrelevant in this regard.
In the following, simulated gravity fields will be assessed first, in order to prove the effectiveness of the OTVCM.
Results already give a good indication of what can be expected in the real world.

Simulated Gravity Field Retrievals
In the simulation environment, we have analyzed a time series of 1 year (2002) of monthly recovered gravity fields in terms of gravity field retrieval error.Four sets of gravity fields have been processed: the first set includes the full spectrum of force-and noise models, which were introduced in Section 3.2, representing a realistic GRACE-FO-like scenario; the second one is identical to the first one but with uncertainty information of ocean tides (OTVCM) included in the gravity retrieval process; the third set only includes errors due to ocean tides (i.e., the non-tidal AO model, as well as the instruments, were regarded as error-free, no stochastic modeling of instruments was applied) without OTVCM included; the fourth set is identical to the third one but with uncertainty information of ocean tides applied in order to validate the effectiveness of the OTVCM.In the case of full signal gravity field recovery (sets 1 and 2), the retrieval error is obtained by removing the static part of the gravity field together with the averaged true HIS signal for the respective retrieved monthly field.For gravity fields which include only signals due to ocean tides (sets 3 and 4), the retrieval error is obtained by removing the static part alone since all time-variable non-tidal gravity field signals are the same in the forward and backward steps and thus removed in the retrieved gravity field solutions.All simulated results discussed in this study are presented as unfiltered, non-post-processed residuals.
The validity and improvement from the OTVCM are shown by the change from set 3 to set 4, illustrated in Figure 3.The graphic shows the ensemble of all 12 monthly retrievals (thin lines), together with the averaged degree amplitudes (bold lines).The true monthly HIS fields as well as the true ocean tide signal, represented through the EOT11a model, and the ocean tide error, represented through the difference between the EOT11a and FES2014 models, serve as reference signals.In the case of tidal retrieval error (sets 3 and 4), the signal amplitudes are dominated by temporal aliasing effects due to ocean tides, exclusively.The results demonstrate significantly reduced tidal aliasing errors when the OTVCM is applied.Despite the temporal aliasing effect induced by the observation geometry and the monthly retrieval period, the corresponding retrieval error (set 4, red line) even outperforms the error represented by the ocean tide model difference (brown line) in the longer wavelength spectrum up to degree 30.This indicates that the co-estimated tidal amplitudes can additionally improve the a-priori tide model.The latter is further investigated in Section 5.
The potential of the OTVCM is further demonstrated by computing the mean over the 1-year simulation period of the difference in magnitude for each estimated dimensionless gravity field coefficient between two gravity field sets (cf. Figure 4d).Similarly, Figure 4b depicts the mean values computed in the spatial domain.Both metrics reveal significantly reduced tidal retrieval errors over the whole spectrum due to the OTVCM and the co-parameterization of ocean tide parameters.It is assumed that a large part of the erroneous ocean tide signal is captured by the co-estimated tidal constituents, which is demonstrated in the scope of the evaluation of the real GRACE-FO gravity field retrievals in Section 4.2 (cf. Figure 10).The aliasing error characteristics due to ocean tides visible in Figure 4b also arise above the continents.This is caused by the along-track observation geometry of the satellites leading to smearing effects transmitted to the continents due to tidal aliasing errors close to the coasts.The averaged ocean wRMS based on tidal retrieval with and without OTVCM (cf. Figure 4e) is 2.95 cm EWH and 0.75 cm EWH, respectively, resulting in a reduction of tidal aliasing errors of about 75%.
In contrast to the solutions including only ocean tide signal, results representing the full retrieval error show a different picture.While the spectral analysis in terms of degree amplitudes shows hardly any differences, the view on the coefficient triangle (cf. Figure 4c) as well as on the global grid (cf. Figure 4a) depicts differences between retrieval errors at specific bandwidths of the time-varying gravity field due to the application of the OTVCM.
Clearly, visible reduced aliasing errors arise for the prominent resonance order bands (m = 15, 30, 45, etc.) and neighboring orders resulting in a reduced error pattern for certain areas at higher latitudes in the spatial domain, such as around Antarctica and in the Arctic ocean.This was proven by performing the SH synthesis only for coefficients, including the prominent resonance orders.Furthermore, significantly reduced error characteristics are seen east of Patagonia.This region is known for an artifact in tidal induced oceanic mass flux changes due to eddy phenomena (Tonini & Palma, 2017).The regions mentioned are typically known for largest errors in ocean tide models (cf. Figure 1) (Knudsen & Andersen, 2002;Ray & Luthcke, 2006).For results including the full retrieval error by considering non-tidal AOHIS mass variations, the retrieval error is dominated by temporal aliasing errors due to the mismodeling effects of the AO de-aliasing model.Therefore, the beneficial effect of the OTVCM is hidden behind the remaining dominating error signals and becomes visible only in a limited sense.Consequently, ocean wRMS values related to results of full retrieval errors show, on average, only small differences of 9.53 cm EWH and 9.22 cm EWH, respectively, and finally results in a reduction of temporal aliasing errors of about 3% when applying the OTVCM.

Gravity Field Retrievals From GRACE-FO
The assessment of the effectiveness of the OTVCM applied to real data was analyzed by means of KRR and LRR GRACE-FO observations of the year 2019.February and March solutions based on LRR observations could not be processed due to missing data caused by diagnostic LRI scan operations and problems with the on-board computer.In the case of the KRR February solution, we oriented to the official SDS solution span and combined the last 5 days of January and the first 5 days of March together with the only 14 available days of February.All gravity field retrievals were processed using the updated ACC transplant product for GF2.Results obtained according to the GFZ RL06.1 scheme serve as reference solutions.We would like to mention that GFZ RL06.1 includes only KRR observations, which is to be considered in the following discussions and illustrations.Besides, Level-2 products were generated, including realistic stochastic modeling of instruments without the application of the OTVCM, hereinafter referred to as V1, as well as products including both, realistic stochastic modeling of instruments and the introduction of uncertainties of ocean tide models, hereinafter referred to as V2.The error assessment was performed in terms of formal errors as well as in terms of residuals with respect to a GRACE/GRACE-FO climatology model (Murböck et al., 2023, see Equation 4).This model was derived by estimating six parameters (bias, linear trend and sine/cosine amplitudes of the annual and semi-annual signal) for each SH coefficient from the whole unfiltered International Combination Service for Time-variable Gravity Field (COST-G) GRACE and GRACE-FO time series (Meyer et al., 2020)   The residual SH degree and order amplitudes (cf. Figure 5) demonstrate significantly reduced error signals for KRR/LRR V1 and V2 residuals, compared to GFZ RL06.1.The differences in magnitude are somewhat more pronounced in the order amplitudes due to the order-dependent vertical structure of the error signal improved by the OTVCM (cf.Figures 6a-6d), while the degree amplitudes smooth this out across all orders.The application of a more realistic stochastic model for instruments leads to clearly reduced errors for the short wavelength spectrum (above d/o = 70, cf. Figure 5a), which is also confirmed in Murböck et al. (2023).Additionally, gravity field retrievals benefit through the introduction of the OTVCM for the mid-to-high degree bandwidth (between d/o = 30 and 70).We assume that temporal aliasing errors due to ocean tides are reduced in the long wavelength spectrum as well, as it has been demonstrated by the simulations (cf. Figure 3, red and blue lines, Figures 4b  and 4d), but the benefit is limited by the dominating error sources from AO and the ACC.The corresponding coefficient triangles (cf.Figures 6c and 6d) confirm this fact by showing the mean of absolute differences of dimensionless SH coefficients among V1 and V2 solutions.They show that the benefit of the OTVCM is mainly related to the prominent resonance order bands (15, 30, 45, etc.) and neighboring orders.
Further, the triangles reveal larger differences between V1 and V2 at the high frequency spectrum (above d/o = 70), showing increased as well as reduced error amplitudes for specific coefficients.Another aspect is the slight degradations for specific coefficients of the very long wavelength spectrum (very upper part of the triangles).They are somewhat more pronounced for the differences between GFZ RL06.1 and V2, compared to V1 and V2, due to the parameterization of empirical accelerations.The main differences between V1 and V2 as well as between GFZ RL06.1 and V2 arise from the relative weighting of the different observation and pseudo-observation components (hl-SST, ll-SST and OTVCM).In this study, we applied identical weighting factors for every month and for different measurement groups, respectively, as introduced in Section 3.3.Tests using slightly larger or smaller weighting factors (equivalent to the up-or down-weighting) have shown that the resolved gravity field coefficients are very sensitive to the applied weighting factors.Therefore, the applied weighting factors might not lead to the best possible gravity field recovery for specific months, depending on   the number of observations available and their quality.However, the factors introduced in Section 3.3 define a first compromise between the different observation components in order to obtain gravity field solutions of high quality.The application of a method for relative weighting, such as variance component estimation (VCE), can help to find more optimized weighting factors for the respective observation components.The implementation of VCE is currently being investigated in the scope of a future GFZ RL07.
The formal errors represent uncertainty estimates for the retrieved SH coefficients.Thereby, they include the pure effect due to the sensitivity of gravity field observations together with the underlying observation geometry and do not include any information from time-varying geophysical signals.Thus, they are directly dependent on the stochastic model applied.The formal errors of KRR/LRR V1 and V2 results, depicted in Figures 5a and 5b (red and green dashed lines) and in Figures 6e-6h, show a more realistic behavior compared to GFZ RL06.1 (cf.Figures 5a and 5b magenta dashed lines and Figure 6i) over the complete spectrum due to the consideration of variance-covariance information for instrument noise and ocean tide background model errors.The expression "more realistic" is related to larger signal amplitudes and a better approximation of the formal errors toward the characteristics of the residuals showing smoother amplitudes around SH 15th-and 30th -order resonances, for instance.The coefficient triangles related to formal errors of V1 and V2 solutions reveal higher amplitudes around the prominent resonance order bands and neighboring orders (Murböck et al., 2023).Characteristics of formal errors due to the application of the OTVCM are related to the upper part of the respective triangles (up to SH d/o = 30), showing slightly increased signal amplitudes compared to V1.This effect is also reflected in the very low degree spectrum of the degree amplitudes (cf.Figures 5a and 5b).We would like to emphasize that a more realistic stochastic modeling of GPS code-and phase observations is not yet included and therefore, the formal errors might still appear too optimistic.
Aside from this, the availability of both, KRR and LRR observations, enables the comparison of their respective residuals.Our results demonstrate high consistency among KRR-and LRR-based retrievals with benefits for LRR-derived residuals related to zonal (m = 0) and near-zonal coefficients of high degrees due to the increased measurement accuracy of the LRI instrument.This fact is illustrated by the comparison of Figures 5c and 5d (lower order spectrum) as well as Figures 6a and 6b, Figures 6e and 6f and Figures 6g and 6h.However, the full potential of the LRI cannot be exploited due to the dominating error sources originating from temporal aliasing errors and ACC errors.
For the assessment in the spatial domain (cf. Figure 7), we used post-processed residuals using a Gaussian filter with a 300 km radius, which were split into two bandwidths, the full wavelength spectrum from SH d/o = 2 until 90 and the short wavelength spectrum from SH d/o = 30 onwards.While the latter includes gravity field signals mainly affected by temporal aliasing effects as well as instrument noise, the full wavelength spectrum still includes, besides aliasing errors and ACC errors, time-varying mass change signals, which have not been removed by the climatology.These signals are mainly related to hydrological and ice mass fluxes, which is the part of the signal that scientists are interested in.Similar to the coefficient triangles, the spatial representations were computed as the mean of differences in magnitude of respective global grids representing GFZ RL06.1, V1, and V2 residuals.Grids depicting the short wavelength spectrum (cf.Figures 7e-7h) clearly demonstrate reduced aliasing errors at low-to-mid latitude areas for V2 residuals due to the consideration of ocean tide model inconsistencies by the OTVCM and the co-estimation of tidal parameters.Increased error amplitudes at the very high-frequency spectrum in the coefficient triangles (cf.Figures 6a-6d) are mirrored as slightly increased error signals in the spatial grids (cf.Figures 7e-7h).However, V2 results outperform those of V1 and demonstrate greatly reduced aliasing errors compared to GFZ RL06.1.Grids representing the full spectrum (cf.Figures 7a-7d) are dominated by mass signals that have not been removed by the climatology.
Consequently, the benefit of the OTVCM is less pronounced here.The impact of the reduced tidal aliasing errors in the shorter wavelength spectrum discussed before is also reflected in the full spectrum grids, mostly in low-to-mid latitude areas.However, reduced aliasing error effects are also visible in higher latitude areas, especially around the Antarctic Peninsula.The spatial grids also reveal some increased error signals at very high latitudes (cf.Figures 7c  and 7d) which are related to the very low degree and order coefficients (cf.Figures 6c and 6d, very upper parts).
The estimation performance of these coefficients is largely dependent on the relative weighting of the individual observation components when solving for the NEQ system and needs to be further investigated in the future.
Similar to the simulation results, the quantitative assessment of GRACE-FO recovered gravity field models was performed by means of ocean wRMS values (cf. Figure 8).They demonstrate, on average, significantly smaller values for KRR V2 residuals of about 23.3%, and about 22.0% for LRR V2 residuals, compared to GFZ RL06.1.The comparison between V1 and V2 residuals shows an average reduction of about 6.6% for KRR V2 and 5.3% for LRR V2.The results indicate that months including fewer observations, such as February, benefit from the consideration The introduction of ocean tide uncertainty information into Level-2 gravity field processing is a further step toward an optimized gravity field processing.This is also demonstrated by the a posteriori, or post-fit residuals   v (cf.Equation 16).The representation in the spatial domain (cf. Figure 9) depicts the KRR and LRR post-fit residuals for GFZ RL06.1, V1 and V2 along the respective ground tracks in terms of bandpass-filtered accelerations for the example of June 2019.The illustrations show mainly sub-monthly hydrological signals as well as temporal aliasing errors and instrument noise.The comparison between V1 and V2 residuals demonstrates reduced error signals due to the application of the OTVCM, especially around Antarctica and in the Arctic ocean.Due to the fact that co-estimated ocean tide parameters include most of the erroneous tidal mass signals, as will be shown in the following discussion, the amplitudes of the V2 post-fit residuals are smaller.The higher sensitivity of the LRI compared to KBR is depicted by post-fit residuals which are much less affected by instrument noise.This is also reflected by their respective standard deviations in terms of nm/s 2 , reading 0.82 (RL06.1),0.81 (V1 KRR), 0.74 (V2 KRR), 0.48 (V1 LRR), and 0.37 (V2 LRR).
Figure 10 depicts the spatial characteristics of the monthly co-estimated tidal parameters as the sum of the eight major tidal constituents in terms of residuals with respect to the reference model FES2014 for the example of V2 KRR retrievals.The illustrations show mainly erroneous signals due to the imperfection of the applied ocean tide background model appearing at locations where these are to be expected (cf. Figure 1), such as around Antarctica, in the Arctic ocean and at dedicated coastal regions which are known for larger tidal error signals.Residual ocean tide parameters for February show larger amplitudes originating from generally larger temporal aliasing errors due to the 2 weeks of missing observations.However, the results also reveal some smaller signals over the continents at specific months, such as in the Amazon region and in central Africa.This issue is related to the contribution of the S 2 tidal component, whose aliasing period (∼161 days) correlates with the semi-annual hydrological signal.This correlation can be explained by hydrological signals erroneously entering into the adjusted parameters of the S 2 tide.The analysis of the individual monthly V2 KRR/LRR gravity field solutions (not shown here) indicates no notable artifacts in regions where semi-annual hydrological signals appear, compared to V1 fields.Nonetheless, this effect should be kept as small as possible since hydrological mass anomalies define one of the major target signals to be recovered.

Comparison of Simulated and GRACE-FO Real Data Gravity Field Retrievals
Contrasting real data solutions with simulated products is always a tightrope act since the truth is perfectly known in the simulation environment and, for example, remaining mass change signals not being removed correctly by the climatology do not exist in the simulations.This issue is the major reason for the differences when comparing simulated full retrieval errors, including OTVCM (set 2) with V2 residuals.Furthermore, gravity field observations from GRACE-FO include errors from the ACC transplant data, which are reflected in the adjusted gravity fields as error signals but do not occur in the simulated solutions.The simulations depict the largest aliasing error reduction in the Antarctic and Arctic regions, which are related to decreased resonance order effects around SH orders of 15, 30 and 45, and some smaller error reductions in lower latitude areas, which are related to less aliasing errors at the higher frequency spectrum (cf.Figures 4a and 4c).The real data assessment demonstrates the largest improvements at the mid-to-high frequency bandwidth (cf.Figures 5a and 5b, Figures 6c and 6d  and Figures 7c, 7d, 7g, and 7h), which are related to mid-to-low latitude areas in the spatial domain.Compared to the simulations, the larger error reduction of ocean tides induced temporal aliasing effects at the mid-to-high frequency spectrum in real data solutions is caused by a larger amount of temporal aliasing error signals, which are likely underestimated in the simulation world.Similar characteristics of error reduction exist in real data solutions as well as simulated solutions, such as less pronounced resonance order effects.
Further similarities are seen around the Antarctic Peninsula in terms of a reduced error pattern.The simulations have shown that temporal aliasing errors due to ocean tides are reduced in the long wavelength spectrum as well when excluding instrument noise and non-tidal AO signals (cf. Figure 3, red and blue lines and Figures 4b  and 4d).Therefore, we assume that tidal aliasing reductions also exist at the longer wavelength spectrum in simulation results, including the complete noise and force model ensemble, and in real data solutions.Especially for the latter, the long wavelength part is dominated by errors due to the ACC transplant and due to non-tidal AO signals, which does not allow the influence of the OTVCM to become visible.

Determination of an Alternative Ocean Tide Model
In Section 1, the possibility of improving ocean tide models for satellite gravity field recovery by means of multiyear GRACE/GRACE-FO gravity field observations has already been addressed (Kvas et al., 2019).Following this idea, we estimated ocean tide parameters for the eight major tidal constituents over a time span of only 1 year (as a first test scenario) using the OTVCM.Due to the well-controlled simulation environment, we investigated this method by means of simulations, as introduced in Section 3.2, with the difference that the target signals were changed to mass variations due to ocean tides instead of non-tidal time-variable hydrological and ice redistributions.As it was already the case in the previous simulated gravity field recovery, we again assumed the true ocean tide signal to be the EOT11a model as part of the ensemble of background models needed for forward modeling, and the FES2014 model has been considered as the reference for the backward modeling.All processing steps remained identical as for the nominal monthly gravity field recovery in the simulation environment.To prevent time-varying monthly gravity field signals from entering into the ocean tide parameters, the daily NEQ systems were accumulated to monthly NEQs in the first step, followed by the pre-elimination of the nominal gravity field parameters.In the second step, the monthly NEQs were then accumulated to the full period of 1 year.Finally, the total system was solved for the tidal constituents as well as for arc-specific parameters and ACC biases.Since the estimated ocean tide parameters are part of the linearization process, the reference model FES2014 was added back after solving for the unknown parameters.The resulting total tidal signal defines a new ocean tide model, which can be directly validated against the true signal from EOT11a.
Figure 11 illustrates the results in the spatial domain.In order to get a comprehensive overview of the characteristics of ocean tides, we show the total signal amplitudes from EOT11a of the eight major tidal constituents, respectively, together with the two sets of residual signal amplitudes, which are EOT11a minus FES2014 serving as initial ocean tide error, and EOT11a minus the estimated ocean tide parameters as the final ocean tide error.The latter represents the ocean tide model error due to the co-estimation of the tidal constituents from satellite gravity field observations.The results demonstrate reduced error amplitudes compared to the initial ocean tide errors for all estimated tides at locations known for larger uncertainties of a specific tide.This can be achieved by global gravity field observations, including very high-latitude areas and, by combination with the OTVCM, which prescribes less of a constraint in those areas due to a larger spread in the input model ensemble.Error characteristics of diurnal tides arise mainly at high-latitude areas around Antarctica (especially in the Filchner-Ronne Ice Shelf along the Antarctic Peninsula) and in the Arctic ocean (especially west of Greenland), while erroneous signals of semi-diurnal tides also appear at low-latitude coastal areas (especially along the northern coast of Australia).We can show that tidal estimates based on the OTVCM reduces predominantly the strongly pronounced error amplitudes originating from tidal constituents having the largest signal amplitudes (M 2 , S 2 , O 1 , and K 1 ).However, errors from weaker tidal constituents, which contribute less to the total ocean tide error budget, can be reduced with this method.In reality, the estimation of ocean tides with observations from the GRACE and GRACE-FO missions requires a sufficiently long observation period (ideally the whole GRACE/GRACE-FO time series) due to long-periodic aliases for K 1 and K 2 constituents.
Simulation results show that the application of the OTVCM allows for the reasonable estimation of all eight major tidal constituents already after 1 year, showing reduced error amplitudes.This finding even holds true for K 2 , which has a long alias period and a small signal amplitude which makes it challenging to de-correlate the tidal signal from other tidal constituents and non-tidal time-varying gravity field parameters.However, the error characteristics in S 2 show signatures over the continents (Amazon region, central Africa and northern Siberia), which couple into estimated parameters of the S 2 tide due to the correlation with the semi-annual hydrological signal.A longer observation period of several years and a different individual weighting of the S 2 component of the OTVCM can help mitigate this effect, which needs further investigation.A quantitative overview of the tidal estimates with respect to the reference ocean tide error in terms of absolute values is given in Table 2.In total, tidal error signals could be reduced by about 15%.Further reductions are expected when extending the observation period up to several years.We emphasize that the estimated tides can be used as an alternative ocean tide background model, embedded in the usual ocean tide de-aliasing model, for Level-2 gravity field recovery.The method does not aim to improve current ocean tide models at small-scale coastal areas.This would require much higher spatial resolution, which, in turn, would result in an increasingly unstable and extremely large NEQ system.In the real world, the general idea would be to exploit the whole GRACE/GRACE-FO time series and reprocess it in a first step with the goal of refining ocean tides using discussed updated parameterization methods toward GFZ RL07.Such an alternative ocean tide model can be merged with an actual tide model, which is used as a background field for de-aliasing purposes and is potentially optimized for GRACE/GRACE-FO corrections.

Conclusions and Outlook
In this study, we have investigated the potential of introducing uncertainty information from ocean tide models for Level-2 gravity field recovery from satellite gravimetric observations by means of realistic simulated data as well as real KBR and LRI range-rate GRACE-FO observations.The ocean tide error information has been derived in the SH domain from inter-model differences between five different ocean tide models in terms of regularized empirical variance-covariance information for the eight major tidal constituents resolved up to SH d/o = 30.We have introduced this uncertainty information in terms of a constrained NEQ system, together with realistic stochastically modeled instrument noise from the key instruments (ACC and ll-SST), into the functional model of GFZ's EPOS software by extending the parameter space for additional tidal parameters to be adjusted.
In order to demonstrate the effectiveness of the OTVCM, we have simulated monthly batches of GRACE-FO-like gravity field retrievals over 1 year, including mass signals due to ocean tides, exclusively.We have shown that the gravity field retrieval error caused by tidal aliasing effects is reduced by about 75% in terms of ocean wRMS when applying the OTVCM and co-estimating tidal constituents.The largest improvements are seen around Antarctica as well as in the Arctic Ocean which are typical regions known for the largest error signatures in ocean tide models.The reduced error signals in these regions are reflected in reduced error amplitudes for resolved gravity field coefficients of the prominent resonance order bands (15, 30, 45, etc.) and neighboring orders.In contrast, results including the full error spectrum, that is, temporal aliasing effects due to non-tidal time-varying signals and instrument noise, show reductions in the retrieval errors at a smaller magnitude (∼3% ocean wRMS) due to the dominating error sources from mismodeled non-tidal atmospheric and oceanic mass distributions and the ACC.However, simulations including the full ensemble of disturbing forces still demonstrate reduced error signatures in the regions mentioned before and even for the well-known eddy anomaly east of Patagonia in the south Atlantic Ocean.
In the scope of real Level-2 data processing, we have assessed the impact of the OTVCM by recovering monthly gravity fields from KBR and LRI GRACE-FO range-rate measures for the year 2019.We have shown that this signal is mostly related to erroneous signatures due to imperfect ocean tide models arising at typical locations on Earth, for example, polar regions.This proves that the application of the OTVCM during Level-2 gravity field processing is capable of absorbing a large amount of signal which is usually found in gravity field solutions in terms of temporal aliasing errors.Results have shown that the error in terms of ocean wRMS can be reduced on average by about 5%-7% for LRR and KRR, compared to retrievals without the OTVCM applied and can be reduced by up to more than 20% with respect to GFZ RL06.1.However, the assessment based on synthetic observations as well as on real GRACE-FO measures has revealed that the visibility of the beneficial effect of the OTVCM at the level of recovered gravity fields is limited since non-tidal high-frequency mass variations from atmosphere and ocean as well as errors from the ACC still dominate the error spectrum of the time-varying gravity field.
Furthermore, we have estimated tidal parameters including the uncertainty information for ocean tide modeling using satellite gravimetric observations in order to improve tide models for Level-2 gravity field processing.Simulated results have demonstrated that error amplitudes of the eight major tidal constituents could be reduced by about 15% when estimating ocean tides only for a period of 1 year, despite the fact of very long-periodic aliases of several years for K 1 and K 2 constituents.Even the errors of tidal constituents, which have small signal amplitudes making estimation from shorter time series challenging, such as Q 1 , P 1 , N 2 , and K 2 , and which makes it difficult to decorrelate them from other tidal constituents as well as from non-tidal mass variations, can be reduced.The estimation of S 2 shows some artifacts at specific continental regions due to the correlation among the aliasing period of S 2 and the semi-annual hydrological signal.This issue needs to be further investigated in future.
The incorporation of variance-covariance information for ocean tides into Level-2 satellite gravity field processing constitute one step further toward GFZ RL07.Similar to what has been done with ocean tides, we plan to stochastically model the errors from high-frequency non-tidal atmospheric and oceanic mass variations (more specifically, the AOD Level-1B de-aliasing error product, Dobslaw et al., 2017) in order to reduce the dominating impact of these errors.This opens the door for further reduction of tidal aliasing errors by means of the OTVCM.Additionally, we are working on the stochastic modeling of GPS code and phase observations which complete the realistic modeling of uncertainty information for all different GRACE/GRACE-FO observations and enables the proper execution of the VCE.The latter will help to find more optimal relative weighting factors for the different observation-and pseudo observation components.The combination of these changes in parameterization and stochastic modeling in a future RL07 might help to further improve GRACE/GRACE-FO gravity field recovery.

Data Availability Statement
The regularized empirical variance-covariance matrices for stochastic modelling of the eight major ocean tides are freely available under https://doi.org/10.5880/nerograv.2023.003.
to the initial ocean tide model to define an updated ocean tide de-aliasing model showing improvements in areas, such as at the Poles or in the Indonesian Sea, where the GRACE and GRACE-FO gravity solutions often show deficiencies.

Figure 1 .
Figure 1.RSS of standard deviations from mutual differences between the ensemble of five different ocean tide models with respect to the mean including the eight major tidal constituents in centimeters of complex ocean tide amplitude ζ.

Figure 2 .
Figure 2. Correlation matrix r c (n 1 , m 1 , n 2 , m 2 ) of the regularized OTVCM up to SH d/o = 4 for the partial tide S 2 , which corresponds to correlations between SH coefficients   {1,1} and   {2,2} (a).The SH functions are ordered from low to high degree n.(b) shows the mean complex amplitude over the five different ocean tide models for the partial tide S 2 , and (c)shows the standard deviation with respect to the complex amplitude for the partial tide S 2 in the spatial domain.

Figure 3 .
Figure3.SH degree amplitudes of the monthly residuals based on simulations including the full spectrum of force and noise models (orange and green lines) and simulations including only ocean tide signals (blue and red lines).Thin lines are the monthly residual amplitudes, bold lines are the averaged residual amplitudes.Simulations with variance-covariance of ocean tides applied are labeled with OTVCM.As references, the monthly averaged HIS amplitudes are displayed (black lines), together with the total ocean tide signal (purple line) and the differential ocean tide signal (brown line).The latter were computed over a time period of 24 hr for the first of June in 2019, respectively.Amplitudes are displayed in terms of cm EWH.
from 2002 to the year 2020.The COST-G Level-2 products can be regarded as high-quality gravity field solutions derived by combining the solutions of several Level-2 processing centers.Therefore, we assumed this climatology model to be representative of the climatological signal, despite the fact that the true climatological signal is not known.Since COST-G products and the resulting climatology model are available only up do SH d/o = 90, the assessment of gravity field retrievals was done up to this maximum resolution.

Figure 4 .
Figure 4. Mean of differences in magnitude of full retrieval errors based on simulations with and without OTVCM applied (a), (c), and tidal retrieval errors based on simulations with and without OTVCM applied (b), (d).Results are displayed in terms of global grids resolved up to SH d/o = 120 in cm EWH (left) and in terms of dimensionless coefficients (right).RMS weighted with the cosine of latitude over the ocean using a 500 km buffer zone of monthly global grids in terms of retrieval errors resolved up to SH d/o = 50 in cm EWH (e).

Figure 5 .
Figure 5. Residual amplitudes computed per SH degree (a), (b), and per SH order (c), (d) of monthly recovered gravity field models with respect to a GRACE/GRACE-FO climatology in cm EWH.Thin lines are related to the monthly retrievals, bold lines represent the averaged degree/order amplitudes.For reasons of clarity, the formal errors (dashed lines) are displayed only in terms of averaged degree amplitudes.

Figure 6 .
Figure 6.Mean of differences in magnitude of dimensionless residual coefficients between GFZ RL06.1 and V2 solutions (a-b), between V1 and V2 solutions (c-d), and mean of a posteriori formal standard deviations of estimated coefficients of V2 (e-f), V1 (g-h), and GFZ RL06.1 (i) solutions.Left column is related to KRR, right column is related to LRR.

Figure 7 .
Figure 7. Mean of differences in magnitude of 300 km Gaussian filtered residual global grids in cm EWH for the full spectrum (a-d) and the spectrum > SH d/o = 30 (e-h).Left column is related to KRR, right column is related to LRR.Information about which solution version is depicted can be found in the sub-labels.

Figure 8 .
Figure 8. Ocean wRMS using a 500 km buffer zone computed from 300 km Gaussian filtered monthly residual global grids in terms of cm EWH.

Figure 10 .
Figure 10.Residual ocean tide parameters for the sum of the eight major tidal constituents of monthly recovered V2 KRR gravity field models with respect to the reference ocean tide model FES2014 in terms of half peak-to-peak amplitudes over 24 hr in cm EWH computed for the first of June 2019.The monthly grids are ordered row-wise starting with January 2019 at the top left and ending with December 2019 at the bottom right.

Figure 11 .
Figure 11.Half peak-to-peak amplitudes over 24 hr of diurnal tides (a), semi-diurnal tides (b) and the sum of the eight major tidal constituents (c) in cm EWH computed for the first of June 2019.Amplitudes are computed in terms of the total signal (first row of each tidal constituent group and left plot in the accumulated group), the difference between the EOT11a and FES2014 models (second row of each tidal constituent group and middle plot in the accumulated group) and the difference between the EOT11a model and the estimated ocean tide parameters (third row of each tidal constituent group and right plot in the accumulated group).
Gravity field retrievals related to both KBR and LRI observations show consistent results with an improved retrieval performance for coefficients corresponding to the mid-to-high frequency bandwidth (∼from SH d/o = 30 up to d/o = 70) of the time-varying gravity field.Similar to the results of the synthetic world, coefficients corresponding to the prominent resonance order bands experience the largest reduction of aliasing errors.In the spatial domain, reduced error signatures are seen, especially around the Antarctica Peninsula and at low latitude areas in the high-frequency spectrum.Additionally, the consideration of ocean tide variance-covariance information leads to more realistic formal errors in the long wavelength spectrum up to SH d/o = 30.The spatial analysis of post-fit residuals, which are based on the retrieval with OTVCM applied, demonstrate reduced error amplitudes compared to residuals which are based on the retrieval without using the OTVCM.The smaller error characteristics are especially pronounced around Antarctica and in the Arctic Ocean resulting in reduced standard deviations of about 9% for KRR and 23% for LRR.The extension of the parameter space for ocean tides implicates the fact that the same signals, which usually enter into the nominal gravity field parameters or post-fit residuals, are incorporated into the co-estimated tidal constituents.

Table 2
Global wRMS in cm EWH for the Eight Major Tidal Constituents and the Sum of These Derived From Half Peak-To-Peak Amplitudes Over 24 hr Computed for the First of June 2019