Global Geopotential Models assessment in Ecuador based on geoid heights and geopotential values

1 Introduction

The Global Geopotential Models (GGMs), according to Pavlis (2013), are mathematical approximations of the external gravitational potential of an attractive body (the Earth in the case of geodesy). A GGM consists of a set of numerical values for certain parameters, statistics of errors associated with these values, and a collection of mathematical expressions, numerical values, and algorithms that allow calculating magnitudes related to gravitational potential (synthesis) and associated errors (error propagation based on a covariance matrix). The development of a high-resolution GGM involves the appropriate combination of varied information from orbital platforms, ground surveys, ocean surveys, airborne sensors, and altimetry. Combined models associate various sources of gravity information (satellite orbits, satellite gravimetry missions, terrestrial, aerial and oceanic gravimetry, satellite altimetry, and Digital Elevation Models) with reduced omission errors (or truncation), but with great dependence on commission errors due to the diversity of references involved.

GGMs have both operational and scientific uses. Among the most important, we can mention the following: determination of satellite orbits; determination of the trajectory of missiles and aircraft (inertial navigation); calculation of geoid heights; the estimate of the Dynamic Ocean Topography, or even of the Mean Ocean Dynamic Topography; study of the ocean circulation for marine applications; realization of the Global Vertical Datum (GVD), considering as a reference the geopotential surface materialized by a GGM; and geophysical prospecting (determination of mass distribution) (Pavlis, 2013).

Among the applications of the GGMs and emphasizing the establishment of the geodetic reference frame, we highlight in this manuscript the realization of the GVD as an equipotential surface of the Earth’s gravity field with geopotential W ₀. The GVD represents one of the key elements for the establishment of a Global Height System, as indicated by the International Association of Geodesy (IAG) in its resolution for the definition and realization of the International Height Reference System (IHRS) (Drewes et al., 2016).

According to the conventions established by the IAG (Drewes et al., 2016) and disseminated at the regional level by Geodetic Reference System for the Americas (SIRGAS), in relation to the establishment of a Global Vertical Reference System, the heights belonging to the fundamental leveling networks must be obtained according to the calculation of geopotential numbers (C) and be linked to a GVD (materialized by an equipotential surface of the Earth’s gravity field with a geopotential value W ₀ = 62636853.4 m² s⁻²) (Ihde et al., 2017).

Furthermore, as stated by Ihde et al. (2017), the IHRS realization (International Height Reference Frame [IHRF]) requires the combination of a global station network, a GGM, and a set of parameters. Thus, the long and medium wavelengths of the gravity field, recovered from a GGM, are one of the key elements for determining geopotential values for the fundamental stations constituting the IHRF.

The evolution towards a GVD has originated in several countries with the beginning of activities aimed at the realization of reference stations linked to the GVD (IHRF stations), considering mainly the calculation of the discrepancy δ W = W 0 − W 0 i . In the SIRGAS region, the countries' members have started activities aimed at the realization of a Global Vertical Reference System for the determination of physical heights and referred to GVD. These activities are currently focused on the determination of the optimal positions (according to Ihde et al. (2017), “the stations must be at least continuously monitored to detect deformations, referred to the International Terrestrial Reference System/International Terrestrial Reference Frame (ITRF), and be connected by levelling with the local vertical datum”) and the geopotential calculation for the IHRF stations, which serve as fundamental geodetic marks for the IHRS realization (Blitzkow et al., 2017; Carrión et al., 2023; Carrión Sánchez et al., 2018; Guimarães et al., 2021, 2022; SIRGAS WG-III, 2021; Tocho and Vergos, 2016; Tocho et al., 2022, 2023).

