Modeling the Recovery of the Earth’s Gravitational Field from Satellite Measurements Using Parallel Computations

The


Introduction
Depending on the granularity of the representation and measurement methods of the Earth's gravity field (EGF), three components of the gravity field can be distinguished: local, regional and global.The local level of the EGF is provided by land and sea gravimetric surveys, as well as by airborne gravimetry.The regional component is represented by satellite altimetry of the Earth's surface and satellite laser ranging.The global gravity field is measured by satellite gravimetry.Of course, such division of the gravity field into three components is not strictly defined.
While the CHAMP and GOCE missions were single spacecraft, the GRACE and GRACE Follow On missions consist of a pair of identical spacecraft in the same near-polar Earth orbit, following each other at a distance of about 200 km.The main measured quantity in such satellite configuration is the inter-satellite range, which varies depending on the gravitational field of the terrain (oceans, mountains, etc.) over which the satellites of the constellation are currently flying.The measurements of the inter-satellite range provide the necessary information for the solution of the inverse problem of recovering the Earth gravity field.The main instrument that measures the inter-satellite range in the GRACE and GRACE-FO missions is a high-precision microwave interferometric system, which operates at a micrometer accuracy level [21,22].In parallel with the microwave system, the GRACE-FO mission satellites are equipped with a laser interferometric system, which allows to raise the measurement accuracy to a nanometer level [1].
Twin satellite configurations have proven to be highly effective in providing high-resolution global EGF models, so there is a high probability of launching new missions of that type in the future.In particular, as a part of the Chinese TianQin space gravitational wave detector project [12], there are plans to launch the TianQin 2 constellation, similar to GRACE-FO, at the end of 2025 [25].
With the beginning of data collecting from near-Earth gravity missions, static global models of the EGF began to be developed.In addition, data from the GRACE and GRACE-FO missions allowed the production of monthly models, i.e., with a temporal resolution of one month.A complete and regularly updated list of such models is published on the website of the International Center for Global Earth Models (ICGEM) [8].EGF models are represented as a series expansion of spherical harmonics, whose degrees and orders characterize the spatial resolution of the model, or, in other words, the detalization of the EGF representation.The first static EGF model produced using the CHAMP mission data, the EIGEN-1 model, was published in 2002 [19].EGF models can be based either on data from a single mission or on a combination of data from different missions.As an example, one of the most recent (2023) EGF models is the WHU-SWPU-GOGR2022S model built from GOCE and GRACE data; this model is complete to degrees and orders of 300 [24].
There are static models of the EGF built from data from the GRACE missions only.The first such models are the ITG-Grace03 [14] and ITG-Grace2010s [15].Both models are represented as a decomposition up to 180 spherical harmonics.The most recent model, which uses all GRACE observations since April 2002, is the ITSG-Grace2018s model [10,16].
Monthly models of the EGF built from GRACE and GRACE-FO data are created in three international centers: Center for Space Research at University of Texas, Austin (CSR), Helmholtz Centre Potsdam German Research Centre for Geosciences (GFZ), and Jet Propulsion Laboratory (JPL).The solutions obtained in these centers are also available on the ICGEM website.
Designing a similar Russian space mission to measure the Earth's gravitational field is also being considered.Some research on this subject is being carried out in a number of scientific and industrial organizations of the country, not only on the development of measurement systems and scientific background, but also on the modeling of the EGF recovery from satellite measurements.The latter is a rather resource-consuming computational process.The number of unknown values in the recovery of the EGF parameters depends quadratically on the degree of decomposition.Therefore, when recovering the EGF up to degree N max = 200, it will be necessary to operate with matrices of 40000 × 40000 and more.This requires hundreds of gigabytes of RAM to store such amount of data if double precision numbers are being used.In this regard, parallel computation in the EGF recovery problem plays a crucial role in optimizing the computational process.This paper presents the results of such modeling using parallel computing based on real data from the GRACE Follow On mission.
The paper is organized as follows.Section 1 presents the mathematical model of the gravity field recovery from satellite measurements.Section 2 describes the parallelization algorithm of the computational process.The results of the parallelization are given in Section 3. Section 4 contains the results of the EGF recovery.The conclusion summarizes the study.

