Earth’s Time-Variable Gravity from GRACE Follow-On K-Band Range-Rates and Pseudo-Observed Orbits

: During its science phase from 2002–2017, the low-low satellite-to-satellite tracking mission Gravity Field Recovery And Climate Experiment (GRACE) provided an insight into Earth’s time-variable gravity (TVG). The unprecedented quality of gravity ﬁeld solutions the gravity ﬁeld recovery processing strategy applied compare the obtained ﬁeld results to the solutions of compare the GRACE-FO performance to that of the preceding GRACE mission in terms of post-ﬁt residuals. use the in-house-developed MATLAB-based GRACE-SIGMA software compute unconstrained solutions based on the generalized orbit determination of 3 h arcs. K-band range-rates (KBRR) and kinematic orbits are used as A comparison of the obtained solutions to the results of the GRACE-FO Science Data System (SDS) and Combination Service for Time-variable Gravity Fields (COST-G) ACs, reveals a competitive quality of our solutions. spatial noise solutions out comparison GRACE-FO post-ﬁt an improvement of the GRACE-FO K-band ranging system performance. overall GRACE-FO solutions from June 2018 until December 2020. The regularly updated LUH-GRACE-FO-2020 time series of monthly gravity ﬁeld solutions can be found on the website of the International Centre for Global Earth Models (ICGEM) and in LUH’s research data repository. These operationally published products complement the time series of the already established ACs and allow for a continuous and independent assessment of mass changes in Earth’s system.


Introduction
The two identical satellites of the twin satellite mission Gravity Recovery and Climate Experiment (GRACE) [1] were orbiting the Earth from March 2002 until December 2017/March 2018 (GRACE-B/GRACE-A) in a near-circular and near-polar low Earth orbit Association of Geodesy (IAG). Within COST-G, monthly gravity field solutions of participating ACs and partner ACs are combined on solution [46] and normal equation levels [47], in order to provide consolidated products with improved quality and higher robustness. Since January 2021, LUH/IfE has been an official COST-G GRACE-FO AC and is contributing with the LUH-GRACE-FO-2020 solutions to the COST-G operational GRACE-FO time series of monthly gravity fields [48].
Here, we present and discuss the processing strategy of the LUH-GRACE-FO-2020 time series of monthly gravity field solutions. We give an overview on the underlying theory (see Section 2.1), as well as review and summarize the corresponding mathematical framework in a compact form (see Appendices A.1 and A.2).The processing strategy of the LUH-GRACE-FO-2020 time series is summarized in Section 2.2. Utilized GRACE-FO data products are listed in Section 2.3. In Section 3, the obtained monthly gravity field solutions are evaluated by comparison to solutions of the GRACE-FO Science Data System (SDS) and COST-G ACs. We compare the mean spectral noise level of the solutions in terms of difference degree standard deviations (see Section 3.1). The spatial noise level is assessed by analyzing residual equivalent water height (EWH) signal over the oceans (see Section 3.1). To evaluate the signal content of the solutions, we compare the annual and semi-annual EWH amplitudes of major river basins, and annual mass loss trends in Greenland's drainage basins (see Section 3.2). Finally, in Section 3.3, GRACE-FO KBRR post-fit residuals are compared to those of the GRACE mission in frequency and spatial domains.
The presented results confirm that the outlined processing strategy is suited for obtaining monthly gravity field solutions with a quality competitive to that of the established ACs. The noise level assessment points out processing strategy related differences among the separate ACs, although the signal content is mostly not affected by these differences and is very similar for all ACs. The carried-out comparison of GRACE and GRACE-FO KBRR post-fit residuals highlights the overall improvement of the GRACE-FO K-Band ranging system performance. Nevertheless, several common systematic effects can be identified in GRACE and GRACE-FO post-fit residuals, e.g., higher noise related to shadow transitions, or spikes at multiples of the orbital and daily frequencies.

