GRACE gravity field recovery

This cookbook chapter describes an example of estimating a gravity field solution using GRACE observation data. For the respective month a set of spherical harmonic coefficients up to a maximum degree is determined. An example scenario for this task can be found at https://ftp.tugraz.at/outgoing/ITSG/groops/scenario/scenarioGraceGravityfieldRecovery.zip including the required GROOPS scripts and data sets for the gravity field recovery process. The background models are provided at https://ftp.tugraz.at/outgoing/ITSG/groops/data/.

Background models

The following background models were used during the data processing:

These models were reduced during the analysis process and are not present in the solution. The GOCO06s model was used as the static gravity field as well as for the trend component and annual oscillation. In the script 000groopsBackgroundModels.xml a monthly mean of the GOCO06s including the time-variable components is determined in form of time splines using Gravityfield2TimeSplines. This model is later added back to the final gravity solution.

Instrument data preparation

The ITSG gravity field solutions are computed from the official GRACE L1B JPL (2018) and GRACE-FO L1B JPL (2019) observation data. The data sets for this example are provided in GROOPS file format in the scenario folder.

The satellite-to-satellite-tracking (SST) data consists of:

Additional observation data required for the processing comprises:

The determination of

is depicted in Kinematic orbit determination of LEO satellites. These data sets are also provided in the scenario folder.

Data preparation is handled in the script 010groopsInstruments.xml. The approximate orbits (initial dynamic orbits) of the satellites, the star camera observations, the accelerometer data and the thruster data are resampled with a 5s sampling and small gaps in the data are filled using InstrumentResample. Gross outliers are removed using InstrumentRemoveEpochsByCriteria and the data is synchronized using InstrumentSynchronize.

The approximate orbits are later used as a priori information for the dynamic orbit integration. In addition to the observed orientation of the spacecrafts (star camera observations), the nominal orientation is computed using SimulateStarCameraGrace. The difference between observed and simulated orientation is determined using InstrumentStarCameraMultiply and is employed in the outlier detection.

The accelerometer data is initially calibrated by estimating a bias using InstrumentAccelerometerEstimateBiasScale with respect to simulated data created with SimulateAccelerometer. For simulating accelerometer data a satellite model implying the satellite's mass and surfaces is required. Such a model can be created with SatelliteModelCreate. Models for the GRACE and GRACE-FO satellites are also provided at https://ftp.tugraz.at/outgoing/ITSG/groops/data/satelliteModel/. Non-gravitational forces comprising atmospheric drag, solar radiation pressure and albedo have to modeled when simulating the accelerometer data. The acceleration bias parameters are determined as degree 3 time splines with 6h nodes. When determining these parameters the thruster events are excluded from the estimation.

The SST observations, the light time corrections and the antenna center corrections are synchronized with a 5s sampling together with simulated SST data created with SimulateSatelliteTracking. Simulated data is used for the outlier detection of the original SST observations.

The sampling of the kinematic orbits is reduced to 60s using InstrumentReduceSampling and an outlier detection is performed using the approximate dynamic orbits.

The approximate orbits, the star camera observations and the accelerometer data are divided into 24h arcs (variational arcs). The kinematic orbits, its 3x3 epoch covariances, KBR observations, light time corrections, antenna center corrections and star camera observations are divided into 3h arcs per day (short arcs). Additionally the approximate orbits and the star camera observations are also synchronized to short arcs.

Further information on instrument data preparation can be found in Instrument data handling.

Variational equations

In this processing step dynamic orbits are computed for a complete 24h orbit arc by integrating the forces acting on the GRACE/GRACE-FO satellites. Additionally, the state transition matrix is set up. The dynamic orbits are then fitted to kinematic orbits and SST observations in a least squares adjustment by co-estimating additional accelerometer calibration parameters together with the initial state vector. The newly estimated parameters are then used to re-estimate the dynamic orbits and setting up the new state transition matrix.

The script 020groopsVariational.xml in the scenario folder implements the required processing steps. Time splines from a time-variable gravity field are estimated using Gravityfield2TimeSplines. In this step the static gravity field (GOCO06s) is combined with the following time-variable components:

maxDegree=220 and sampling=10/1440 is sufficient.

In PreprocessingVariationalEquation the variational equations comprising the integrated orbit together with the state transition matrix are stored in outputfileVariational.

This program has to be executed for both GRACE or GRACE-FO satellites and it is recommended to use LoopPrograms.

The force models include:

In PreprocessingVariationalEquationOrbitFit the integrated orbit (inputfileVariational) is fitted to the kinematic orbit (inputfileOrbit) by least squares adjustment. The additional accelerometer calibration parameters can be defined by

The observation equations (parameter sensitivity matrix) are computed by integration of the variational equations (inputfileVariational) using a polynomial with integrationDegree=7. PreprocessingVariationalEquationOrbitFit has to be executed per satellite.

PreprocessingVariationalEquationSstFit fits two dynamic orbits inputfileVariational1/2 to the SST observations and the kinematic orbits.

The estimated accelerometer calibration parameters from PreprocessingVariationalEquationOrbitFit and PreprocessingVariationalEquationSstFit are determined as corrections and stored in outputfileSolution. Both correction estimates have to be summed up using FunctionsCalculate.

The dynamic orbit and the resulting accelerometer calibration parameters are now used to re-integrate the orbit once more using PreprocessingVariationalEquation and introducing parametrizationAcceleration as estimatedParameters. This step usually ensures convergence. If the maximum orbit difference is still not sufficient this step can be repeated again.

Preprocessing

The script 030groopsPreprocessing.xml implements the following steps and settings. The program PreprocessingSst processes SST observations and kinematic orbit data and performs a complete least squares adjustment for gravity field determination by computing the observations equations. Force model parameters (gravitational potential coefficients and accelerometer calibration parameters) are computed by integrating the parameter sensitivity matrix from the variational equations. Parameters describing effects due to the SST observation system and geometry (KBR antenna phase center variations) are computed using the dynamic orbits as a Taylor point. Short time gravity variations can be co-estimated together with the monthly mean gravity field. The autoregressive model sequence constraining the short time parameters is provided in the data folder. See Kvas 2019 for more information about this co-estimation.

ParameterSelection2IndexVector and MatrixCalculate with matrix:reorder can be used to extract the desired spherical harmonic coefficients from outputfileSolution and the respective standard deviations from outputfileSigmax up to a certain degree.

In the program Gravityfield2PotentialCoefficients the estimated spherical harmonics coefficients are read with gravityfield:fromParametrization. The monthly mean gravity field can be added back by additionaly selecting the time splines created in 000groopsBackgroundModels.xml using gravityfield:timeSplines. The preprocessing solution is saved as a spherical harmonics file.

Setting up normal equations

Normal equations are set up in the script 040groopsMonthlyNormals120.xml using the program NormalsBuildShortTimeStaticLongTime. The time intervals which the normal equations are divided into are defined in inputfileArcList. The normal equations are based on observation including the SST data, the kinematic orbits and the variational equations. The parametrization of the gravity field can be set with observation:parametrizationGravity (e.g. spherical harmonics up to degree and order 120). Accelerometer calibration parameters and KBR antenna phase center variations can be parameterized using parametrizationAcceleration and parametrizationSst. With estimateShortTimeVariations short time variations of the gravity field can be co-estimated. The parameters selected by parameterSelection (e.g. linear splines with 6h nodes) are constrained by an autoregressiveModelSequence. Additional temporal variations (e.g. trend and annual oscillation) could be estimated with estimateLongTimeVariations.

Solving normal equations

The desired spherical harmonic coefficients are determined in the script 050groopsMonthlySolve.xml. NormalsSolverVCE accumulates normalEquation and solves the total combined system. Variance component estimation is used to determine the relative weighting of the individual normals. The estimated parameter vector (outputfileSolution) the estimated accuracies (outputfileSigmax) and the full covariance matrix (outputfileCovariance) can be saved. Using Gravityfield2PotentialCoefficients the final solution can be saved as a spherical harmonics file by adding back the monthly mean gravity field to the estimated spherical harmonic coefficients.