Mathematical Model
The GRACE and GRACE-FO missions have provided the scientific community with a significant amount of inter-satellite measurements.In order to use this data to recover the EGF, it is crucial to have prior knowledge of the exact positions of the spacecraft in orbit, as well as their relative locations.The main objective is to assess the impact of all potential effects, including both gravitational and non-gravitational forces, on the motion of the spacecraft.Therefore a mathematical model of spacecraft motion has been developed to numerically integrate the equations of their orbital motion.
The main forces considered in modeling such motion include the non-sphericity of the geopotential, the attraction of the Moon, Sun and planets, solid tides in the Earth's body, ocean tides, resistance of the atmosphere, and solar radiation pressure.Each of these forces can cause on a GRACE-type spacecraft an acceleration greater than 10 −9 m/s 2 .Accelerations below this value will be compensated for by a non-gravitational acceleration compensation system and can be disregarded.It is worth noting that, when forming a complete mathematical model of inter-satellite measurements, it is necessary to take into account the relativistic effects of signal passage from one spacecraft to another.
As mentioned in the introduction, recovering the EGF is a resource-consuming procedure.Therefore, the idea of using parallel computations in structures that can be parallelized arises naturally.Examples of such structures include the non-sphericity of the EGF and ocean tides.
The nonsphericity of the Earth's gravitational potential is the main factor to be considered when modeling spacecraft motion with high precision in low-Earth orbit.Mathematically the non-sphericity of the EGF is determined by the Stokes coefficients.These coefficients form the base of its decomposition by Legendre polynomials.In this work, we utilize the Belikov and Taibatorov algorithm for decomposing the EGF [5].The main advantage of the algorithm is its novel normalization of Legendre polynomials, which eliminates the need to recalculate coefficients during the recurrent procedure.
The expansion of the EGF in this algorithm is achieved using the following expressions: where U is the Earths gravitational potential; GM is the Earths gravitational parameter; R is the mean equatorial radius of the Earth; r is the geocentric distance to a given point; θ is the polar angle (0 ≤ θ ≤ π); λ is the longitude (0 ≤ λ < 2π); n and m are degree and order of the EGF expansion; C nm and S nm are dimensionless harmonic coefficients (the Stokes coefficients); N max is the maximal degree of the expansion.The coordinate derivatives of the gravitational potential give the acceleration of the satellite along the corresponding axes as functions of the new variables T m n : which simplifies the calculation of the accelerations.
The following are the expressions for the derivatives of geopotential with respect to spherical coordinates: The Cartesian derivatives of the EGF, which are used in the numerical integration of the satellites equations of motion, are obtained from (3) by the following relations: It is also worth noting that the Belikov and Taibatorov algorithm is a recurrent "column" type algorithm.This means that the right-hand side contains functions with the same index n − 1, as opposed to a "row" type algorithm.This property makes the algorithm fast and stable even for calculations with large values of n and allows it to be efficiently parallelized, since the summation limit N max can reach the values N max = 200.
The expression used to calculate the impact of ocean tides also employs double summation.The corrections to the Stokes coefficients take into account the effect of ocean tides in the same way as for solid tides.The formula for calculating these corrections is as follows: where C ± f,nm and S ± f,nm are the harmonic amplitudes of the geopotential for the tidal component at frequency f [17].These amplitudes are model-dependent.We used the FES2014 ocean tide model [4,6,13], which considers degrees and orders up to 100 for all tidal waves, except for long-period tidal waves, which are calculated up to and including the 50 th degree.
Finally, the third part, which can be significantly optimized by introducing parallel computations, is directly related to the procedure of the EGF recovery performed with the least-squares method.During the recovery process, all observations made on the spacecraft are divided into small 30-minute intervals (arcs), and the same observation accuracy is assumed for each arc.Thus, the following normal equations can be formed for each arc: with and where the components are given by } is the radius-vector of each satellite in the GCRS inertial frame; r1,2 = {a x1,2 , a y1,2 , a z1,2 } is the acceleration vector of each spacecraft; is the distance between spacecrafts; ω i is the observation weights, such that √ ω i = 1/σ i , where σ i is the standard deviation of the observations for the i-th arc; ∆ξ i = rO i − rC i is the difference between the observed and modeled inter-satellite accelerations.
The normal equations for each arc (with a total duration of about 30 days for a monthly model) can be combined into a general system: where and I is the total number of arcs.This combination is done under the assumption that the covariance matrices of each group of observations are independent.The value σ k is the measurement error of the observation group k.The solution of the system can be easily found in the form Therefore, the calculation and summation can be performed in parallel for individual arcs.

Parallelization Algorithm
The mathematical model of the EGF recovery, as described above, includes models for calculating geopotential non-sphericity and ocean tides.Measurements of the execution times of the program modules corresponding to these models have shown that these two algorithms are the most time-consuming.Since the EGF recovery algorithm is programmed in Fortran language, it was possible to use OpenMP parallelization tools [7] with the help of the PARALLEL DO/END PARALLEL DO construct, which was utilized to compute the ocean tide corrections (Fig. 1).

