GRACETOOLS—GRACE Gravity Field Recovery Tools

: This paper introduces GRACETOOLS, the ﬁrst open source gravity ﬁeld recovery tool using GRACE type satellite observations. Our aim is to initiate an open source GRACE data analysis platform, where the existing algorithms and codes for working with GRACE data are shared and improved. We describe the ﬁrst release of GRACETOOLS that includes solving variational equations for gravity ﬁeld recovery using GRACE range rate observations. All mathematical models are presented in a matrix format, with emphasis on state transition matrix, followed by details of the batch least squares algorithm. At the end, we demonstrate how GRACETOOLS works with simulated GRACE type observations. The ﬁrst release of GRACETOOLS consist of all MATLAB M-ﬁles and is publicly available at Supplementary Materials.


Introduction
The Gravity Recovery and Climate Experiment (GRACE) mission [1] has started a new era in our understanding of the Earth's gravity field temporal variations.There are several methods to recover Earth's gravity field from GRACE inter-satellite ranging observations.Classically, a distinction has been made for gravity field recovery approaches using satellite observations: timewise and spacewise approaches [2].A comprehensive review of timewise and spacewise approaches can be found in [3].Table 1 summarizes these approaches.
Table 1.Methods of Earth's gravity field recovery from GRACE inter-satellite ranging observations.

Observations Spacewise Observations
Classical approach ρ, ρ Energy balance approach ρ e.g., [1,4] e.g., [5,6] Celestial mechanics approach ρ, ρ Acceleration approach ρ e.g., [7,8] e.g., [9] Short arc approach ρ, ρ Line of Sight Gradiometry ρ ρ e.g., [10] e.g., [11] Spacewise approaches model the observables as a function of spatial coordinates, leading to spherical harmonic analysis of data on a sphere.In other words, they project in-situ observations, in the case of GRACE, range (ρ), range rate ( ρ) and range acceleration ( ρ), that are observed at varying orbital heights to the height of a mean orbital sphere.Spacewise methods are mostly conducted in two steps.First, an observable is calculated along the orbit.Then, the geopotential spherical harmonic coefficients are computed using least squares estimation.
On the other hand, timewise approaches treat the data as a time series along the orbit.They solve variational equations, in the case of GRACE, the equation of motion for two orbiting satellites, and observation equations for inter-satellite ranging observations.Consequently, the solution consists of numerical integration of the equation of motion and improving the satellite trajectory and gravity field parameters by least squares solution of the observation equation.This makes gravity field recovery by variational equations conceptually complicated and computationally demanding.
Although the estimated gravity fields based on various methods are similar, the differences in processing strategies and tuning parameters result in solutions with regionally specific variations and error patterns.While previous papers have compared GRACE final gravity solutions by different data processing centres (e.g., [12]), we aim to find out how and why these differences occur.For this reason, we initiated an open source GRACE data analysis platform, where the existing algorithms and codes for GRACE data analysis are transparently shared and improved.
Section 2 summarizes the fundamental mathematical models for solving variational equations.In Section 3, we rewrite variational equations in a matrix format for two GRACE satellites.Section 4 gives step by step GRACETOOLS batch processor algorithm for solving variational equations with two GRACE satellites.In Section 5, we present the results of numerical computations using GRACETOOLS with simulated data.

Mathematical Models
Two main mathematical models are required for the orbit determination and gravity field recovery resulting in variational equations.First, equations describing the satellite dynamics are necessary to represent the current, best knowledge of the space environment, Earth's processes, and the satellite's behavior.In addition, an observation model is required to relate the evolution of the satellite's state to the measured observable(s).This section summarizes a solution to the variational equations that relates the equation of motion to the observation model.The content of this section is mainly borrowed from [13][14][15].

