Gravity field models derived from the second degree radial derivatives of the GOCE mission: a case study

The GOCE satellite mission is one of the main achievements of the satellite geodesy for the Earth’s gravitational field recovery. Three different approaches have been developed for the estimation of harmonic coefficients from gradiometry data measured on board of GOCE-satellite. In this paper a special version of the space-wise method based on the second method of Neumann for fast determination of the harmonic coefficients Cnm, Snm of the Earth’s gravitational potential is given based on the radial gravity gradients of the EGG_TRF_2 product, except of two polar gaps filled by radial gradients from the EGM2008 gravity model. In the pre-processing stage GOCE-based second degree radial derivatives were averaged to the regular grid through Kalman static filter with additional Gaussean smoothing of residual radial derivatives. All computations are made by iterations. As the first step the determination of the preliminary NULP-01S model up to degree/order 220 derived from the Gaussean grid of the GOCE radial derivatives with respect to the WGS-84 reference field was developed based only one of the radial gradients EGG_TRF_2 in the EFRF-frame. In the second iteration the same algorithm is applied to build the NULP-02S gravity field model up to degree/order 250 using the same Gaussean grid with respect to the NULP-01S reference field. The NULP-02S model was verified by means of applying various approaches for the construction of the gridded gravity anomalies and geoid heights in the Black Sea area using processing of datasets from six altimetry satellite missions. Comparison of different models with GNSS-levelling data in the USA area demonstrates the independent verification of achieved accuracy of the constructed NULP-02S Earth’s gravity field model.