Figure 1. Parallelization of the part of the code responsible for calculating ocean tides
A more complex approach was used in the subroutine responsible for calculating the geopotential.In addition to the PARALLEL DO construct, there is a part of the code that adds up the values obtained in different execution threads, and this fragment is executed sequentially by the CRITICAL/END CRITICAL construct.The calculation of each type of derivative (dV_dr, dV_dteta, dV_dlyam) is done in its own execution thread (Fig. 2).
The partial sums dV_dr, dV_dteta, dV_dlyam are obtained separately in each thread (CRITICAL section) and recorded in the variables dV_dr2, dV_dteta2, dV_dlyam2.
The parallelization of the EGF recovery process is performed as follows: matrices with the calculated coefficients can be written to disk for each arc, or they can be summed up during the calculation and written to disk once at the end of the process.From a runtime point of view, it is better to sum the matrices at the end of the calculation, but this requires more memory.The program allows efficient parallelization, since the calculations for each arc are done independently.The program divides the whole range of arcs into N parts and starts N processes simultaneously for each part.After all the parts of the time interval have been processed and Figure 2. Parallelization of the part of the code responsible for the geopotential calculation the coefficient files have been written (this is the 1 st processing stage), the procedure for reading the files with the partial coefficient matrices and their summation in RAM is started (the 2 nd processing stage).
Previously calculated files are used to recover the EGF coefficients.This is done by the least-squares method (LSM) using the routines of the LAPACK library [3].After each run of the LSM, a new file containing the EGF coefficients is written.The coefficients are written up to the degree N max specified in the initialization file.
The parallelization of the EGF recovery process, implemented in the Fortran programming language, is universal for the following three cases: • computing cluster; • multi-core computer running Linux OS; • multi-core computer running Windows OS.The differences between the systems are localized in a single procedure that starts a child process by calling "system()".Different values of environment variables are used for the "computing cluster" and the Windows OS environments, in other cases it is assumed that the program is started under the Linux OS.
Synchronization of parallel computations is performed by creating/deleting flag files (the process of the next computation stage waits for the end of the previous stage).The desired number of parallel processes is always selected at startup using command line arguments.The number of parallel processes depends on the available RAM, because each process running at the 1st stage needs a certain amount of RAM to work, and all these processes are running simultaneously.For example, if the recovery is performed up to N max = 300, then each process of the 1 st stage would require 120 GB of RAM.If the computer only has 512 GB of RAM, then the stage can only be divided into 4 parallel processes.

Results of Parallelization
The performance of parallelization on computing nodes was tested by measuring the time T (1) required to compute on one node and the time T (N ) required to compute on N nodes: Here, the SPD parameter indicates how much faster the parallelized computation is compared to the single node variant.The Eff parameter is a measure of the parallelization efficiency and shows the impact of the overhead on the program execution time.The tests were performed on nodes each having an Intel Xeon E5-2680 v4 2.40GHz/3.3GHzprocessor with 28 cores and 128GB RAM.
Table 1 shows the speed-up of the program complex as a function of N (the number of computing nodes used).It should be noted that the results are influenced by the speed of reading and writing intermediate files, which can vary considerably.The initialization procedures cannot be parallelized and start to take up a larger and larger part of the total computation time with increasing N .This explains the decrease in efficiency with increasing number of processes.When N > 16, the time for reading intermediate files increases drastically.The decrease in efficiency for large N , which has the same origin, can also be observed in Tab. 2.
Table 1.Computation performance of the EGF recovery on the interval of 1 month, up to the maximal degree of the expansion N max = 96 (1) 3 summarizes the computation time of the EGF recovery for a 60-day interval.It shows that the number of processes executing the 1 st stage decreases as N max increases, which is due to the RAM limit.The values are based on a specific computer configuration (Intel Xeon Gold 6148 x 2, 512GB RAM) and are not universal.