Equations of Motion and Observation Models
We begin by representing the dynamic and observation models as functions of the state vector, X(t).As usual in parameter estimation, it not only contains the satellites states position and velocity, but additionally all model parameters to be estimated such as geopotential coefficients.The equations of motion for two GRACE satellites can be generally written as a system of first order ordinary differential equations (ODE) ˙ The observations Y can be expressed with the state vector X(t) and a suitable function where represents the observation errors.The dynamics of orbit mechanics, described by Equations ( 1) and ( 2), are generally highly nonlinear.In orbit determination or gravity field recovery, the equations of dynamic and observation models are linearized at a reference, or nominal trajectory that is close enough to the true trajectory to allow linearization.We are not directly estimating the values of the model coefficients, X(t), but rather the deviation between the true and nominal values.These deviations of the state and observations are defined as where the * denotes the nominal values.Since we have assumed that the nominal trajectory is within close proximity to the true solution, we can expand y and x about X * via Taylor expansion. with The same can be done for the observation equation with

The State Transition Matrix
The general solution to the system of differential equations in Equation ( 7) can be expressed as where t 0 is some specified epoch and Φ is called the state transition matrix.It satisfies the following conditions Φ(t, t 0 ) = A(t)Φ(t, t 0 ), Φ(t 0 , t 0 ) = I.
The state transition matrix maps deviations in the state vector from a time t 0 to t.With a given matrix A, the differential equation can be solved, at least numerically, and Φ(t, t 0 ) can be computed for every epoch.
The state transition matrix is valid as long as it stays within the linear regime or (t − t 0 ) is small enough, i.e., five seconds for GRACE.In the case of GRACE, real data processing staying within linearity is highly dependent on how accurately the background force models are modeled.The accuracy of the state transition matrix depends on the accuracy of the numerical integrator.High accurate multistep numerical integrators are provided in GRACETOOLS (see Supplementary Materials), with a brief description in Appendix A.