Gravity Field Recovery as a Generalized Dynamic Orbit Determination
The movement of an artificial satellite of the Earth is affected by the characteristics of its surrounding environment. By implication, the position (absolute, relative) and velocity coordinates-or the orbit of a satellite-contain information on a considerable number of parameters that are primarily but not exclusively describing physical and geometrical properties of the Earth. The time-variable gravity (TVG) is the main effect governing the motion of a LEO satellite. Modeling of the satellite-environment interaction allows us to derive the parameters describing Earth's TVG in terms of the gravitational potential V. The gravitational potential V at a location (λ, ϕ, r) on or above the Earth's surface can be represented as a synthesis of normalized spherical harmonic coefficients C nm and S nm , e.g., [49,50]: C nm cos mλ + S nm sin mλ P nm (sin ϕ) (1) where λ, ϕ are the spherical longitude and latitude in an Earth-centered and Earth-fixed frame, e.g., International Terrestrial Reference Frame (ITRF) [51]; r is the radial distance; GM ⊕ is the standard gravitational parameter of the Earth or the product of the universal gravitational constant G and Earth's mass M ⊕ ; R ⊕ is chosen as the semi-major axis of the Earth's ellipsoid; P nm are normalized associated Legendre functions; and n, m are degree and order of the spherical harmonic coefficients expansion. The objective of GFRand consequently of this work-is the estimation of the normalized spherical harmonic coefficients C nm and S nm . A set of estimated coefficients until a specific maximum degree n max , including the constants GM ⊕ and R ⊕ , is referred to as a gravity field solution. In the case of GRACE and GRACE-FO, the most common type of products are monthly gravity field solutions. Powerful tools for the analysis of the satellite-environment interaction are the dynamic and reduced dynamic orbit determination methods. These methods make use of orbit modeling, numerical integration and parameter estimation to solve the satellite's equation of motion, e.g., [52,53] where r andr are satellite's cartesian position and acceleration vectors in an inertial coordinate frame, e.g., Geocentric Celestial Reference Frame (GCRF) [51]; the first added − GM ⊕ r 3 r is the central body acceleration based on Newton's law of universal gravitation and the assumption that the satellite's mass is negligible if compared to Earth's mass M ⊕ ; and r P represents the sum of perturbing accelerations of gravitational and non-gravitational nature. For the precise dynamic orbit determination of LEO satellites and for GFR as applied in this work, the sum of perturbing accelerations can be formulated as follows: where the summation ∑ 8 i=1r i considers separate gravitational effects of tidal and non-tidal nature, as numbered (i) and described in Table 1. The accelerationr tvg is caused by Earth's TVG. Applying the Nabla operator ∇ = ( ∂ /∂x ∂ /∂y ∂ /∂z) T to potential V results in the corresponding acceleration in ITRF. The rotation matrix R gcrf itrf transforms the acceleration to GCRF. The accelerationr ng is caused by non-gravitational effects such as atmospheric drag, direct solar radiation pressure, Earth's albedo and thermal emission. A bias vector b and scale matrix S are needed for the calibration of the observed non-gravitational acceleration r ng,srf . The rotation matrix R gcrf srf formed from normalized quaternions transforms the calibrated non-gravitational acceleration from the satellite body-fixed science reference frame (SRF) to GCRF.  [43]; -Non-gravitational GRACE-D: alternative ACT products [57] The equation of motion is a vector form of an ordinary differential equation (ODE) of second order. Since an ODE of order n can be reduced to n first order ODEs, the equation of motion is related to a satellite state y = (r T ,ṙ T ) T containing the position vector r and the velocity vectorṙ via the following two first order ODEs, e.g., [52,53]: The integration of these ODEs results in a satellite state y. The six elements of the initial state vector y 0 = y(t 0 ) at time t 0 can be treated as integration constants. Due to the complexity of the perturbing accelerationsr P , the integration can not be performed analytically, i.e., a satellite state at an arbitrary time t can not be obtained directly. Several numerical integration methods can be used to achieve an approximative solution of satisfactory quality. A dynamic intermediate orbit consisting of several states y at epochs t can be obtained stepwise through the numerical integration of the two above-stated first-order ODEs, as follows: whereẏ = (ṙ T ,v T ) T combines the two first-order ODEs in a column vector. The intermediate dynamic orbit as a series of satellite states (y 0 , y(t 1 ), y(t 2 ), · · · ) is often referred to as a dynamically modeled or numerically propagated orbit. Even when assuming the best possible case for the integration constants in Equation (5), i.e., when the initial state is known very accurately, the dynamically modeled orbit will deviate from a true orbit in the course of time considerably. This is most of all due to uncertainties present in the models describing separate effects of the perturbing accelerationr P . In addition, also the numerical integration method, including specifics like the length and step size of the integration, contributes to a deviation of the dynamically modeled orbit from a true orbit. The concept of dynamic orbit determination is to adjust the dynamically modeled orbit to observations. Speaking generally, optimal values for the initial state and other dynamic parameters have to be found, so that the propagated orbit fits the observations in the best possible way. When additionally empirical parameters are introduced as unknowns, then the dynamic orbit determination approach becomes reduced-dynamic. Making use of these additional parameters allows a better fit of the numerically propagated orbit to observations. If the normalized spherical harmonic coefficients C nm and S nm are also part of the to be estimated dynamic parameters, the dynamic or reduced-dynamic orbit determination becomes a combined orbit determination and GFR.
The adjustment of the numerically propagated orbit to observations is usually accomplished by batch least squares adjustment. The linearized observation equations needed for the estimation of unknown parameters in GFR from GRACE and GRACE-FO sensor data can be summarized in the following simplified form: On the left side of this equation, the reduced observation vectors (observed−computed) can be found. These are formed as differences between the pseudo-observed orbit positions of the two satellites r C , r D or measured KBRRρ; and the dynamically modeled counterparts r C , r D ,ρ. Note that vec() is the vectorization operator, which converts a batch of orbit position vectors to a column vector. The right side of the observation equations consists of the partial derivatives of the dynamically modeled quantities with respect to (w.r.t.) the unknown parameters q. It is very common to divide the unknown parameters into a subset of local (∼) and global (⊕) parameters. The local parameter vector q ∼ consists of n ∼ elements and the global parameter vector q ⊕ is made up of n ⊕ quantities, among them the spherical harmonic coefficients of Earth's gravitational potential. Since dynamic orbit determination and GFR constitute a highly non-linear problem, the final parameters are obtained as the sum of a priori parameters and the estimated correction vectors ∆q ∼ and ∆q ⊕ in an iterative manner.
For the sake of completeness, clarity and a more simple reproducibility, the mathematical framework of the gravity field parameter estimation is treated in detail in the appendix. The appendix is divided into two parts: Appendix A.1 linear algebra and Appendix A.2 analysis. Appendix A.1 summarizes aspects of the parameter estimation, such as the parameter pre-elimination, combination of normal matrices and post-fit residuals computation. Appendix A.2 covers topics such as the linearization of observation equations and the formation of design matrices.