As mentioned by Sánchez et al. (2021), a pragmatic way for the calculation of geopotential numbers referred to an equipotential surface of the gravity field is the direct calculation of W(P), which is accomplished by introducing the ITRF coordinates (Altamimi et al., 2016) of a point into the harmonic expansion corresponding to a high-resolution GGM. In this sense, it should be highlighted the technological advances in satellite gravimetry achieved with the Gravity Recovery and Climate Experiment (GRACE) (Tapley et al., 2004) and Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) missions (Drinkwater et al., 2003), reaching accuracies of ±0.2 m² s⁻² (equivalent to a vertical positional accuracy of ±2 cm) for medium and long wavelengths of the gravity signal, i.e., to a maximum spatial resolution of 100 km (Rummel et al., 2002).

Furthermore, geopotential values W _P can also be calculated by solving the geodetic boundary value problem (GBVP) for regional pure gravimetric geoid/quasigeoid determination, considering the remove-compute-restore scheme. In this case, the information derived from the GGMs is also fundamental. Thus, one of the main inputs for the calculation of geopotential values in IHRF stations is the knowledge of the long (up to degree 70), medium (up to degree 200), and short (about degree 250–200) wavelengths of the gravity field (Sánchez et al., 2021), which can be obtained from the high-resolution GGM. In this sense, it is evident that there is a need to carry out studies that allow its validation in local and regional scenarios. In the present work, geoid heights from GGMs are compared with geoid heights obtained for the level references of the Ecuadorian vertical control network.

Several studies have been carried out with the purpose of validating the GGMs. These works are based on comparing quantities such as geoid heights, gravity anomalies, and gravity disturbances with those derived from global models. The Ecuadorian Vertical Reference System (EVRS) must adapt to the precepts established by the IAG and disseminated in South America by SIRGAS. The establishment of the vertical datum for the EVRS must be based on the calculation of the geopotential ( W 0 i ) at one or more fundamental stations (i.e., IHRF stations).

2 Methods

The GGMs are implemented based on a spherical harmonic expansion; the low and medium frequencies of the series are derived from the analysis of the satellite orbits, the satellite-to-satellite tracking, and gravimetric gradiometry. The high degree and order of expansion are achieved by combining the information from satellites with terrestrial-oceanic gravimetry, Digital Elevation Models, and products derived from satellite altimetry (Torge and Müller, 2012).

Geodesists have used several methods to represent the Earth’s gravity field, and among these methods are mass points (Sünkel, 1981), finite elements (Meissl, 1981), and splines (Sünkel, 1984). Earth’s gravity field modelling using these approaches is limited to some specific applications, while the representation with spherical harmonics has prevailed as the standard method for representing the gravity field globally. The external gravitational potential (V) of a point P defined by its geocentric distance (r _P ), geocentric colatitude (θ _P ), and longitude (λ _P ) can be calculated according to the following expression (Heiskanen and Moritz, 1985):

It is observed that in equation (1), the coefficients of degree and order one (C ₁₀, C ₁₁, and S ₁₁) were assumed to be equal to zero (it refers to the first-order moments of inertia that determine the coordinates of the center of mass), which implies assuming that the system originates from the center of mass. Being GM the geocentric gravitational constant. GM is a constant, and a related to the spheroidal model, that acts as a scale factor associated with the fully normalized spherical harmonic coefficients (C _nm ). The surface harmonic functions S _nm are defined as follows:

where P ̅ n | m | ( cos θ P ) is the fully normalized Legendre function with degree n and order m.

Gravity field modelling allows the estimation (in the most accurate way) of the C _nm coefficients through the optimal combination of gravity information from various data sources. The C _nm coefficients are used for the calculation of various functionals of the Earth’s gravity field, and therefore, the estimation of the errors associated with these coefficients is an important aspect of model generation (Pavlis, 2013).

The zero-degree term, on the geopotential quantities calculation, refers to the difference (W ₀−U ₀) between the geopotential W ₀ adopted as vertical datum (e.g., which was selected as datum for the IHRF; Drewes et al., 2016) and the normal potential U ₀ corresponding to a specific reference ellipsoid. Furthermore, the zero-degree term also takes into account the difference between the GGM’s geocentric gravitational constant (GM_GGM) and the corresponding geocentric gravitational constant associated with the reference ellipsoid (GM_GRS) (e.g., GRS80) (Lemoine et al., 1997).