The Normal Equations
The primary functions of least squares estimation is to fit a model to a set of observations.For example, given the following system we would like to find the parameters, x, which come closest to representing the observed measurements in y.One way to accomplish this is to treat the system as an optimization problem and define a performance index that can then be minimized.For the least squares method, the performance index is chosen to be the sum of squares of the residuals, or errors.Defining the error to be the performance index becomes For the case of weighted least squares, we can introduce the weight matrix into Equation ( 14) to obtain Minimizing the performance index is done by taking the 1 st derivation of the performance index (i.e., setting The result, Equation ( 17), is commonly referred to as the normal equations.If H consists of at least m linearly independent observations, then the normal matrix, H T W H, is both symmetric and positive definite.The condition also implies that the inverse H T W H exists, allowing us to solve for x.A good approximation of the weight matrix W can be derived from post fit residuals [16].The algorithm for implementing the normal equations approach using the unit matrix for W involves the accumulation of m rank-one updates, one for each observation: Once the contributions from all observations have been accumulated, the resulting normal matrix can be inverted via Cholesky decomposition, or other adequate methods, to get the solution for x.

Partitioned Normal Equations
The typical scenario for the GRACE parameter estimation is to have data files for each arc, where an arc is defined as a specific length of time, typically one day [13].The parameters to be estimated are categorized into two different groups, or levels: • Local: Parameters that are valid for only one arc, such as initial position and velocity for each day • Global : Parameters that are valid across all arcs, such as monthly spherical harmonics coefficients Consider the following partitioning of the generalized state vector z, and observation-state mapping matrix where H x is the local contribution and H c is the global contribution.To solve such a system, we need to divide, or partition, the local and global parameters so that each group may be solved for separately.
Taking the matrix from Equation (20) and inserting into Equation ( 17), we obtain From this, we define the following Inserting these expressions into Equation (21), we have Multiplying the above equations, we get Solving for x in Equation ( 23) and inserting this result into Equation (24) gives us an expression for the global estimates: The above scenario applied only to a single arc, but it is not difficult to extend the idea to incorporate any number of arcs, each with their own set of local parameters.H z would look like k is the number of arcs, or days.Again, the local contributions are independent of other parameters, and this explains their location along the diagonal.Inserting Equation (27) into Equation (17), we obtain an expression for the generalized partitioned normal equations: These can be arranged in a similar fashion as the single arc case to solve for the global parameters: While the partitioned normal equation method is slightly challenging to implement, it takes advantage of the structure of the matrix in Equation ( 27) and avoids unnecessary operations with zeros.

Matrix Form
In this section, we present the three main matrices A, Φ and H in their matrix format, and consider their components in more detail.

Matrix A: State Partials
For two GRACE satellites, the state vector X and its derivative Ẋ are X is set of unknown parameters including satellite's position r and velocity ṙ, and the Earth's gravity field spherical harmonic coefficients K lm .K lm includes both sine and cosine related coefficients of the spherical harmonic series.Klm equals to zero during the recovery period of one month, which means, during one month, the coefficients are assumed to be constant.Then, which is a (6 + 6 + n) × (6 + 6 + n) matrix, assuming there are n spherical harmonic coefficients to be estimated for one month.

Matrix Φ: State Transition Matrix
For two GRACE satellites, the state transition matrix is with the dimension of (6 + 6 + n) × (6 + 6 + n).

Matrix H: Observation Partials
The H matrix is given by which is a (6 + 6 + n) vector.The relative position and velocity of GRACE A and GRACE B are calculated as follows: and range and range rate are calculated according to: where • is the vector dot product, ṙAB is relative velocity vector and e AB is the line of sight unit vector of two GRACE satellites.Accordingly, the partials of the leading satellite (A) can be written for range rate measurements as: The partials of the trailing satellite (B) are the same as ones of the leading satellite, but with opposite sign [9].The partials of the range rate with respect to K lm are [14]: where the partial derivatives are components of the state transition matrix.

GRACETOOLS Batch Processor Algorithm
There are several ways to perform orbit determination and gravity field recovery; batch processor and Kalman filter are the most famous algorithms [15].For the first release of GRACETOOLS, we implemented the batch processor approach.In the following, the step by step algorithm for the batch processor is given: (1) Initialize at t 0

•
Read the initial state vector, i.e. position and velocity of two GRACE satellites (( X * 0 ) k ).They are the first rows from two GPS navigation level 1B (GNV1B) daily files; GRACE level 1B orbits are given in an Earth-fixed frame and they need to be transformed to the inertial frame.

•
Read the range rate observation Y from the K-band ranging level 1B (KBR1B) daily files.

•
Read a-priori gravity model in terms of spherical harmonic coefficients ( K * lm ).
• The a-priori deviation values for state and spherical harmonics coefficients and their associated error covariance matrix are: Equations (6.3.22),(6.3.23), and (6.3.25) from [15] give (2) Supply the numerical integrator with the following vector at each time point The first four elements provide the reference orbit X * (t), and the last one yield the elements of Φ(t, t 0 ).The reference orbit is used to evaluate A(t), which is needed to evaluate Φ(t, t 0 ).
Accumulate current observation • Build H i from Equation (33) using the last column of Equation (32), and then • Partition H i into (H x ) i for initial state and (H c ) i for spherical harmonics coefficients.
Repeat for each day and save M xx , M xc , M cc , N x , N c for each day.Save y,H x ,H c for each day, to plot daily post fit range rate residuals.(5) Solve normal equations First, for global parameters, spherical harmonics coefficients and then for local daily parameters, initial state of the two satellites for each day Estimate postfit range rate residuals for each day Update the initial state of both satellites for each arc (or day here) and gravity field coefficients: Repeat this until the least squares estimation is converged or it is below an accepted error tolerance.

Evaluation with Simulated Data
In this section, the validity and quality of GRACETOOLS is tested by simulated data.This offers the possibility to arbitrarily set all input data and their uncertainties.With error free measurement data, it is possible to recover the true gravity field from an initial guess, which proves the functionality of GRACETOOLS.GRACETOOLS can be adapted for gravity field recovery using real GRACE and GRACE Follow-On L1B data.
For this example, we integrate the orbital trajectories with initial states of GRACE A and GRACE B for four days, and with Global Gravity Model (GGM05S) [17], degree and order 10, as the true field.We use GRACE data sampling of five seconds.From these simulation, observations (range rate, position and velocities) are computed in the form of KBR1B and GNV1B data files.As an initial guess for the gravity field parameters in the estimation, we add noise to the true gravity field.We use normally standard distributed random numbers n lm and calculate a deviation to the true gravity field spherical harmonic coefficients by Subsequently, we use the corrupted coefficients K lm,err + K lm as initial gravity field in our estimation with the perfect simulated observations.Because no noise and other unknown disturbances are corrupting the observation data, we are able to recover the initial GGM05S gravity field, with a sufficient number of iterations.According to Equation (52) with L min = 2 and L max = 10, we need to estimate 117 spherical harmonics coefficients.
In Figure 1, the original (true) gravity field (GGM05S), which is used to produce the observations and the disturbed field, which is used as initial guess in the estimation process and the difference between both gravity fields are shown.All gravity fields are shown in the form of square root of degree variances.Additionally, for comparison, the difference between GGM05S and Earth Gravitational Model 1996 (EGM96) [18] is also plotted.This shows that the imposed initial gravity field error is quite big.
The true (GGM05S) and the initial gravity fields are shown spatially in terms of gravity anomalies in Figure 2.Both plots show the gravity fields up to degree and order 10.In these plots, the initial error is visible quite distinctly.The big error on the C 20 coefficient overlays nearly all other features of the gravity field.As comparison the difference between GGM05S and EGM96 is also shown.To recover the gravity field, we use four days of simulated data and a local arc length of one day.The perfect initial states for each arc are used.The difference between the estimated gravity field after each iteration with respect to the true field is shown in Figure 3 for 28 iterations.While the lower degrees (2 and 3) improve with the first iteration, the higher degrees are getting a bit worse and start improving after some iterations because there is a big difference between the initial and true gravity field for the lower degrees (2 and 3).After the fifth iteration, all degrees are improving subsequently with further iterations up to a square root of degree variances of about 2 × 10 −14 .After this value is reached, further iterations are not improving the result anymore.Because the deviations of the spherical harmonic coefficients ( ˆ c) have reached the order of 10 −15 to 10 −16 , which is the border of the double data type precision.Until about Iteration 12, the resulting range rate residuals are decreasing.After Iteration 12, they stay at the same level of about some 10 −10 m/s.Exemplarily, the range rate residuals for Day 1, for the sixth and last iteration (Iteration 28) are shown in Figures 4 and 5.The evolution of the range rate residuals over all iterations is displayed in Figure 6, also for Day 1.The root mean square (rms) of the range rate residuals is computed and plotted over all iterations.This demonstrates that the residuals do not decrease further than 10 −10 m/s.This is reasonable because this is about the numerical integration precision that is possible to archive for velocities (or relative velocities) with double precision data type.However, the gravity field coefficients can still be improved because, at this point, they are not limited by the integration error, but the linearization.The linearization error can be further reduced with subsequent iterations.

Conclusions
In this work, we have presented the first version of our new GRACE gravity field recovery code, written in MATLAB, using variational equations, which we call GRACETOOLS.The current version of our code provides an environment to evolve a gravity field recovery tool using GRACE satellites type observations.This first release of GRACETOOLS has been mainly focused on recovering a gravity field from simulated GRACE range rate observations.The on-going development of GRACETOOLS includes more sophisticated end to end simulations with all input force models and GRACE Follow-On instrument observations.We hope that, with the transparency provided by opening up the code, we can address open questions in GRACE data analysis faster and more efficiently.

Figure 1 .
Figure 1.Square root of degree variances from true gravity field (GGM05S), with which observations were simulated, initial gravity field for the estimation process, and the difference between both fields.As comparison the difference between GGM05S and EGM96 is also shown.

Figure 2 .
Figure 2. True (GGM05S) gravity field (left) and initial gravity field (right) in terms of gravity anomalies up to degree and order 10.

Figure 3 .
Figure 3. Square root of degree variances of difference between true gravity field and estimated fields after each iteration.

Figure 4 .
Figure 4. Time series of range rate residuals for Day 1 after the sixth iteration.

Figure 5 .Figure 6 .
Figure 5.Time series of range rate residuals for Day 1 after the last (28th) iteration.

Figure 8 .
Figure 8. Amplitude spectral density of range rate residuals for Day 1 after the last (28th) iteration.