Analysis of high - resolution global gravity ﬁ eld models for the estimation of International Height Reference System ( IHRS ) coordinates in Argentina

: Following the de ﬁ nition and realization of the International Height Reference System ( IHRS ) , the ver - tical coordinate of a given point at the Earth ’ s surface can be obtained from the computation of the geopotential value from a harmonic expansion of a Global Gravity Model of High - Resolution ( GGM - HR ) or based on the computation of a local or regional pure gravimetric geoid or quasigeoid. Therefore, an evaluation of the accuracy of GGMs - HR and the geoid model available is needed in order to assess its capability to infer IHRS coordinates. In this study, di ﬀ erent GGMs - HR are evaluated against 2287 benchmarks in Argentina. Moreover, the most recent geoid model of Argentina is also evaluated. Geoid undula - tions at these benchmarks are obtained based on ellipsoidal and orthometric heights in the local vertical datum. Results suggest that among the evaluated GGMs - HR, XGM2019e provides the best agreement with the observed geoid heights, but none of them is accurate enough in order to infer vertical coordinates in the IHRS. Similar conclusions are obtained for the local geoid model for Argentina demonstrating the necessity for a more precise geoid model, following the standards and recommendations given for the IHRS.


Introduction
With the definition of the International Height Reference System (IHRS) given by the International Association of Geodesy (IAG) in 2015, many efforts are being made for the realization of the IHRS, the International Height Reference Frame (IHRF).The definition of the IHRS is based on the following conventions: the vertical coordinate of a point P, close to or on the Earth's surface, is given by the geopotential number Cartesian coordinates of P refer to the International Terrestrial Reference Frame (ITRF).W(P) is the potential of the Earth's gravity field at a point P (Drewes et al. 2016).
According to Sánchez et al. (2021), the determination of physical coordinates worldwide referring to the IHRS can be obtained by (i) the direct computation of the geopotential value at any point by introducing the ITRF coordinates into the harmonic expansion equation representing a Global Gravity Model of High-Resolution (GGM-HR), or (ii) solving a Geodetic Boundary Value Problem (GBVP) for regional pure gravimetric geoid (or quasigeoid) determination in the remove-compute-restore (RCR) scheme.In both cases, a GGM has to be included, and therefore, its reliability has to be assessed by, e.g., comparison of geoid heights inferred from the GGM with those derived from levelling co-located with GPS at given benchmarks.
The main objective of this study is to present a qualitative analysis of the latest combined GGMs-HR and the current official geoid model of Argentina, GEOIDE-Ar16 (Piñón 2016), using GPS/levelling benchmarks over the Argentinean territory.For the analysis, 2287 GPS/levelling benchmarks that cover the mainland part of Argentina have been taken into account.This analysis needs a careful treatment of the permanent tide in each involved quantity.To consistently compare geoid undulations, all observations and data should be related to the mean tidal concept/mean crust (see IAG Resolution No. 1 (2015) in Drewes et al. 2016).The zero-degree term for the geoid was carefully analysed using the GRS80 reference ellipsoid (Moritz 2000) and the conventional reference value W .