Thus, the zero-degree term, referring to the calculation of the anomalous potential (T ₀), is given as follows (Hofmann-Wellenhof and Moritz, 2006):

In the present study, the validation of the GGMs shown in Table 1 is accomplished in terms of geoid heights, using Global Navigation Satellite System (GNSS)/lev records collected within Ecuadorian territory (Figure 1) and also considering geopotential values for the Ecuadorian Vertical Datum (EVD).

Figure 1

Study area.

In addition to the analysis carried out with the GNSS/lev records, geopotential values derived from GGM’s height anomalies were evaluated through comparison with the geopotential value obtained for the EVD as a solution to the fixed GBVP (Carrión et al., 2023).

The analyses were carried out with the following combined GGMs: EIGEN6C4 (Förste et al., 2014), EGM2008 (Pavlis et al., 2012), GECO (Gilardoni et al., 2016), SGG-UGM1 (Liang and Reißland, 2018), SGG-UGM2 (Liang et al., 2020), and XGM2019e (Zingerle et al., 2020).

2.1 GNSS/lev records

Geometric leveling does not consider the effect of the non-parallelism among the equipotential surfaces of the gravity field. The determination of orthometric heights (H) requires the availability of gravity records along the spirit levelling lines. Since not all leveling records from Ecuador’s vertical control network have associated gravity observations (gravity observations have been carried out on 3,236 level references from 4,192 – Figure 2), gravity values for the remaining records were obtained through interpolation from existing gravity data.

Figure 2

Observed and predicted gravity data for the Ecuadorian vertical control network.

To predict gravity values, the least-squares collocation method, which is implemented in the PREDGRAV computational package, was used (Drewes, 1976, 1978). Thus, gravity within a radius of 100 km or the nearest 50 records from the interpolated points was considered. PREDGRAV performs interpolation for simple Bouguer anomalies (Δg _B ) (4) and then, following the inverse procedure, calculates the interpolated gravity values.

(4) ∆ g B = g + C FA + C B − γ ,

where g is the observed gravity, C _FA is the free-air correction, C _B is the Bouguer correction, and γ is the normal gravity.

To evaluate the predicted values, we followed a cross-validation procedure, considering 323 gravity records that were randomly selected to cover 10% of the entire gravity dataset. The corresponding statistical values are presented in Table 2.

Table 2

Statistical for predicted gravity from the cross-validation procedure (values in mGal)

Min.	−28.30
Max.	13.34
Mean	−0.10
Std.	2.69

Thus, the orthometric correction was applied according to the expression given by Heiskanen and Moritz (1967):

(5) OC AB = g − γ 0 γ 0 A B ∆ n + g A − γ 0 γ 0 H A − g B − γ 0 γ 0 H B ,

where OC _AB is the orthometric correction for the segment defined between points A and B, Δn is the levelling increment, g is the average gravity along the levelling line, γ 0 is the normal gravity for a standard latitude, and g _A and g _B are the mean gravity values along the plumb line between the physical surface and the geoid.

For the calculation of g _A and g _B , the simplified prey reduction was used considering a normal density of ρ = 2.67 g cm⁻³.

(6) g = g + 0.0424 H .

H _A and H _B are the corresponding orthometric heights for the leveling segment and are calculated iteratively as follows:

(7) H ort = H lev + OC .

Then, orthometric (H _ort) and ellipsoidal heights (h) for Ecuadorian vertical control network records are used to infer geoid heights such as this expression:

(8) N ≈ h − H ort .

Geoid heights obtained from equation (8) are used to validate geoid heights (N) obtained from GGMs. For this study, 3,104 GNSS/lev records from the Ecuadorian vertical control network were used (Figure 3). All the calculations were performed considering the mean-tide concept, thus the ellipsoidal heights from Global Navigation Satellite System (GNSS) positioning are transformed from the tide-free to the mean-tide concept, according to the expression given by Rapp (1989):