Processing Details
For performance reasons, the processing consists of two main steps. In a first step, an orbit pre-adjustment is performed without solving for the gravity field parameters. Estimated arc parameters from the pre-adjustment are used as a priori values in a second step, where the gravity field parameters are estimated along with the orbit in one iteration. The main characteristics of this two-step approach are summarized in Table 2 and outlined in detail below:  [58]. A straightforward implementation of this efficient integration approach is described in [14]. Forces of gravitational and non-gravitational nature affecting the motion of the satellites are modeled according to the information given in Table 1. An exception is made for the forces due to Earth's gravity field. In order to speed up the computation in this step, the gravity field is considered only until degree and order 120. The force modeling implementations were evaluated in a software comparison in the framework of COST-G. Implementations of all separate force effects agree well with the implementations of the COST-G ACs. The differences with regard to the implementations of other ACs are several orders of magnitude below 10 −10 m/s 2 [59]. • Arc length-The 3 h arc length employed in our approach differs considerably from the approaches of other ACs. The usual standard arc length employed by other ACs for numerical integration and GFR from GRACE and GRACE-FO data is 24 h. The aim of the rather short arc length is to allow a more precise orbit fit to the pseudo-observed positions and KBRR measurements, as inaccuracies, e.g., in force modeling, can be compensated by the frequent estimation of local arc parameters. In contrast to the very common arc length of one day, no constrained parameters, e.g., cycle per revolution accelerations, have to be co-estimated in order to achieve an adequate orbit fit. Since a decrease of the arc length increases the amount of arcs that can be processed independently, a considerable amount of processing time can be saved if parallel computing is utilized. • Observations-The reduced observation vectors are formed as differences between kinematic orbit positions or KBRR measurements and the corresponding quantities obtained from the dynamically integrated orbit (see Equations (A10) and (A12)). For KBRR measurements, the original 5 s sampling is kept. Kinematic orbits from the Astronomical Institute of the University of Bern (AIUB) are downsampled to 30 s and used as pseudo observations. Due to the general prevailing noisy character of kinematic orbits, compared to the rather smooth reduceddynamic orbits, a screening is performed. Epochs with a position difference larger than 8 cm w.r.t. the reduced-dynamic GNV1B orbits are not considered during parameter estimation. • Weights-The partials of the numerically integrated orbit w.r.t. unknown parameters are used to set up the design matrices. Then, technique-specific normal matrices are formed and combined. To set up the weight matrices, an initial standard deviation of 0.2 µm/s for KBRR is used. For the pseudo-observed position components, an initial standard deviation of 0.02 m is assumed. Variance component estimation, e.g., [60], is used to improve the technique-specific weights after each iteration of the orbit determination. • Parameters-Corrections to the initial state vectors and accelerometer bias parameters are estimated arc-wise based on least squares adjustment. Accelerometer scale factors, accelerometer shear and rotation parameters needed for a full scale matrix [23] are estimated monthly. For this purpose, the scale factors are fixed to 1 and the accelerometer shear and rotation parameters to 0. The objective of this step is not to obtain a best possible orbit fit, but rather to estimate appropriate initial values for the second step, therefore no empirical parameters are co-estimated. A priori values of the unknowns are corrected iteratively until convergence using the estimated corrections. A convergence is assumed when the mean of absolute KBRR reduced observation differences of two consecutive iterations is smaller than 0.1 µm/s . • Outlier arcs-After the pre-adjustment of all monthly arcs, an inspection is performed in order to detect spurious arcs that might disturb the GFR in the second step. Kinematic empirical KBRR parameters [61] consisting of 90 min biases and bias-rates, as well as 180 min periodic biases and bias-rates are fitted to each KBRR reduced observation vector. The fitted signal is subtracted from the KBRR reduced observation vectors. Then, the root mean square (RMS) of this difference is formed for each arc of a month. Finally, a sigma-based screening is applied to the time series of these quantities. Arcs outside the 3 sigma bounds are not considered in the further processing.
2. Orbit adjustment and gravity field recovery • Orbit modeling-Initial states and accelerometer biases from the pre-adjustment are used as new a priori values for the dynamic orbit modeling and computation of the state and parameter sensitivity matrices. Forces are modeled according to the description given in Table 1. Method specifics for the numerical integration are unchanged from step 1. • Observations-The formation of reduced observation vectors is consistent with the procedure in the pre-adjustment. • Weights-A standard deviation of 0.2 µm/s is used to set up the KBRR weight matrices. In case of kinematic positions, the inertial orbit covariance information is used to form diagonal weight matrices. We divide the elements of the kinematic positions weight matrices by an empirical factor of 25. Without this downweighting of the kinematic orbit covariance information, the quality of obtained solutions is unsatisfactory. A downweighting of GNSS-based observations w.r.t. KBRR measurements is an issue already known from the GRACE processing, e.g., [6,8], and deserves further attention in the future. • Parameters-The local parameters, i.e., initial states and accelerometer biases, are re-estimated in this step. In addition, the set of local parameters is extended by kinematic empirical KBRR parameters to absorb effects due to the possible mismodeling of perturbing accelerations. The set of kinematic empirical parameters is consistent with the definition utilized previously in the outlier arc detection. The normalized spherical harmonic coefficients of the monthly Earth's gravitational potential up to degree and order 96, as well as accelerometer scale factors, rotation and shear parameters, are introduced as global unknowns. The contribution of an arc to the set of global parameters is estimated after pre-elimination of the local parameters from the system of normal equations. The contributions of all 3 h arcs of a month are accumulated in order to obtain the final global parameters (see Equation (A3)). • Outlier screening-KBRR post-fit residuals are computed and a visual screening is performed in time and frequency domains (see Section 3.3). In case of additional outliers, the second step is repeated.