IHRS
The tested GGMs-HR include: EGM2008 (Pavlis et al. 2012), EIGEN-6C4 (Förste et al. 2015), GECO (Gilardoni et al. 2016), SGG-UGM-2 (Liang et al. 2020), and XGM2019e up to degree and order 2159 (Zingerle et al. 2020).The dedicated satellite gravity field missions, mainly the Gravity Recovery And Climate Experiment (GRACE, Tapley et al. 2004) and the Gravity field and steady-state Ocean Circulation Explorer (GOCE, Drinkwater et al. 2003), made possible the improvement of the quality of global gravity field models.The Earth Gravitational Model 2008 (EGM2008) was the first ultra-high GGM.It is complete to degree and order 2159 and contains additional coefficients up to degree 2190 and order 2159.This model was released before the availability of GOCE data and includes only four years of GRACE data combined with a global set of mean free-air gravity anomalies.EIGEN-6C4 was calculated with EGM2008derived gravity anomalies over continents combined with GOCE and LAser GEOdynamics Satellite (LAGEOS) data.GECO combines GOCE with EGM2008 data.SGG-UGM-2 is based on EGM2008-derived continental gravity data combined with GRACE and GOCE data.XGM2019e contains three main data sources: the combined satelliteonly model GOCO06s (Kvas et al. 2019), the 15′ ground gravity anomaly data set provided by the National Geospatial-Intelligence Agency (NGA) augmented with topographically derived gravity information over land (Earth2014, Hirt and Rexer 2015), and the 1′ augmentation data set consisting of gravity anomalies derived from altimetry over the oceans and topography over continents (Zingerle et al. 2020).All analysed models are published at the International Centre for Global Earth Models (ICGEM, Ince et al. 2019).
Argentina accounts for five proposed stations to be established as IHRF stations (Tocho et al. 2020).These stations are, from north to south: UNSA, OAFA, AGGO, UNPA and RIO2.AGGO (Argentinean-German Geodetic Observatory) is a fundamental geodetic station which co-locates all the main geodetic techniques together with a gravity lab.Although Tocho et al. (2020) compute preliminary potential values for the IHRF stations applying both methodologies described earlier, no evaluation of the accuracy of GEOIDE-Ar16 or the GGM-HR is presented.The results of this study allowed us to assess the capability of GGM-HR and GEOIDE-Ar16 to derive geopotential values in the IHRS.
2 Data sets and methods