(9) H mean − tide = H tide − free − ( 1 + k − h ) ∆ W zero − tide g ,

Figure 3

GNNS/lev data.

where h and k are the love numbers (h = 0.605; k = 0.30190), and the term ∆ W zero − tide g is obtained with sub-centimeter accuracy as follows (Heikkinen, 1978):

(10) ∆ W zero − tide g ≈ − 0.198 3 2 sin 2 φ ̅ − 1 2 [ m ] ,

where φ ̅ is the geocentric latitude of the calculation point.

2.3 Geoid height discrepancies and outlier detection

The calculation of discrepancies (ΔN, equation (14)) between geoid heights for the level references (N) and those from the GGMs (N _GGM) was carried out initially for an outlier detection’s stage, considering only the geoid heights coming from the GGM XGM2019e at its maximum degree (Förste et al., 2014).

(14) ∆ N = N GNSS / lev − N GGM .

The criterion for the detection of outliers was considered as expected values that were all within the interval [−2σ +2σ]. The 2-σ criterion was established empirically, considering a normal distribution for which the probability that a value (N _GNSS/lev − N _GGM) belongs to the set of expected values is at least 95.45%.

The statistics for ΔN, before and after outlier exclusion, are shown in Table 3.

Table 3

Outlier elimination N _r statistics

	Before outlier exclusion	After outlier exclusion
Min. (m)	−3.620	−0.967
Max. (m)	3.531	1.071
Mean (m)	0.047	0.018
Std. (m)	0.509	0.379
N°	3104	3000
%	100	96.65

High maximum/minimum and standard deviation values for ΔN statistics (before outlier exclusion) shown in Table 3 are due to gross errors in heights (leveled/ellipsoidal), registered in the observation, processing, or handling of data collection’s stages.

Furthermore, Table 3 shows the extreme values close to 1 m and the standard deviation of 0.3791 m calculated after eliminating 104 records detected as outliers. These values are graphically represented in the frequency distribution of Figure 4, in which, additionally, it can be seen that approximately 61% of the calculated ΔN are in the range [−0.4 to 0.2] m.

Figure 4

Frequency distribution for N _GNSS/lev − N _GGM (m).

Figure 5 shows the spatial distribution of the resulting discrepancies (N _GNSS/lev − N _XGM2019e). An increase in N _r values is observed in the Andes mountain range region, whose highest heights in Ecuador are located approximately between longitudes 79°W and 77°W (Figure 5). This is explained by the fact that in mountainous regions, the omission error of the GGM tends to increase due to deficiencies in the representation of the topographic masses effect for gravity field modeling (high frequency of the gravity field).

Figure 5

N _GNSS/lev − N _GGM.

The calculated ΔN are interpreted through the degree of adherence of the two series observed in Figure 6, which shows the values of N and N _XGM2019e for all points considered in the analysis. Discrepancies between geoid heights from GNSS/lev records and those from the GGM XGM2019e are mainly due to errors that affect the accuracy of leveled and ellipsoidal heights, omission/commission errors of GGMs, inconsistencies in the local datum, and systematic effects.

Figure 6

N _GNSS/lev and geoid heights from XGM2019e.

Model omission errors tend to increase in mountainous regions due to deficiencies in the representation of the effect of topographic masses on the gravity field. Figure 7 shows the variation of geoid heights and ΔN (N – N_XGM2019e) as a function of orthometric height (a) and longitude (b). The influence of the mountainous region’s topographic masses (the Andes mountain range) can be interpreted in Figure 7(a), in which a slightly increasing tendency of geoid heights associated with the increase in orthometric heights can be observed. In Figure 7(b), it is observed that the computed and modeled N values and also the corresponding ΔN values, increasing in the rank of longitudes corresponding to the Andean region.