GRACE-FO Sensor Data and Products
For the computation of the gravity field solutions presented in this study, GRACE-FO Level-1B data products are used [43], which are generated by the SDS. The data products can be obtained from NASA's Physical Oceanography Distributed Active Archive Center (PO.DAAC) [62] and from GFZ's Information System and Data Center (ISDC) [63]. The following Level-1B data products are used in this work: • KBR1B: biased K-band ranges, as well as their first and second time derivatives K-band range-rates and K-band range-accelerations given in 5 s sampling. KBRR measurements are used as the main observations in the estimation process. The light time correction and antenna center offset correction as given in the KBR1B product are applied. • GNV1B: main data in these products are 1 s satellite positions and velocities in ITRF obtained from a reduced-dynamic orbit determination approach. In this work, the positions are used for modeling satellite accelerations caused by different forces (see Table 1). Instead of evaluating the accelerations at intermediate positions during every iteration of orbit determination, a major part of the accelerations is pre-computed using the precise GNV1B orbits. Only the acceleration caused by the Earth's gravita-tional potential is evaluated at every intermediate position during orbit determination and GFR. • SCA1B: 1 s normalized quaternions describing the rotation between SRF and GCRF.
Since the numerical orbit propagation is accomplished in an inertial frame, the quaternions are needed for transforming calibrated non-gravitational accelerations to GCRF (see Equation (3)). • ACT1B: main data in these products are 1 s linear accelerometer measurements given in SRF. The measurements represent the sum of acceleration variations caused by non-gravitational effects. Accelerometer calibration parameters have to be estimated during orbit determination and GFR. The ACT1B products replace ACC1B products that were formerly used for GFR from GRACE data. The main feature of ACT1B is a so-called transplantation of GRACE-C accelerometer measurements to satellite GRACE-D. The necessity for this transplantation [24] arises because of a severe degradation of GRACE-D measurements, e.g., [40]. With the exception of June 2018, only GRACE-C ACT1B products are used in this work.
In addition to the listed products, the following GRACE-FO data products are utilized: • Alternative ACT products: For GRACE-D, for all months except June 2018, alternative ACT products [57] from the Institute of Geodesy at Graz University of Technology are used instead of the official ACT1B products. • AIUB kinematic orbits: These kinematic orbits are produced at the Astronomical Institute at University of Bern. Processing details are summarized in [64]. The kinematic orbits do not contain any information from dynamic models. The positions are treated as pseudo-observations during parameter estimation. In addition to the positions, cofactor matrices are also available. These matrices are used to form the corresponding weight matrices.