Results of the EGF Recovery
The developed program allows to create models of the EGF with different temporal resolutions using real measurements of the inter-satellite range of the GRACE and GRACE-FO missions.These models can cover long time intervals, such as a year or more (static models), as well as short time intervals, such as monthly models.The latter can track global processes that change over time, such as mass movements in the atmosphere and oceans, melting of glaciers, and more.
Various methods can be used to graphically present the results of the EGF recovery.One such method is visualizing the distribution of the gravity field over the Earth's surface.In this case, the EGF is represented either in acceleration units or in geoid heights.Milligals (mGal) are typically used as the acceleration units, defined as 1 Gal = 1 cm/c 2 .The EGF expressed in geoid heights is the difference between a generalized surface of the Earth defined by the actual gravitational field and the ellipsoid of rotation defined by the C 20 coefficient in the expansion of the gravitational potential (1).Accordingly, the height of the geoid is determined by the where P nm (sin ϕ) are the associated Legendre polynomials of degree n and order m; ϕ, λ are the angular spherical coordinates of a given point at the Earths surface; R is the mean radius of the Earth; C nm and S nm are the dimensionless Stokes coefficients.Geoid heights are usually expressed in meters.
The following presents the results of the Earth's gravity field recovery using real data from the GRACE-FO mission.These results were obtained using the program described in this paper.
Figure 3 shows the static model of the Earth's gravitational field MSU2024-1 built from GRACE-FO data for the whole year 2021.The model is decomposed to degrees and orders of n = m = 120 and presented here in terms of geoid heights.
It represents the difference between two EGF models per degree, i.e., for each degree, the orders are summed.It is given in terms of geoid heights.The quantities ∆C nm and ∆S nm in (15) represent the difference between the Stokes coefficients of the reference (REF) and recovered (REC) models: ∆C nm = (C nm ) ref − (C nm ) rec and ∆S nm = (S nm ) ref − (S nm ) rec .Figure 4 shows a comparison between the EGF obtained from GRACE-FO data in this study and the reference models of CSR, GFZ, and JPL.As mentioned in the introduction, monthly temporal models are created in three international centers: the Center for Space Research at the University of Texas, Austin (CSR), the Helmholtz Centre Potsdam German Research Centre for Geosciences (GFZ), and the Jet Propulsion Laboratory (JPL).Figure 4 displays the 'gravity' curve for the cumulative gravity signal per degree (in terms of geoid heights) of the reference EGF, excluding the second zonal harmonic C 20 .The degree difference variance between the recovered and reference models shows a characteristic behavior of the errors of the spherical harmonic coefficients.Specifically, the differences decrease and reach a minimum value in the low-frequency and early mid-frequency parts of the spectrum (as seen, for example, in [23]), before increasing again.These differences result in minimum geoid heights of about 0.1 to 0.2 mm for degrees n = 30 to 50.

Conclusion
Global EGF models, built from data collected by space geodetic missions, play a very important role in studying global processes across Earth's various geospheres.This paper presents the results of the EGF recovery using parallel computing based on real data from the GRACE Follow On mission.As degree of decomposition N max increases, the recovery of EGF parameters becomes more time consuming and expensive in terms of memory requirements.This is because the number of parameters increases quadratically with N max .For example, recovering the EGF up to a large degree N max ≈ 200 requires hundreds of gigabytes of RAM.Therefore, parallel computations are crucial for optimizing the computational process for the EGF recovery.
Parallelization is especially effective for code sections that involve summation.In this study, we utilized OpenMP parallelization tools to compute the geopotential non-sphericity and ocean tides models.Furthermore, we applied parallelization directly to the EGF recovery procedure, which can be significantly accelerated by dividing the calculation of the normal equation matrix into parts and then summing them up.
Testing showed that parallelizing the program for recovering the EGF over a one-month interval with N max = 96 resulted in an 11-fold reduction in computation time.Similarly, recovering over a one-year interval with N max = 180 resulted in a 22-fold improvement.It was found that the speed of access to the cluster's common filesystem was critical to recovering the EGF with large N max .The optimal number of nodes for simultaneous processing strongly depends on the speed of reading files with matrix coefficients.Other factors determining the optimal number are the maximum degree N max and the duration of the time interval to be processed.
With the help of the developed program the static model of the Earth's gravitational field MSU2024-1 was built using the GRACE-FO mission data for the whole year 2021.The model is decomposed to degrees and orders of n = m = 120 and presented in terms of geoid heights.We also compared the EGF recovery on a monthly interval using the GRACE-FO data obtained in this work with the CSR, GFZ, and JPL temporal models built at other world centers.
Further research will focus on software modifications to increase the spatial and temporal resolution of EGF models.

Figure 3 .
Figure 3. Results of the EGF recovery up to n = m = 120 using the GRACE-FO data for a one-year interval (the static model MSU2024-1).Colors correspond to geoid heights in meters

Table 2 .
Computation performance of the EGF recovery on the interval of 1 month, up to the maximal degree of the expansion N max = 180 (1), time in seconds

Table 3 .
The computation time of the EGF recovery with parallelization for a 60-day interval N max Number of threads 1 st stage (hrs) 2 nd stage (hrs) Total time (hrs)