Figure 7

(a) N _GNSS/lev, N _XGM2019e, and N _res vs H _ort. (b) N _GNSS/lev, N _XGM2019e, and N _res vs longitude.

3 GGMs assessment based on geoid heights

Once outliers had been eliminated from the GNSS/lev dataset, the analysis oriented to the assessment of several combined GGMs was carried out. In this study, the combined GGMs shown in Table 1, considering their maximum degree, were evaluated. The analysis was performed in terms of geoid height discrepancies (ΔN) in equation (14), following for its calculation the same approach that was used for the outlier detection stage. The statistics for ΔN are shown in Table 4.

The standard deviation, in all cases, exceeds 40 cm, with the smallest value corresponding to the GGM XGM2019e. The expected accuracy of the GGM for the study region is in the range of [40, 47] cm, with mean values close to [20, 24] cm. As can be seen in the frequency distribution in Figure 8, for all cases, approximately 70% of ΔN are in the range [−50 50] cm, while the interquartile range, as shown in Figure 9, is approximately [−20 30] cm.

Figure 8

Frequency distribution for N _GNSS/lev − N _GGM.

Figure 9

Boxplot for N _GNSS/lev − N _GGM.

As an additional element, regarding the statistical analysis, correlations (Corr.) close to 1 between the calculated geoid heights and those derived from the GGMs (Table 4) are shown graphically in Figure 10.

Figure 10

N _GNSS/lev and N _GGM correlation.

Moreover, based on the approach proposed by Kotsakis and Sideris (1999) and following the analysis carried out by Tocho et al. (2022), the discrepancies between the calculated and modeled geoid heights are minimized by adjusting the parameters of the following formula:

where v i are the residuals, and a i T x can be calculated as follows (Kotsakis and Katsambalos, 2010):

where φ i and λ i are the latitudes and longitudes for the GNSS/lev records, and a, b, c, and d are the parameters calculated by the least-squares method.

Thus, by the least-squares adjustment, a correction surface, established by the computed parameters (four-parameter model), was estimated for each GGM, minimizing the residuals ( v i ) of equation (15) and theoretically minimizing local vertical datum inconsistencies and systematic errors of the datasets (Tocho et al., 2022). Analogous to the statistics presented in Table 4, new statistics were calculated this time considering the geoid heights affected by the correction surface.

The statistics for ΔN (Table 5), after applying the correction surface, show for all the cases a reduction in the standard deviations (Std.), with a value of 0.355 m for XGM2019e, which best fits the geoid heights from the GNSS/lev records. Furthermore, for all cases, the null mean values denote the influence of an offset on the reference surface.

Figure 11 shows the standard deviations for the adjusted and unadjusted geoid height discrepancies (ΔN). For all cases, there is a decrease in the Std. of about 5 cm. Note also that, according to the Std., the worst agreement between GNSS/lev and GGMs geoid heights corresponds to EGM2008, for which Std. of 0.469 and 0.440 m were calculated for the unadjusted and adjusted geoid heights, respectively.

Figure 11

Std. [N _GNSS/lev − N _GGM].

According to Tocho et al. (2022), the previously described procedure generates a deformed surface (geoid) that better fits the GNSS/lev records, and only a pure gravimetric geoid can be used for the determination of IHRS coordinates. Therefore, this surface cannot be used to determine geopotential values in the context of the IHRS realization. However, the described procedure allows for estimating the effects related to datum inconsistency and other systematic effects contained in the data sets.

4 GGMs assessment based on geopotential values

As a second stage in the evaluation of the GGMs, the analysis was carried out in terms of geopotential values, which were compared with the corresponding value calculated for the EVD (La Libertad tide gauge) through the fixed solution of the GBVP (Carrión et al., 2023).

According to Barthelmes (2013), the geopotential is calculated as follows:

where r, λ, and φ are the spherical geocentric coordinates of the computation point; and C l m W and S l m W are the fully normalized spherical harmonic coefficients.

The geopotential value for the EVD (equation (17)) was obtained from the ICGEM calculation service, considering the maximum degree for the evaluated GGMs, the zero-tide as permanent tide concept, and finally the zero-degree term was included.

Because W ₀ is determined using the conventions defined by the IAG, the zero-degree term must be considered in the context of the unification of vertical reference systems and the establishment of the IHRF. However, the calculation of quantities originating from GGM is independent of this term and should only be considered to rescale these quantities in relation to a given geoid (equipotential surface).

As recommended by Sánchez et al. (2021), the calculations were performed in the zero-tide concept, and later, the geopotential values were transformed into the mean-tide concept. According to Mäkinen (2017), this transformation is performed by including the latitude depending on the term W _T .

The differences between the calculated and the modelled values of the GGMs vary between approximately 0.01 and 1.25 m² s⁻², with the XGM2019e (as in the case of the analysis for geoid heights) being the model that best approximates the geopotential value calculated through the fixed GBVP solution. Eigen6c4, on the other hand, presents a considerably higher difference in relation to the other models. The corresponding differences in terms of geopotential values are shown in Table 6.

Figure 12 shows the differences found in this analysis graphically.

Figure 12

Calculated/modeled geopotential differences.

5 Conclusions

As a result of the GGM validation based on the comparison of geoid heights with GNSS/lev records in Ecuador, it was found that the expected accuracy of the GGM for the study region is in the range of [40, 47] cm and the mean values are in the range of [20, 24] cm.

The comparison of GGM geopotential values with the corresponding value obtained through the GBVP solution denotes that the differences between the calculated and the modelled values of the GGMs vary approximately between 0.01 and 1.25 m² s⁻².

For the two analyses carried out (with geoid heights and geopotential values), the model with the best performance was the XGM2019e.

Authors as Rummel (2012), Sánchez and Sideris (2017), have concluded that the expected mean accuracy for the GGMs in well-surveyed areas is ±0.4 to ±0.6 m² s⁻², and ±2 to ±4 m² s⁻² for sparsely surveyed areas. In this sense, in relation to the potential differences in Table 6 [∽0.01 ∽1.25] m² s⁻², we can mention that, although the values obtained are close to the range that corresponds to well-surveyed areas, it is not possible to carry out an analysis in terms of precision because these values are estimated for a single point.

According to Ihde et al. (2017), although the target accuracy for the IHRS geopotential numbers is 0.01 m² s⁻², in practice, this value is in a range of 0.1–1 m² s⁻². This is due to factors such as non-unified standards for: determination of potential values, gravitational field modeling, estimation of position vectors, heterogeneity in the distribution of references that make up the geodetic infrastructure, etc. Thus, the estimated accuracies for the GGM in Ecuador, which are at least 2 orders lower than the target value, allow us to conclude that the determination of geopotential values, exclusively from GGMs, for IHRF fundamental stations in Ecuador is not feasible.

The zero-degree term considered in the calculation of geoid heights has an approximate value of −17 cm. Its influence on the calculation is not negligible, and it was one of the main causes of discrepancies in the results obtained for the Colorado experiment, carried out for the standardization of geopotential calculation procedures in the context of the establishment of the IHRS (Wang et al., 2021).

To be consistent, and following the recommendations given by the IAG on the IHRS realization, the calculations were performed in the zero-tide concept, and at the end of the calculation procedure, the geopotential values were transformed to the mean-tide concept.

Orthometric heights were calculated with the aim of giving a physical meaning to heights used in the analyses carried out in this study.

Through an adjustment based on the least-squares criterion, a surface that best fits the geoid heights of the GNSS records/lev was established. However, the adjustment reduces the discrepancies between the GNSS/lev geoid heights and the GGM’s values, and this procedure generates a deformed surface that cannot be used to determine geopotential values in the context of the IHRS realization.