Results and Evaluation
To evaluate the quality of the computed monthly gravity field solutions, the LUH-GRACE-FO-2020 time series is compared to solutions of the SDS and COST-G ACs. In Section 3.1, the noise level of the solutions is assessed and compared in spectral as well as spatial domains. Annual and semi-annual EWH amplitudes of river basins and annual mass loss trends in Greenland's drainage basins are compared in Section 3.2. Typical GRACE-FO and GRACE KBRR post-fit residuals are shown and compared in Section 3.3. Note that because of the so far short GRACE-FO operation time, the below presented comparisons constitute only a very preliminary image. The LUH solutions are compared to the following publicly available solutions: All these solutions are available on the website of the International Centre for Global Earth Models (ICGEM) [72]. The regarded solutions allow a fair comparison, since the processing is based on variants of the generalized dynamic orbit determination with variational equations. In addition, all compared solutions are unconstrained, i.e., computed without applying regularizations.

Noise Level
First, the overall noise levels of the solutions are compared in the spectral domain in terms of mean difference degree standard deviations (DDSDs) with regard to a reference model. Since in general, the normalized spherical harmonic coefficients of a gravity field solution may be scaled with different standard gravitational parameters of the Earth GM ⊕ and Earth reference radii R ⊕ , solutions of all ACs are re-scaled to the standard values from the IERS Conventions 2010 [51]. One common reference model is computed as the mean of all monthly GRACE-FO solutions available until now (June 2018 until December 2020). Before computing the reference model, the C 20 and C 30 coefficients of all solutions are replaced, with more accurate values obtained from satellite laser ranging (SLR) [73,74]. The reference model is then subtracted from all solutions before computing monthly DDSDs. Note that the DDSDs of the low degree coefficients are dominated by signal (approximately until degrees 20-30), while higher degrees are dominated by noise.
The averaged DDSDs are illustrated in Figure 1. A high level of consistency can be observed for the low degree coefficients that are of major importance for mass change estimations due to a large signal content. Larger differences are observed for degree 2, among which the LUH and AIUB solutions have the smallest DDSDs. The smaller DDSDs are primarily caused by less noisy C 20 coefficients, meaning that, in general, the C 20 coefficients are closer to the more precise SLR coefficients. Since the GRACE-FO C 20 coefficients are usually replaced with SLR coefficients, this smaller noise is not of great relevance for the estimation of mass variations. Starting at around degree 12, the different processing strategies of the ACs are causing noticeable deviations of the DDSDs. The LUH solutions slightly outperform the GFZ solutions and have a noise level similar to that of the JPL solutions along major parts of the spectrum. AIUB's solutions based on the celestial mechanics approach, e.g., [75], have slightly less noise in the coefficients between degrees 30 and 70, compared to GFZ, JPL and LUH. The ITSG solutions that incorporate an advanced stochastic modeling and additional tidal estimates, e.g., [7,76], have the overall smallest noise level. Moreover, the solutions from CSR suppress noise better than most ACs. This might be due to the fact that in CSR's GFR strategy, the local parameters are estimated beforehand. During the estimation of spherical harmonic coefficients of Earth's gravitational field, the local parameters are not re-estimated [4]. The noise level of the solutions in the spatial domain can be assessed by analyzing residual mass variations in suited regions. A suited region is defined as an area where only relatively small mass variations are present, e.g., oceans. For comparison of the noise characteristics in the spatial domain, the normalized spherical harmonic coefficients are converted to 2 • × 2 • global gridded mass variations in terms of EWHs, e.g., [77,78]. The EWH values are computed with regard to the mean model that was already used for the mean DDSDs. A Gaussian filter [77] with an averaging radius of 400 km is applied in order to mitigate the typical meridional striping. A model consisting of a bias, trend, annual and semi-annual variation is fitted to each grid cell of the EWH time series to absorb time-variable signal. After reduction of this simplified model, the residual signal over the oceans is rather homogenous, while several regions on land show the presence of unabsorbed hydrological signal, e.g., in Africa and South America (see Figure 2). Since solely the areas over oceans are of interest, only values inside an ocean mask are considered. The comparison of the monthly RMS noise over the oceans can be seen in Figure 2. The monthly RMS noise over the ocean values highly correlate with the earlier presented DDSDs. Particularly striking are two ocean noise RMS values in the ITSG time series, i.e., October 2018 and February 2019, where the ITSG values are slightly larger, compared to the very small noise level during other months. This can be explained by larger sensor data gaps during these months. In contrast to ITSG, the other ACs utilize additional sensor data from neighboring months in order to stabilize their solutions during these months.