Introduction
Very often different applications in natural sciences may have the same mathematical tool. Tensor analysis represents one of such general approaches: "A particularly important role in applications (physics, engineering) is played by tensors represented by 3×3 symmetric matrices: strain tensor, Eötvös tensor (or Marussi tensor), the Earth's inertia tensor Müller 1987, Moritz 1990]" [Moritz and Hofmann-Wellenhof 1993]. It should be noted that the Eötvös tensor can transform to the Marussi tensor after dividing its components by gravity. Therefore the Marussi tensor can also be simply measured (calculated), taking into account the possibility of such observation (or computation with high resolution) of the magnitude g of the gravity vector g among various functionals of the gravitational potential V. In this paper (hereafter) the gravitational part of the Eötvös tensor consisting of the second derivates of the gravitational potential V [Seeber 2003] will be of special significance.
The Gravity field and steady state Ocean Circulation Explorer (GOCE) mission was launched in 2009 by the European Space Agency (ESA) in the nearly Sunsynchronous orbit at 96.5 degrees inclination and at an altitude of about 255 km with the GOCE onboard satellite gravitational gradiometer (SGG), measuring components of the gravitational tensor V. This leads to the increase of spatial resolution within certain frequencies and the corresponding improvement of the accuracy of the Earth's gravity field due to the measurements onboard of GOCE-satellite. Therefore, we get a significant upgrading of the middle and long wavelength components of V by means of data processing. Such experiment is the first satellite mission in which the 3D gradiometry in space was applied. Components of the tensor V of the potential V were measured worldwide excepting for both polar caps. So, nowadays one of the old Eötvös [1896] problems has a new development. According to a lot of publications [Schuh 1996, Seeber 2003, Koop et al. 2007, Bruinsma et al. 2009, Reguzzoni and Tselfes 2009, Schuh et al. 2010, Pail et al. 2011, Bruinsma et al. 2013, Brockmann et al. 2014, Förste et al. 2014, Gatti et al. 2014, Barzaghi et al. 2015 an expected improvement of the Earth's static gravity field models can be achieved up to 1-3 cm RMS in the accuracy of geoid heights with a half wavelength of about 80 km and longer. According to ICGEM (International Center for Global Gravity Field Models) information, the improved GOCE only gravity fields with the resolution up to degree/order 250-280 for the harmonic coefficients − C nm , − S nm of the gravitational potential V are usually obtained.
This study focuses on the algorithm for fast estimation of the harmonic coefficients − C nm , − S nm by means of the space-wise approach, which has been used for regular processing of GOCE data. The mathematical treatment of suggested methodology goes back to ideas of Gauss and Neumann [Sneeuw 1994] which can be applied to the GOCE Level 2 products. ESA provides different products of the GOCE mission. We will restrict this paper only by the study of the radial gravity gradient. GOCE derived gradients (EGG_TRF_2 files) are rotated [Rummel and Gruber 2012] from GRF to TRF system (terrestrial reference frame: north, west, up). These gradients are given in the local north-oriented frame (LNOF) through initial GRF system (EGG_NOM_2 files) having the exact connection with TRF where axis Z in GRF and LNOF almost equal. The X-axis has the north direction, the Z-axis is the perpendicular to this axis in the radial direction, and the Y-axis is cross-track completing a right-handed frame. Initial gradiometer data are measured in the gradiometer reference frame GRF and then transformed to the system LNOF via quaternion representation of rotation matrices [Migliaccio et al. 2004, Pail 2005, Bouman et al. 2009, Rummel and Gruber 2012. According to Bouman et al. [2015] the computation of potential derivatives in the local north-oriented frame (LNOF) requires the first and second derivatives of the gravitational potential V with respect to the polar coordinates given in the Earth fixed reference frame (EFRF) that leads to the relationship V zz = V rr between the GOCE gradients in LNOF and EFRF and to the following expressions for V and V rr : where i is the co-latitude, m is the longitude, r is the distance from the origin of the EFRF system to the current point P, GM is the product of the gravitational constant G and the planet's mass M, a is the semi-major axis of the ellipsoid of revolution, ( − C nm , − S nm ) are the fully normalized harmonic coefficients in the EFRF system where one axis points to the North pole, one to Greenwich and the third completes the system.
Because the second Neumann method [Sneeuw 1994] was chosen (based on the Gaussean grid for initial gradiometry data), contrary to other papers we will accept only one kind of data, in particular, second degree radial derivatives V rr , having a linear connection with unknowns − C nm , − S nm (see Equation 2). These derivatives V rr possess largest values among the tensor components in LNOF, high accuracy (about 1-10 mE; Müller and Wermut [2003], Bouman and Fuchs 2012, Bouman et al. [2009), and higher accuracy with respect to other gradients. Therefore, in this paper we treat only one component of gradient tensor as given data with a selection of a suitable fast method of data processing. According to [Sneeuw 1994, Sneeuw andBun 1996] the second Neumann's method, based on the Gauss' quadrature formula, includes a certain simplification. Mathematically speaking this contains "(a) special gridding assumptions, (b) band limitation that means no signal beyond the Nyquist frequency and (c) stochastic nature without noise". In reality we have usually average and corrupted by noise data [Sneeuw 1994]. There are different estimators that could be suitable for gradiometry data processing. The normalized coefficients − C nm , − S nm can be determined by one of numerical quadrature methods, which will provide essential accelerate computations. In contrast to the traditional approach of the least squares adjustment on the first iteration we come in this way to the intermediate NULP-01S gravity field model (up to degree 220) using WGS84 as reference field. In the second iteration the coefficients − C nm , − S nm of the NULP-02S gravity field model up to degree/order 260 were improved from the same radial gradients and NULP-01S, taken as reference field. It should be noted that the NULP-02S field is used in the following for the determination of the gravity anomalies and geoid heights in the Black Sea basin through processing of satellite altimetry data from six altimetry missions. Further verification of this and other gravity field models are made by means of comparisons with GNSS/levelling data in the USA area. Surprisingly, in all comparisons of NULP-02S and GOCE-only models a similar accuracy for the constructed gravity model, which was developed in this paper only from the radial derivatives V rr , is obtained.

Data pre-processing
Thus in this study the GOCE Level 2 data gradients as initial data are applied to get estimations of the global (2) gravity field. Selected data sets represent the second degree radial derivatives V rr which adopted equal to GOCE gradients in LNOF and EFRF coordinate systems as stated by [Koop et al. 2007, Bouman et al. 2015 for all three years of EGG_TRF_2 data (V rr with the polar coordinates i, m, r and accuracy) where the signal below MBW is replaced by the TIM R3 model gradients [Pail et al. 2011, Bouman et al. 2015. Further use of well-known methods (direct approach, spacewise or time-wise approaches) to data processing leads to the harmonic coefficients − C nm , − S nm of the gravitational potential.
In the first step, all the outlying values of the V rr data were deleted using the criterion 3d with respect to the EGM2008 model up to degree/order 360 [Pavlis et al. 2008]. Using datasets with the extension 0001.trf, about 1 million points were deleted; on the contrary only several hundred points were removed when datasets with the extension 0101.trf are analyzed. So, the datasets 0101.trf were applied for the determination of − C nm , − S nm . In this way the total amount of data is decreased insignificantly. GOCE data EGG_NOM_2 do not have a Gaussean white noise and are limited in MBW bandwidth usually from 5 mHz to 100 mHz in the GRF system [Müller et. al. 2010, Schuh et al. 2010, Polgár et al. 2013]. Due to error propagation the mixed GOCE and model data of the EGG_TRF_2 gradients V rr have lower accuracy in the LNOF system. Recent GOCE Level 2 products include also regular grids of the filtered gradients in the LNOF system limited in MBW bandwidth [Bouman et al. 2016]. However below we use another approach for V rr gridding. According to Novak et. al. [2012] the gradients V rr were reduced to the sphere with the mean radius of GOCE altitude applying a Taylor series expansion and the model TIM R3 for derivatives of higher orders. Since Level 2 product is obtained by processing Level 1 data, we assume Level 2 data partly filtered in MBW bandwidth. In the papers of Canuto and Massotti [2010], Sun et al. [2016] etc., the dynamical and extended Kalman filters were accepted for orbit determinations through gravity gradients. We use another approach for sequential calculation of the weighting average values at the centres of the spherical grid based on the static Kalman filter because EGG_TRF GOCE datasets are appeared also sequentially. According to Strang and Borre [1997], the "key to the Kalman filter is the idea of updating after reading new observations". Mathematically, a sequential least squares for a static problem are used: we express the new estimate of average value of the radial derivative V (k) rr as a linear combination of the old estimate V rr and new observations V rr (k) . It is nothing else but recursive weighted least squares method (3) where K k is kth gain Kalman matrix updated in this case via variance of the average value V (k) rr and the variance v 2 k of independent measurement. This algorithm was applied to the the EGG_TRF_2 datasets for different spherical cells and additional Gaussean smoothing of residual gradients dV rr for higher quality of solutions because of high frequency measurement errors. In this way Gaussean smoothing of residual gradients dV rr is considered as the high-pass filter. The spherical grid 15´×15´became sufficient for further use.
Since GOCE data have two polar gaps they were filled by V rr produced from EGM2008 up to 360 degree/order. So, all initial data were pre-processed and due to the measurements errors additional smoothing by Gauss filtering was applied to get rid of the low and middle frequency of the gradients V rr content. The building of Gaussean grid of V rr was based on the interpolation by modified local non-smooth splines [Marchenko et al. 2005], which were used to predict V rr from the spherical blocks 15´×15´at the nodes of Gaussean grid. It should be noted that the nodes of Gaussean grid connect with the following split along N circles of latitude (0 ≤ j ≤ 2N) with Dm = r/N ⇒ m j = Dm/2 + j ·Dm. Along the meridians the points are located at N parallels at the N zeros i i of the Legendre polynomial P N (cos i i ) = 0 of degree N. Therefore, the total number I of grid points is equal to I = 2N 2 .

The gravity field models derived from GOCE EGG_TRF_2 data
Then the difference between pre-processed V (k) rr and computed gradients (using a reference gravity model) at the Gaussean grid nodes was applied for the enhanced model construction by means of the removerestore technique. Consequently in the first step we get where according to (3) V rr is nothing else but the gridded V rr GOCE-based second order radial derivative, V M rr is the computed gradients from a gravity field model, dV rr -difference between measured and modelling gradients although the coefficients − C nm , − S nm connect with V rr in the linear form. In this way we come to the relationship (4) and the following equation coefficients of a selected reference gravity field model, d − C nm , d − S nm are the differences between coefficients of the new and old gravity field models. It is evident that the components V M rr can be calculated from the refer- The remove-restore procedure was applied for the evaluation of both of the global gravity models or iterations NULP-01S and NULP-02S. In the first step the second order gradients V M rr were computed using adopted WGS-84 normal field and according to Equation (6) V M rr were subtracted from initial gradients. After restoring through Equation (5) the NULP-01S gravity field model was constructed, with the harmonic coefficients up to 220 degree/order. In the second stage the gradients V M rr (produced by the NULP-01S model) were subtracted from measured radial gradients V rr . After restoring through Equation (5), Equations (8)-(10) the harmonic coefficients − C nm , − S nm of the NULP-02S model up to 260 degree/order were constructed.
Determination of the contribution into harmonic coefficients d − C nm , d − S nm was obtained from the second Neumann method [Sneeuw 1994] based on the Gauss quadrature formula written especially for the residual second order radial derivatives dV rr : where w i are the weights introduced for the discrete orthogonal functions system on the regular Gauss-Legendre grid, d pq is the Kronecker delta function.
Thus, in the first step using formulas (8)-(10) and remove-restore procedure the harmonic coefficients model NULP-01S up to degree/order 220 was developed based on the V rr gradients from GOCE and V M rr normal field adopted in this stage according to the WGS84 system. In the second step Equations (8)-(10) lead to the spherical harmonics model NULP-02S up to degree/order 260 based on the same V rr observed gradients and V M rr taken in accordance with the NULP-01S gravity field radial derivatives. Thus Equation (9) provides in both iterations the harmonic coefficients d − C nm , d − S nm or the contribution dV rr given only at satellite altitude.
Additional corrections were adopted according to Marchenko and Schwintzer [2003] by means of EGM2008 coefficients and the mean pole coordinates for the fixed coefficients − C 2m , − S 2m at epoch 2010 yr. Finally the NULP-02S model was truncated up to 250 degree/order. Geoids comparison was made with respect to the 1´× 1´EGM2008 geoid heights. Differences between geoids of NULP-02S (250,250) and EGM2008 up to 2190 degree/order are shown in Figure 1.
Thus Figure 1 illustrates in general a good agreement of the geoid surfaces NULP-02S and EGM2008 (with the mean deviation = −0.42 m and standard deviation = 0.51 m), excluding some scattered places observed in Antarctic, Greenland, Southern America, Himalaya, and Africa regions. This corresponds to fact presented in articles by several authors [Mayrhofer et al. 2010, Rummel and Gruber 2012 and indicates that GOCE data can improve the harmonic coefficients in the middle and long wavelength parts only.
It should be mentioned that initial data have a direct connection with the second degree radial derivatives V rr of gradient tensor (expressed in the EFRF polar coordinates) and the parameters of curvature of geoid surface. For better comparisons the influence of the centrifugal force was removed to avoid additional systematic effects. Therefore, the mean curvature J of the surface V = const based on the NULP-02S model (Figure 2) was computed according to Marchenko [2003] and compared with the mean curvature of the EGM2008 level surface V = const (Figure 3) reflecting a quality of both models. Figure 4 illustrates the models comparison in terms of spectral characteristics of geoid heights called as (PSA) "power signal amplitude" and "PSA differences" (see ICGEM, http://icgem.gfz-potsdam.de/ ICGEM/): The power variances of the geoid heights differences between two models are where D − C nm , D − S nm are the differences of spherical harmonics coefficients between two geopotential models, R is the mean Earth's radius (R=6371KM). The total    power variance of the geoid heights differences of two geopotential models is also used for the varying order The coefficients of the NULP-02S model starting from 200 to 250 degree/orders has increasing differences in power variances with respect to EGM2008 (Figure 4). At the same time these models have a good agreement up to 200-250 degrees. Small deviations about 1-3 cm are shown in Figure 4 up to degree/order 250, which may be caused by insufficiently smooth solutions related to processing of such satellite observations as GOCE data. These deviations correspond to small differences between geoid heights ( Figure 1) and irregularities of their curvature (Figure 2).

Validation of the NULP-02S model using GNSS/ levelling data
According to the seminal papers of Heiskanen and Moritz [1967], Krarup [1969], Tscherning [1976], Moritz [1980], Sansò et al. [1986], Sansò and Sideris [2013] etc., such linear geodetic functionals as geoid heights, free air anomalies, deflections of the vertical, and gravity disturbances can be used simultaneously for the approximation of the anomalous potential in the frame of the least squares collocation or Tikhonov regularization as operational approaches, providing processing of heterogeneous discrete data. It is evident that the NULP-02S model was built using only one type of homogeneous discrete data. To achieve an independent verification of NULP-02S with respect to the EGM2008 solution, we can compare GNSS/levelling geoid heights with the corresponding ones based on gravity field models. As stated by ICGEM a largest number of GNSS-derived geoid heights can be found for the USA area. Differences between observed and computed geoids were checked according to Milbert [1999] "Documentation for the GPS Benchmark Data Set" for 5379 GPS/levelling points (http://www.ngs.noaa.gov/GEOID/GPSonBM99/gps bm99.html). So, finally we demonstrate the quality of NULP-02S by means of comparisons with independent GPS/levelling points It should be noted that the determination of the normal heights H c or orthometric heights H ort using GNSS measurements requires subtracting the (quasi)geoid from geodetic heights (Equation 14). For further comparison, we use the basic relationships of  and in the first iteration the wellknown Bruns formula was applied. GOCE and gravimetric models were also tested on these GNSS-levelling data, which are consistent with other solutions. Generally this is no longer considered sufficient as new recommendations of ICGEM, such as Barthelmes [2013], propose that topography models are necessary for comparison of geoids. Thus, 5379 GPS-levelling points are compared with various models and resulting statistics are shown in Table 1. (13)

Validation of the NULP-02S model using satellite altimetry data in the Black Sea area
The method of satellite altimetry used in the last 20 years provides various Earth's sciences with information about a state of the ocean and its changes over time. With this satellite technology, the ocean surface is mapped in very simple way, which is based on altimetry satellite measurements of different missions with a level of accuracy 1-5 cm. Hence, recent altimetry satellites give more and more accurate determination of the sea surface heights (SSH), which in combination with an ocean circulation model can be considered as a direct determination of the geoid undulations in different ocean areas.
On the other hand, appropriate data sets from various altimetry satellites are also collected for such relatively small internal marine basins of the globe such as the Black Sea and Sea of Azov. As a consequence, the inversion of SSH data set into the gravity anomalies Dg becomes possible within the mentioned internal marine areas, assuming filtered SSH data as "measured" geoid undulations N. The NULP-02S gravitational field was applied in the frame of the remove/restore procedure for additional verification of this model in the Black Sea basin. The methods used are: least-squares collocation [Moritz 1980, Tscherning 1994, Kilicoglu 2005; Tikhonov regularization [Moritz 1980, Marchenko et al. 2001; FFT solution (fast Fourier transformation) in the studied area as stated by Schwarz et al. [1990].
According to Marchenko and Yarema [2006], our data set 1 represents 643128 TOPEX/POSEIDON, ERS-1, ERS-2, JASON-1, ENVISAT and GFO SSH (sea surface heights) points taken in the time-period 1992-2005 and corrected by CLS AVISO for different geophysical phenomena and instrumental effects (declared accuracy ~2-5 cm). Because of data gaps in the corrected SSH, the additional set of point gravimetry data was used to support the SSH-only solution for Dg. So, our data set 2 represents 20207 values of BGI and national marine point gravimetry measurements in addi-tion to land gravimetry surrounding the Black Sea area (with mean accuracy about ~5-8 mGal).
In this study two types of the empirical covariance functions (ECF) are used: in the first step the ECF of residual geoid heights with possible systematic parts when satellite altimetry data are inverted into residual gravity anomalies by means of adopted method. In this step we come to residual systematic effects since AVISO SSH data are used before in the time-period of 15 years with various time-dependent models for various reductions. Then, the regularization parameter is accepted as equal to one (=1) and applied for the simultaneous processing the improved by fitting covariance function of residual gravity anomalies, including their cross-covariance function by means of collocation method only taking into account the covariance propagation. Mathematical explanation could be approved, if a statistical treatment of the least squares collocation is not considered. According to Moritz [1980] this method is nothing else but the method of the stable solution of the ill-posed problem. Hence, a stability of solution can be achieved via introduction of a kernel function which defines a certain Hilbert H q 2 (∑) (Sobolev) space of harmonic functions, i.e. a possible set of solutions for the anomalous potential T. Sobolev spaces H q 2 (∑) are complete subspaces of the general Hilbert space and therefore Hilbert spaces by themselves. They are well-known in physical geodesy as Hilbert spaces of harmonic function with reproducing kernels which have the evident classification Lelgemann 1998, Marchenko 1998]. Figure 5 illustrates AVISO corrected SSH (referred to the GRS80 ellipsoid) from six satellite missions in the Black Sea area, which are filtered by Gauss smoothing procedure for residual temporal effects. It has to be pointed out, that the surface in Figure 5 differs from the geoid due to small deviations of sea surface topography (SST) and various residual effects averaged in time at the level about 10 cm [Marchenko and Yarema 2006]. Figure 6 demonstrates the result of recovery of the gravity anomalies Dg from the combination of the above two datasets into the 2´× 2´grid of Dg. In spite of different gravity field models (EGM96, EIGEN-CG01C, EGM2008, NULP-02S) applied in these tests, all methods give a good accordance for geoid solutions in terms of the standard deviation within 10 cm which coincides with estimated level of SST in Black Sea basin.
The reproducing kernels (analytical covariance functions) applied in both cases were described by singular point harmonic functions [Marchenko 1998]. The optimal degree n = 1 was chosen that corresponds to the modified Poisson kernel or the so-called dipole kernel. Observed and predicted gravity anomalies were compared at control-scattered points over the Black Sea area that gives a similar agreement of the collocation, regularization, and FFT solutions in terms of a stan-dard deviation of about 8 mGal of the difference "observed -computed" by various methods. Nevertheless, in all comparisons we should note certain disagree-  ments in the areas near seashores, which is a weak point of the solution in Figure 6. Similar results of the comparison were obtained for the observed Dg in this area with the standard deviation about 5-8 mGal in the "East and West sub-regions" [Kilicoglu 2005], based on the direct using of the well-known GRAVSOFT package [Tscherning 1994] and the kernel function of line singularities. The latter has a certain theoretical significance since the modified Poisson kernel (dipole kernel) or reproducing kernel of point singularities was used successfully in the frame of such instable problem as the inversion of SSH data into gravity anomalies.

Conclusions
The Gauss' quadrature formula makes available a simplest approach for fast computations of the harmonic coefficients of the global Earth's physical fields without degradation of accuracy. Hence, contrary to previous studies the proposed version of the space-wise approach is based on the second Neumann method (allowing a fast estimate of the harmonic coefficients of the gravitational potential) using radial gravity gradients of the EGG_TRF_2 datasets, excepting the polar regions filled by radial gradients from the EGM2008 gravity model up to 360 degree/order. All computations are made by iterations. In summary we would like to emphasise several important aspects of presented solutions.
-In the pre-processing stage GOCE-based EGG_ TRF_2 second degree radial derivatives were averaged sequentially at the regular grid through Kalman static filter (3) with additional Gaussean smoothing of the residual radial derivatives dV rr . This crucial stage leads to the Gaussean grid filled by the radial derivatives dV rr and permits the corresponding version of space-wise approach via Gauss' quadrature formulas (8-10).
-The resulting NULP-02S (as second iteration) global gravity model provides a good agreement with the independent GNSS-levelling data in the USA and other areas with the standard deviation about 50 cm.
-Satellite altimetry (SSH) or/and gravimetry data can be used for recovery of the shortwave part of the geopotential either using GOCE models or the NULP-02S solution. The regularization method, the leastsquares collocation, and FFT technique give in the Black sea basin a similar results in terms of the standard deviation about 10 cm. In all cases the inversion of SSH data into gravity anomalies is required for further comparison with independent data of control-scattered marine gravimetry.
-The gravity model NULP-02S referred to epoch 2010 demonstrates a similar accordance in terms of the mean and standard deviations with other GOCE-only models. This fact illustrates that 2nd order derivatives from EGG_TRF_2 data can improve the middle and long waves of the geopotential. At first sight these radial derivatives alone provide the simplest approach of the estimations of the global gravity models. Nevertheless, to confirm such a possibility we need to repeat a data processing with the full gradient tensor in the GRF system which requires a change of the adopted algorithm for the determination of − C nm , − S nm .