GPS/levelling-based geoid undulations
The accuracy of a geoid model can be evaluated by an external comparison with geoid undulations derived from GPS and spirit levelling.Benchmarks with GPS-derived ellipsoidal heights and orthometric heights with respect to a local datum are significant types of data to include in databases because they enable the precise determination of geoid undulations on discrete sites.
Based on the known ellipsoidal and orthometric heights, geoid undulations ( / N GPS Levelling ) can be computed at GPS/ levelling benchmarks on the land according to the following equation: where h is the GPS-derived ellipsoidal height and H is the orthometric height.
Equation (1) has some limitations due to systematic and random errors in the derived heights h and H . First of all, / N GPS Levelling cannot be derived at the sea, which represents a limitation for the evaluation of geoid models along the oceans.On the other hand, there are systematic and gross errors in levelling, especially at higher altitudes.Levelling points are often difficult to access, and they are sometimes covered by vegetation or destroyed.Other limitations are datum inconsistencies (h and H refer to different reference surfaces), assumptions and theoretical approximations made in the orthometric correction, and the effects of not taking into account the deflection of the vertical, which can cause an error in the / N GPS Levelling geoid determination (Tocho 2006).
/ N GPS Levelling can be derived with very high relative and absolute accuracy, but one of the main disadvantages is its poor resolution.The errors that affect the accuracy of the ellipsoidal heights are originated mainly from three sources: satellite or orbit errors, signal propagation and receiver errors (Fotopoulos 2003).Some of these errors are atmospheric effects produced by the ionosphere and troposphere, the Earth's body tides, the ocean tides, the atmospheric loading, orbital errors of the GPS satellites, the antenna phase centre, set-up effects of the antenna, and multipath effects.
For this study, 2287 benchmarks across Argentina were taken into account (Figure 1).Geocentric Cartesian coordinates ( ) X Y Z , , were given in the terrestrial reference frame Posiciones Geodésicas Argentinas 2007 (POSGAR07, Cimbaro and Lauría 2009).POSGAR07 is based on ITRF2005 with epoch 2006.632 and it was adopted by the Argentinean Geographic Institute (IGN) as the official terrestrial reference frame for Argentina.( ) X Y Z , , were transformed to geodetic coordinates ( ) φ λ h , , using the GRS80 reference ellipsoid (Moritz 2000).At these points, Mader orthometric heights in the Argentinean National Vertical Reference System 2016 (SRVN16, Piñón et al. 2016) are available.As orthometric heights are defined in the mean tide concept, ellipsoidal heights were transformed from the tide-free concept (h TF ) to the mean tide concept (h MT ) following Petit and Luzum (2010), Sánchez and Sideris (2017): (2 where φ is the geodetic latitude.In this way, all quanti- ties given in equation (1) are expressed in the same tidal concept: mean tide.

Global gravity field models of highresolution-based geoid undulations
Geoid undulations have also been computed at the 2287 GPS/levelling benchmarks using several GGMs-HR up to degree 2190 ( = n 2190 max ).All GGMs-HR are available at the International Centre for Global Earth Models (ICGEM, Ince et al. 2019).Table 1 enumerates the GGMs-HR used in this study and its main characteristics.For all GGMs-HR, the Earth's gravitational constant ( -GM GGM HR ) and its equa- torial radius (r e ) remain the same with the exception of the equatorial radius of EIGEN-6C4, which is slightly higher.No significant impact was found in the derived geoid undulations due to the difference in the radius.A denotes altimetry, S satellite (e.g., GRACE, GOCE, and LAGEOS), G ground data (e.g., terrestrial, shipborne, and airborne measurements), and T topography.
Evaluation of GGMs-HR in Argentina  133 The GGM-HR-based geoid undulations (N GGM-HR ) have been determined through the formula (Rapp 1997): where ζ GGM HR is the height anomaly and g FA is the free- air gravity anomaly, both of them computed from the corresponding series expansions based on the spherical harmonic coefficients of each model and the GRS80 reference system parameters.H 0.1119 is the attraction of an infinite Bouguer plate of height H assuming a constant density of / 2.67 g cm 3 (Heiskanen and Moritz 1967), N 0 is the zero-degree term and γ ¯is the mean normal gravity between a point on the ellipsoid and the corresponding point on the telluroid.
The computation of the GGM-HR height anomaly ( ζ GGM HR ) was performed on the original tidal concept of the GGM-HR excluding the zero-degree term, following (Barthelmes 2013): where n and m are the degree and order of the spherical harmonic, respectively, and n max is the maximum degree, R is the mean radius of the Earth, P ¯nm are the fully nor- malized associated Legendre functions, and C nm T and S nm T are the coefficients of the disturbing potential.
The zero-degree term (N 0 ) in equation (3) takes into account the difference between the Earth's and reference ellipsoid's gravitational constant (GM) and also the difference between the reference potential W 0 adopted for the IHRS (Drewes et al. 2016) and the normal potential U 0 on the reference ellipsoid.N 0 has been computed from the following equation (Rapp 1997):

Wo
, , where GM GRS80 is the gravitational constant of the GRS80, r is the geocentric radial distance of the point on the geoid, γ is the normal gravity at the reference ellipsoid computed using Somigliana's closed formula and the GRS80 ellipsoidal parameters, W 0 IHRS is the potential value adopted by the IHRS, and U 0 GRS80 is the normal potential on the reference ellipsoid GRS80.The last term in equation ( 5) is about cm 76 , and the complete N 0 is approximately −17 cm.All geoid undulations computed with equation (3) were transformed from the original tidal concept of the GGM-HR (Table 1) to the mean tidal concept following (Ekman 1989): where = k 0.30190, being consistent with the Love num- bers proposed in Petit and Luzum (2010).Note that the first equation of equation (6) stands for the transformation from the tide-free concept to the mean tidal concept, while the second corresponds to the transformation from the zero tide concept to the mean tidal concept.

Gravimetric geoid model of Argentina:
GEOIDE-Ar16 for all gravity observations.Within the RCR scheme, the residual gravity anomalies were computed using Helmert's Second Method of Condensation for the treatment of the masses above the geoid.The residual gravity anomalies were gridded by the Kriging method.The resultant grid was applied in the Stokes' integral using the spherical multiband Fast Fourier Transform (FFT) approach and the deterministic kernel modification proposed by Wong and Gore (1969).Finally, the GGM contribution to the geoid and the indirect effect due to Helmert's Second Method of Condensation were restored.
The GBVP has been solved with gravity field observables in the tide-free concept.The evaluation of the obtained gravimetric geoid model was based on its comparison with geoid undulations derived from 1904 GPS-Levelling benchmarks.The statistics of the absolute differences shows a maximum, minimum, mean, and standard deviation of 2.34, −0.71, 0.68, and 0.24 m, respectively.
The final solution of the geoid model (GEOIDE-Ar16) is the result of the computed pure gravimetric geoid fitted to the Argentinean vertical datum through the determination of a corrective trend surface.This corrector surface was computed using the co-located GPS/levelling benchmarks.Figure 2 shows the obtained Argentinean geoid model.

Ellipsoidal, orthometric, and GGM-HR geoid heights statistics on benchmarks
The statistics of the individual height data sets that will be used in this study is given in Table 2.The statistics of the GGMs-HR geoid undulations refers to the values computed from equation (3) at the GPS/levelling benchmarks using the full spectral resolution of each model are also shown in Table 2.
Mean values given in Table 2 evidence that there exists a large discrepancy (>5 cm) between the reference surface of SRVN16 and the equipotential surface that is specified by the IHRS conventional value and realized by the various GGMs over Argentina.

Comparisons between geoid undulations computed from different GGM-HR and geoid undulations derived from GPS/ levelling benchmarks
The quality of the geopotential model can be assessed by comparing the geoid undulations obtained from GPS/ levelling ( / N GPS Levelling ) to those obtained from a global geopotential model ( -N GGM HR ) at i benchmarks: In practice, there are many effects that make equation (7) not equal to zero (e.g.Fotopoulos et al. 1999, Kotsakis andSideris 1999).Some of these effects are (i) random errors in the values of h, H , and N , (ii) datum inconsistencies as each related quantity (h, H , and N ) may refer to a different reference, (iii) the applied geoid may not coincide with the pure gravimetric geoid if it was fitted to the GPS/levelling benchmarks, (iv) systematic effects and distortions, such as long wavelength errors in N GGM-HR , poorly modelled GPS errors, and vertical datum distortions caused by overconstrained levelling networks; (v) geodynamics effects (post glacial rebound, land subsidence, monuments instabilities and mean sealevel rise), and (vi) theoretical approximations to obtain orthometric heights.
The statistics of the differences between the GPS/ levelling and the GGM-based geoid undulations are given   3. From the results in Table 3, it is clear that XGM2019e offers an improvement in the agreement among ellipsoidal, orthometric, and geoid heights throughout Argentina.The improvement achieved matches the GPS-levelling benchmarks within 22 cm, while EGM2008 does not perform better than 25 cm.We can also see that the mean values are about 4 cm, which is consistent with the difference between W 0 IHRF and the W 0 (local vertical datum estimated by Tocho and Vergos 2015) taken into account for the establishment of SRVN16.
The differences between different geoid undulations for GGGMs-HR used in this study and those based on GPS/ levelling points over Argentina are depicted in Figure 3.
Datum inconsistencies and systematic effects are the most important contributions to the discrepancies in equation (7).order to minimize these deviations, the following formula is usually applied (e.g.Kotsakis and Sideris 1999) where x is an n × 1 vector of unknown parameters, a is a n × 1 vector of unknown coefficients and v is the residual random noise term.The parametric a x i T part is supposed to describe all datum inconsistencies and other systematic errors in the data sets.The usual four-parameter model is often used (Kotsakis and Katsambalos 2010).
where φ and λ are the latitude and longitude of the GPS/ levelling benchmarks, respectively and the parameters a, b, c, and d are calculated by the least-squares method.
The adjusted parameters provide a corrector surface to be applied to the geoid undulations given by the GGM-HR in order to minimize their differences with the geoid undulations observed at the GPS/levelling benchmarks.
In the context of this study, it is important to mention that the obtained geoid model after applying a corrector surface cannot be used to estimate IHRS coordinates.As mentioned in the introduction, only a pure gravimetric geoid can be used for the determination of IHRS coordinates, and the geoid model after the application of the corrector surface is not.The geoid model with the four-parameter model applied provides only a deformed geoid, which minimizes the differences with the GPS/levelling benchmarks.
A corrector surface was estimated for each GGM-HR applied in this study.The statistics of the differences between the GPS/levelling and the GGM-based geoid undulations after the four-parameter model was applied are given in Table 4.As expected from the adjustment, all mean values shown in Table 4 are zero, as the corrector surface provides a shift of the reference surface.The standard deviations in Table 4 are lower in comparison to those in Table 3 in all cases.Nevertheless, XGM2019e still offers the best agreement between the GPS/levellingbased and the GGMs-HR-based geoid undulations.This value rises up to 19 cm; meanwhile, EGM2008 gave the worst results after the fitting of about 23 cm.
As the main objective of this study was to investigate the capability of GGMs-HR to derive IHRS coordinates,   geoid undulations from each GGMs-HR without the fourparameter model were used to derive potential values for all GPS/levelling benchmarks in Argentina following Tocho et al. (2020).In order to analyse the discrepancies between the potential values derived from different GGMs-HR, the standard deviation for each benchmark was computed, taking into account the five GGMs-HR considered.
Figure 4 depicts the standard deviation as a function of the ellipsoidal height (h) of each benchmark.From this analysis, the expected precisions proposed by Rummel et al. (2014) and Sánchez and Sideris (2017) are confirmed.For benchmarks located in flat areas up to about 1,000 m height, the standard deviation stays, in general, below 1 m 2 s −2 , while for higher altitudes, the standard deviation considerably increases given the topography variability and the poor distribution of data.
The mean standard deviation of all benchmarks rises up to about 0.5 m 2 s −2 , i.e. 5 cm, value that is well above the expected precision for the IHRS coordinates.This result confirms that it is not possible to use the GGMs-HR involved in this study to derive potential values in the IHRS.
3.2 Comparison between geoid undulations computed from N GEOIDE-Ar16 and geoid undulations derived from GPS/levelling benchmarks As we have mentioned in Section 2.3, GEOIDE-Ar16 is not a pure gravimetric geoid since it has been adjusted to the local levelling network using a corrector surface, and therefore, it cannot be used to estimate the IHRS coordinates.To do so, we must work with the pure gravimetric geoid prior to fitting to the local network, which we call N Grav .Moreover, both geoid models GEOIDE-Ar16 and N Grav include the zero-degree term applying equation (5), but using the reference potential adopted by the IERS Technical Notes Table 5 shows the statistics of GEOIDE-Ar16 and the pure gravimetric geoid undulations (N Grav and N Grav new ) interpolated at the GPS/levelling benchmarks.
The geoid undulations from N Grav and N Grav new were compared with those derived at the GPS/levelling benchmarks.Statistics from the differences are summarized in Table 6.6 show discrepancies of about 26 cm, being consistent with the difference between the reference potential values given by IERS and the IHRS.The standard deviation presents the same variability as only the reference is changing between N Grav and N Grav new .From the results of the standard deviation shown in Table 6, GEOIDE-Ar16 cannot be used to recover IHRS coordinates, and therefore, a new gravimetric geoid model must be determined consistent with the definition and recommendations of the IHRS.

Conclusions
A qualitative analysis of different GGMs-HR is presented in order to derive geoid undulations in Argentina and compare them with GPS/levelling benchmarks.In this article, five GGMs-HR up to degree 2190 are evaluated: EGM2008, EIGEN-6C4, GECO, SGG-UGM-2, and XGM2019e.2287 GPS/levelling benchmarks are taken into account for the analysis where ellipsoidal and orthometric heights are available.Orthometric heights are given in the local vertical datum of Argentina: Argentinean National Vertical Reference System 2016 (SRVN16).Together with the GGMs-HR, the current official geoid model of Argentina (GEOIDE-Ar16) is also evaluated.
In order to estimate the geoid undulations from the GGMs-HR and GEOIDE-Ar16, special care is taken into account for the treatment of the permanent tides, the use of the reference ellipsoid, and the zero-degree term, including the difference between the reference surface defined for the IHRS and the ellipsoid (GRS80).To ensure the consistency of the analysis, all quantities are obtained using the mean tide concept.
From the absolute differences between the geoid undulations based on the GGMs-HR and those derived at the GPS/levelling benchmarks, XGM2019e shows the best agreement, with a standard deviation of 23 cm.After fitting a corrector surface for each GGM, this model still provides the best fit with a standard deviation of 19 cm.
In order to assess the capability of the analysed GGMs-HR to infer IHRS coordinates in Argentina, potential values are obtained for all benchmarks.The mean standard deviation of the obtained potential values from the GGMs-HR rises to 0.5 m 2 s −2 , meaning about 5 cm in terms of height.This result allows concluding that the GGMs-HR evaluated in this study are not precise enough to derive IHRS coordinates.GEOIDE-AR16 is also compared to the GPS/levelling benchmarks.Since it is a pure gravimetric geoid fitted to the vertical datum of Argentina by using a corrector surface, all analyses were performed with the geoid model prior to the fitting.The pure gravimetric geoid model of Argentina fits slightly better than XGM2019e, showing a standard deviation of 21 cm.This result also suggests that it cannot be used for the estimation of IHRS coordinates and emphasizes the need for a new more precise computation of the gravimetric geoid in Argentina.

Figure 1 :
Figure 1: Distribution of the GPS/levelling benchmarks in Argentina.The topography is also included, noting its variability along the country.

Figure 4 :
Figure 4: Standard deviation of potential values derived from each GGMs-HR as a function of the ellipsoidal height of each benchmark (green dots).The red-dashed curve indicates the mean standard deviation from all benchmarks.
potential of IERS and the one established for the IHRS, a new geoid undulation N Grav new was obtained.

Table 1 :
GGMs-HR used for the evaluation tests at the 2287 GPS/levelling benchmarks (Piñón 2016)is the official geoid model for Argentina(Piñón 2016).It is obtained by solving the scalar-free GBVP using Stokes' theory.GEOIDE-Ar16 was developed by the IGN, together with the Royal Melbourne Institute (Jarvis et al. 200813-University of Australia, using the RCR technique(Schwarz et al. 1990), on a grid with a spatial resolution of 1′ × 1′.The GOCO05S Global Gravity Model (Mayer-Guerr 2015) completed up to degree and order 280, together with 671,547 gravity measurements observed in the International Gravity Standarization Net 1971 (IGSN71,Morelli et al. 1972) on the Argentine continental territory, its neighbouring countries, Islas Malvinas and the coastal (marine) area were used for the computation of the gravimetric geoid.Regions with a lack of gravity observations were filled with data from the DTU13(Andersen et al. 2013) world gravity model.Terrain corrections were calculated using a combination of the SRTM v4.1(Jarvis et al. 2008) and SRTM30 Plus v10(Becker et al. 2009) Digital Elevation Models (DEMs)

Table 2 :
Statistics of ellipsoidal, orthometric, and geoid heights on benchmarks in meters

Table 4 :
Statistics of the differences between the GPS/levellingderived geoid undulations and those computed from the different GGMs-HR after the four-parameter transformation model has been applied in meters

Table 3 :
Statistics of the differences between geoids undulations derived from different GGMs-HR and GPS/levelling-derived geoid heights in meters

Table 5 :
Statistics of GEOIDE-Ar16 and the pure gravimetric geoid models, N Grav and N Grav new , in meters

Table 6 :
Statistics of the differences between undulations N Grav new with GPS/levelling-derived geoid in meters / N N − Figure 5: Differences between geoid undulations derived from GPS/ levelling and N Grav at different benchmarks.