Signal Content
Annual and semi-annual amplitudes in terms of mean EWHs of about 180 major river basins are compared. The geographical boundaries of the river basins were taken from the Total Runoff Integrating Pathways (TRIP) product [79,80]. The global EWHs are computed according to the same procedure earlier used for the assessment of the spatial noise. Note that the C 20 as well as C 30 coefficients were replaced with SLR-derived values. To each 2 • × 2 • cell, a model consisting of a bias, trend, annual and semi-annual variation was fitted. Then, the mean annual and semi-annual EWH amplitudes in the boundaries of each river basin region were formed (see Figure 3a   Next, mass loss trends in six drainage basins of Greenland (see Figure 3d) are compared. The definitions of the drainage basin region boundaries were taken from the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE) dataset [81]. The mean EWH trends are converted to mean mass trends expressed in Giga tons per year (Gt/yr). Degree 1 coefficients were not added. Note that the mass trends are neither corrected for glacial isostatic adjustment nor the leakage effect. The obtained annual mass loss trends are illustrated in Figure 3c. The trend estimates of all ACs show a high level of agreement for all six regions of Greenland. The highest discrepancy in the trend amplitudes of about 3.35 Gt/yr can be found in the SW basin. On average, the maximum discrepancy is about 2 Gt/yr.

KBRR Post-Fit Residuals
In this section, three months of exemplary GRACE KBRR post-fit residuals (January 2008 until March 2008) are compared to three months of exemplary GRACE-FO KBRR post-fit residuals (September 2019 until November 2019). The utilized GRACE gravity field solutions were computed according to the processing strategy summarized in [13,14]. The GRACE-FO KBRR post-fit residuals are routinely examined and are an important component of the outlier screening strategy of the LUH-GRACE-FO-2020 solutions (see Section 2.2). The KBRR post-fit residuals are computed according to Equation (A7). It should be emphasized that for GRACE-FO, a revised processing is applied. The most important differences concern the utilized background models for the ocean tides (GRACE: EOT11a [18], GRACE-FO: FES2014b) and non-tidal short-term variations of the oceans and the atmosphere (GRACE: AOD1B RL05 [20], GRACE-FO: AOD1B RL06). Moreover, in the case of GRACE-FO, a full accelerometer scale matrix is estimated.
In the following, we inspect spatial and temporal characteristics of KBRR post-fit residuals in time-argument-of-latitude (TAL) diagrams, the argument of latitude being the geocentric angle in the orbital plane between the ascending node and the position of the spacecraft. In TAL diagrams, systematic effects that depend on the orbital configuration (e.g., the position and orientation of the satellites with regard to Earth and Sun) tend to stand out as distinct and coherent structures. We focus on two bands of band-pass filtered KBRR post-fit residuals. The first band (A) extends from 0.6 mHz to 5 mHz. This is above the frequencies where empirical parameters directly absorb errors and unmodeled signal contributions. In band A, we expect the geophysical aliasing of unmodeled tidal and non-tidal mass variations, as well as slowly drifting instrument systematics as main contributors to residuals. The second band (B) extends from 5 mHz to 20 mHz. In this band, the influence of mass variations to residuals is fading out due to the limited sensitivity of the GRACE-FO constellation, however, it is a band that is well suited to detect systematics due to changes in the orbital configuration and/or the instrument operation.
The TAL diagrams of these two bands can be seen in Figure 4. The post-fit residuals in the 0.6 to 5 mHz frequency band differ in their intensity. The GRACE-FO residuals are distinctively smaller. Several geographical patterns of different amplitudes and periodicities are visible in this band. Some of these features are, for example, daily patterns related to the Earth rotation and patterns with a period of approximately 30 days, which are more distinct around degrees 90 and 270, i.e., at the poles. The features in the 0.6 to 5 mHz frequency band might be caused by uncertainties in models of rapid tidal and non-tidal mass variations. The modeling of the rapid tidal and non-tidal mass variations is still one of the major error contributors in GFR. The smaller GRACE-FO residuals in this band can be explained to a certain extent by earlier mentioned updates in background force modeling.
The appearance of the residuals in the 5 to 20 mHz frequency band differs considerably for the two missions. Several systematics are present in GRACE post-fit residuals, e.g., during the entering and exit phase into and from the penumbra region (1.) [82], frequency-related KBR system signal-to-noise ratio drops (2.) [83,84], baffle-related KBR system signal-to-noise ratio drops (3.) [83,84], systematics possibly related to the star camera assembly (4.). In addition, a vertical line can be seen (5.) that is caused by bad sensor data in the involved arcs. In the computation of the LUH-GRACE-FO-2020 solutions, such artifacts are detected and removed during outlier detection. Except the patterns related to the entering and exit phase into and from the penumbra region (1.), the GRACE-FO residuals in the 5 to 20 mHz frequency band do not explicitly reveal any of the previously listed systematic errors. Nevertheless, very slight systematics which manifest as horizontal lines are recognizable at approximately 90 and 270 degrees (2.). These horizontal lines correlate with systematics in the 0.6 to 5 mHz frequency band.
Noticeable is the overall amplitude difference of the post-fit residuals. The corresponding unfiltered RMS during the three examined GRACE-FO months (8.0 × 10 −8 m/s) is approximately three times smaller than the RMS during the three GRACE months (3.0 × 10 −7 m/s). These findings are consistent with the results presented in [40].  Figure 5 illustrates the power spectral density (PSD) of the KBRR post-fit residuals in the frequency band from 10 −4 to 10 −1 Hz. GRACE-FO amplitudes are distinctively smaller along the whole spectrum, compared to GRACE. Very prominent is the difference in the frequency range between 10 −2 and 10 −1 Hz. The smaller GRACE-FO amplitudes in this frequency band are related to the lower noise level of the GRACE-FO KBR system, e.g., [40]. Highlighted with red vertical lines are the first 15 orbital frequency multiples (1 cycle per revolution (cpr) to 15 cpr). The PSD of the residuals is dominated by spikes located at orbital frequency multiples. As expected from Figure 4a,b, spikes are also present at multiples of the daily frequency, although they only become visible when an appropriate zooming is applied to the PSD figure. In contrast, the smaller residuals around 1 cpr can be regarded as exceptions, since they are absorbed by the co-estimation of arc-wise empirical KBRR parameters (see Section 2.2). The estimation of additional empirical KBRR parameters in order to absorb further frequencies is currently under investigation.

Conclusions
The Institute of Geodesy at Leibniz University Hannover produces unconstrained GRACE-FO monthly gravity field solutions. The solutions are computed with the MATLABbased GRACE-SIGMA software package recently developed at Leibniz University Hannover. The solutions are obtained from GRACE-FO Level-1B products, alternative accelerometer products produced at Graz University of Technology, and kinematic orbits from the Astronomical Institute at University of Bern. The regularly updated time series named LUH-GRACE-FO-2020 is accessible on the website of the International Centre for Global Earth Models [85] and in LUH's research data repository [86].
The quality of the solutions is competitive with those of the GRACE-FO Science Data System and Combination Service for Time-variable Gravity Fields analysis centers. While the spectral and spatial noise levels of the separate analysis centers slightly differ, the signal content of the solutions is very similar among all analysis centers. The C 20 and C 30 coefficients were excluded from the comparison of the signal content and therefore deserve further attention in the future. The GRACE-FO K-band range-rate (KBRR) post-fit residuals are about three times smaller, compared to GRACE. Most pronounced systematics in GRACE-FO KBRR post-fit residuals are related to the entering and exit phase into and from the penumbra region. The power spectral density of the post-fit residuals is mainly dominated by spikes located at multiples of orbital frequency. The analysis and further understanding of the systematics in the post-fit residuals are important for identifying factors that limit the quality of gravity field recovery.  [85]; and in the LUH Research Data Repository at https://data.uni-hannover.de/ dataset/luh-grace-fo-2020, accessed on 24 March 2021 [86].

Acknowledgments:
We would like to thank the German Space Operations Center (GSOC) of the German Aerospace Center (DLR) for providing continuously and nearly 100% of the raw telemetry data of the twin GRACE-FO satellites. Majid Naeimi is acknowledged for developing the initial version of GRACE-SIGMA. The colleagues from the Combination Service for Time-variable Gravity Fields (COST-G) are acknowledged. The COST-G meetings helped us to improve our gravity field recovery strategy. The International Space Science Institute (ISSI) in Bern/Switzerland is acknowledged for hosting the COST-G team meetings in 2019, 2020 and 2021. Participation in these meetings was financially supported by ISSI. We are thankful for the valuable comments of the anonymous reviewers. The publication of this article was funded by the Open Access Fund of the Leibniz Universität Hannover.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used repeatedly in this manuscript: One of the most important parts of orbit determination and gravity field recovery is parameter estimation. According to weighted least squares adjustment, the parameter vectorx containing a set of unknowns can be obtained as follows, e.g., [87,88]: where A is the design matrix; P is the weight matrix; l is the observation vector; N is the normal matrix; b the right hand side vector. Because of the non-linearity of the orbit determination, approximate values of the unknowns x 0 are introduced and corresponding corrections ∆x are estimated iteratively, so that the parameter vector is defined asx = x 0 + ∆x. The presence of different types of unknowns, i.e., local parameters such as the initial states; and global parameters like the spherical harmonic coefficients of a gravity field solution, allows one to divide the parameter correction vector into two parts: ∆x = (∆x T ∼ , ∆x T ⊕ ) T , where local and global parameters are denoted with the subscripts ∼ and ⊕, respectively. In case of a set of arcs i = 1, 2, · · · , j, the parameter correction vector can be extended to ∆x = (∆x T ∼1 , ∆x T ∼2 , · · · , ∆x T ∼j , ∆x T ⊕ ) T . Applying this separation of parameters to Equation (A1) leads to the below stated system of normal equations that has to be solved, e.g., [60]: The inversion of this system of normal equations according to Equation (A1) can become heavily memory intensive and unsolvable. However, the final estimate of the global parameter corrections ∆x ⊕ can be obtained after pre-elimination of all local parameter corrections, e.g., [60]: where KBRR observations have to be combined with GNSS tracking data, e.g., kinematic orbits, in order to obtain a gravity field solution of satisfactory quality. The normal matrices N ∼i , N ⊕i , N ∼⊕i and the right hand side vectors b ∼i , b ⊕i can be formulated as techniquespecific weighted superimpositions of observations, e.g., [89]: where the contribution of satellite-specific positions is denoted with the subscripts C, D and the contribution of KBRR with the subscript K; the vectors ∆l Ci , ∆l Di and ∆l Ki are reduced observation vectors; these vectors as well as the design matrices A ∼Ci , A ∼Di , A ∼Ki , A ⊕Ci , A ⊕Di and A ⊕Ki are defined in Appendix A.2. P Ci , P Di and P K are the weight matrices representing the stochastic model of the observations. From the subscripts of the weight matrices, it can be seen that the weight matrices for kinematic positions are different for every arc i, whereas no arc-specific weights are applied for KBRR (see Section 2.2). After the final global parameter corrections ∆x ⊕ are estimated, the correction values for the local parameters of an arc i can be obtained, e.g., [60]: Finally, the KBRR post-fit residuals v Ki of an arc i are obtained as follows:

. GFR Parameter Estimation: Analysis
The relationship between a set of measurements l and a set of unknown parameters x can be established via a functional model f . Errors are accumulated in the residual vector v, resulting in the following observation equation for a linear functional model: The functional model of the dynamic orbit determination is not linear with regard to the unknowns, so a linearization of f (x) is mandatory. The linearization is carried out by means of a first-order Taylor series expansion that is then evaluated with the a priori values of the unknown parameters x 0 , e.g., [87,88]: The partial derivatives of Equation (A8) can be summarized in a design matrix A: In the case of position components, the reduced observation vectors of satellites C and D are defined as ∆l Ci = vec(r Ci − r Ci ) and ∆l Di = vec(r Di − r Di ).r C ,r D contain inertial kinematic positions and r C , r D the dynamically modeled positions, i.e., the outcome of the numerically integrated equation of motion (see Equation (4)). The subscripts i denote that several positions of an arc i are considered. The reduced KBRR observation vector of an arc i is accordingly defined as ∆l Ki =ρ i −ρ i . Vectorρ i contains the observed KBRR measurements that are part of the KBR1B products. The counterpartρ i is calculated from the dynamically modeled positions r Ci , r Di and velocitiesṙ Ci ,ṙ Di as a projection of the relative velocity vectorsṙ CDi to the line-of-sight connection e CDi . One computed KBRR is defined as, e.g., [52,90]:ρ =ṙ CD · e CD (A10) Let us define two vectors containing a set of unknowns related to the concurrent arcs i of satellites C and D as: q ∼i = (y T 0Ci , b T Ci , y T 0Di , b T Di , e T i ) T and q ⊕i = (C T nmi , S T nmi , s T Ci , s T Di ) T . Vector q ∼i contains local parameters that are related to both satellites. Herein, y 0 represents initial states, b accelerometer biases and e empirical KBRR parameters. Since local parameter corrections ∆x ∼i can only be obtained when the corrections to the global parameters are estimated, the vector q ∼i represents only intermediate local parameters. The vector q ⊕i contains the spherical harmonic coefficients of Earth's gravitational potential: C nmi = (C 20 , C 21 , C 22 , · · · , C n max n max ) T , S nmi = (S 21 , S 22 , S 31 , · · · , S n max n max ) T . s are the elements of the full scale matrix. Relating observations and unknown parameters according to the general linearization Equation (A8) and omitting the residuals results in: vec(r Ci − r Ci ) = where ∆q ∼i and ∆q ⊕i are the corrections to the a priori values of parameters in vectors q ∼i and q ⊕i ; and n ∼ and n ⊕ are the number of parameters in vectors q ∼i and q ⊕i , respectively.
The partial derivatives of the positions r Ci , r Di w.r.t. the corresponding initial states y 0Ci and y 0Di are part of the so called state transition matrices Φ. A state transition matrix contains the partial derivatives of a satellite state y at a specific time t w.r.t. the initial state y 0 at time t 0 , e.g., [53,91]: The partial derivatives of the local and global dynamic parameters, namely: p Ci = (b T Ci , s T Ci , C T nmi , S T nmi ) T for satellite C, and p Di = (b T Di , s T Di , C T nmi , S T nmi ) T for satellite D, are part of the so called sensitivity matrices S. A sensitivity matrix contains the partial derivatives of a satellite state y at time t w.r.t. the dynamic parameters p, e.g., [53,91]: In case of precise dynamic orbit determination, these partials can not be obtained analytically. An elegant way is a combined numerical integration of the two ODEsΦ anḋ S with the satellite's equation of motion using the integration constants Φ 0 = I and S 0 = 0, e.g., [53,91] The partial derivatives of one computed KBRRρ w.r.t. the initial states y 0Ci , y 0Di and dynamic parameters p Ci , p Di can be paraphrased as: