Combining Global Geopotential Models, Digital Elevation Models, and GNSS/Leveling for Precise Local Geoid Determination in Some Mexico Urban Areas: Case Study

Abstract

This work shows improvements of geoid undulation values obtained from a high-resolution Global Geopotential Model (GGM), applied to local urban areas. The methodology employed made use of a Residual Terrain Model (RTM) to account for the topographic masses effect on the geoid. This effect was computed applying the spherical tesseroids approach for mass discretization. The required numerical integration was performed by 2-D integration with 1DFFT technique that combines DFT along parallels with direct numerical integration along meridians. In order to eliminate the GGM commission error, independent geoid undulations values obtained from a set of GNSS/leveling stations are employed. A corrector surface from the associated geoid undulation differences at the stations was generated through a polynomial regression model. The corrector surface, in addition to the GGM commission error, also absorbs the GNSS/leveling errors as well as datum inconsistencies and systematic errors of the data. The procedure was applied to five Mexican urban areas that have a geodetic network of GNSS/leveling points, which range from 166 to 811. Two GGM were evaluated: EGM2008 and XGM2019e_2159. EGM2008 was the model that showed relatively better agreement with the GNSS/leveling stations having differences with RMSE values in the range of 8–60 cm and standard deviations of 5–8 cm in four of the networks and 17 cm in one of them. The computed topographic masses contribution to the geoid were relatively small, having standard deviations on the range 1–24 mm. With respect to corrector surface estimations, they turned out to be fairly smooth yielding similar residuals values for two geoid models. This was also the case for the most recent Mexican gravity geoid GGM10. For the three geoid models, the second order polynomial regression model performed slightly better than the first order with differences up to 1 cm. These two models produced geoid correction residuals with a standard deviation in one test area of 14 cm while for the others it was of about 4–7 cm. However, the kriging method that was applied for comparison purposes produced slightly smaller values: 8 cm for one area and 4–6 cm for the others.

1. Introduction

Global Geopotential Models (GGM) have become an essential tool in geodesy as well as in other Earth Sciences and engineering field applications. For instance, for surveying, mapping, and engineering projects GGM can be combined with precise Global Navigation Satellite System (GNSS) positioning to efficiently yield orthometric heights to benchmarks over the Earth’s surface, avoiding the costly and time demanding spirit leveling. Several GGMs have been produced during the last decades. At present there are more than 170 GGMs available at the International Center for Global Gravity Field Models (ICGEM) (Potsdam Germany), in the form of fully normalized spherical harmonic coefficients that can be used to compute geodetic and Earth’s gravity field quantities (http://icgem.gfz-potsdam.de). The different GGMs vary with respect to accuracy and resolution. The accuracy mainly depends on the type quality and amount of data used for the estimation of the model’s parameters, and the error associated is called commission error. With respect to resolution, which is the wavelength of the frequency related to the highest degree and order of the spherical harmonic expansion of the model, and since the mathematical representations of these models are truncated series expansions, there is an omission error involved [1,2,3,4]. Some of the most relevant GGM models with respect to both accuracy and resolution are the EGM2008 [5], GECO [6], and XGM2019e_2159 [7] models; they have the highest resolution among all available modes being of about 2190 with respect to degree and order. The related spatial resolution is around 9–10 km depending on latitude. Some previous studies have found that the omission error of high resolution GGMs could reach up to 15 cm or more in mountainous areas and the mean value is of about 5 cm [8]. This error is mainly due to the topographic masses corresponding to the wavelengths of higher frequency content that is not present in the GGM. It can be reduced by employing a Residual Terrain Model (RTM) through a Gravity Forward Modeling (GFM) procedure [3,4,9,10,11,12,13,14,15,16,17,18,19]. The RTM can be obtained by subtracting, from a Digital Terrain Model (DTM), the frequency content common to a certain GGM. Thus, the RTM will be composed of all the higher frequency content that the GGM is lacking, being part of the omission error [2,3,4,13,16,20,21,22,23,24,25,26,27]. The removal of the common masses from the Digital Terrain Model (DTM) can be performed either by subtracting a reference field also called topographic reference model or Topographic Gravity Field Model (TGFM), which contains the common frequency content or, by applying a low pass filter, keeps the RTM higher frequencies.

With respect to the DTM, some of the most commonly employed models are produced by the Shuttle Radar Topography Mission (SRTM), which is a joint project between the U. S. National Imagery and Mapping Agency NIMA and the National Aeronautic and Aerospace Administration NASA (https://www2.jpl.nasa.gov/srtm). A 3-arc-second SRTM DEM (STRM3) for many parts of the world has been compiled and released [24], SRTM V4.1, with SRTM 1 arc-second (30 m) global elevation data of the world-wide coverage. As mentioned above, the RTM technique requires the use of a reference field topography through a TGFM that contains the frequencies of the topographic masses included in the GGM. A TGFM is produced by the gravitational potential generated by the attraction of the Earth’s topographic masses. These masses are composed (constituent) mainly by rock, sand, water, and ice. TGFM models such as DTM2006.0 [22] and EARTH2014 [28] are given in terms of spherical harmonic coefficients and can be evaluated at any location on the Earth’s surface [2,4,18,21,22,24,29,30,31]. It is recommended that this reference topography model be the same as the one used in the generation of the GGM model. For instance, in the establishment of EGM2008 model, the DTM2006.0 topography model was employed [22].

The effect on the geoid produced by the masses of an RTM can be computed by the so called Forward-Modelling (FM) method for gravitational potential. This method denotes the effect of the associated topographic masses assuming certain distribution and density of them. The FM techniques are based on Newton’s law of universal gravitation. Thus, a volume integration has to be applied. This one can be analytically executed similar to prism methods ([11,17,32,33,34,35]) or numerically such as with the tesseroid methods [15,36,37,38,39]. The prism scope is the most precise and rigorous but at the same time is more time consuming due to the several logarithmic and arc tan functions that have to be evaluated [17,32]. A significant reduction in computation time can be achieved by applying fast Fourier transform techniques [12,40,41].

In the treatment of the topographic masses either for TGFM establishment or for FM application, assumptions have to be made about their densities, considering their different types like rock, water, and ice. For this mater, two main approaches have been proposed: One is the rock water and ice Rock–Water–Ice (RWI) approach [2,13,29]. With this technique, the Earth’s topography is decomposed in three layers each one with its associated mass density. The other approach uses the concept of rock-equivalent topography, rock-equivalent topography (RET), or rock-equivalent heights (REQ) [13,16,29,30,42]. With this technique, the areas of water and ice are condensed so that they form a layer with thickness equivalent to the rock density, producing a change in mass distribution, but having a unique density. That is, the topography is replaced by a layer of constant density and with the same mass as the original layer. The mass is equivalent but with rock density. The differences between the RWI and the REQ-based model reach maximum amplitudes of about 1 m globally [29].

Once geoid undulations from a GGM are corrected by the topographic masses effect, a validation scheme can be performed by using a set of GNSS/leveling stations distributed in a certain area or region. From them, we can get independent and precise geoid undulation values, which can also be used to improve the GGM geoid by estimating biases, trends and other systematic errors resolved through a corrector surface and a suitable least squares adjustment model [2,4,23,27,43,44,45,46,47]. This process can also be applied to a gravimetric geoid. Toward that end, we also employ the Mexican gravimetric geoid GGM10 [48].

The main purpose of the present work with local scope is to combine a GGM with a DTM and a local network of GNSS/leveling stations to improve geoid undulations as produced by the GGM. The analyzed models are EGM2008 and XGM2019e_2159. The numerical analysis is made on five Mexican urban areas.

2. Numerical Technique

2.1. Gravity Forward Modeling

Using the volumetric integral of mass elements to compute the gravitational potential $V_P$ at a point P with spherical geocentric coordinates, i.e., radius $r_p$, colatitude ${\theta }_p$, and longitude ${\lambda }_p$, from [49]:

$$V_P=G\iiint \frac{{\rho }_Q}{{\ell }_{PQ}}{r}_Q{}^{2}d{\sigma }_Q$$

(1)

where $Q$ is the integration point with corresponding spherical coordinates, $r_Q, {\theta }_Q, {\lambda }_Q$, which is also the position of the atracting mass volumetric element $d{\sigma }_Q=\sin {\theta }_{Q}drd\theta d\lambda$. The integration is run over the attracting mass or body volume, $G$ represents the universal gravitational constant, ${\rho }_Q$ is the mass density at integration point, and ${\ell }_{PQ}$ is the Euclidean distance between evaluation and attraction point.

The gravitational potential can also be formulated for a spherical mass layer, ibid. p5.

$$V_P=G\iint \frac{dm}{{\ell }_{PQ}}=G\iint \frac{{\tau }_Q}{{\ell }_{PQ}}d{s}_Q$$

(2)

where

$${\tau }_Q=\frac{dm}{d{s}_Q}={\rho }_Q{H}_Q$$

(3)

then

$$V_P=G\iint \frac{{\rho }_Q{H}_Q}{{\ell }_{PQ}}d{s}_Q$$

(4)

With H the thickness of the mass element and with density $\rho =\frac{dm}{Hds}=\frac{dm}{d\sigma }$ and differential area element $d{s}_Q=\sin {\theta }_{Q}d\theta d\lambda$. When applied to the topographic masses, we obtain the disturbing potential $T_P$, and further considering Bruns equation [49], produce the effect for geoid undulation $d{N}_P$.

$$d{N}_P=\frac{T_P}{\gamma }=\frac{G}{\gamma }\iint \frac{{\rho }_Q{H}_Q}{{\ell }_{PQ}}d{s}_Q$$

(5)

where γ is the normal gravity at the evaluation point. Equation (5) is called the condensation approximation [12]. In order to use a constant mass density value in Equation (5), we use the Rock Equivalent Topography concept (RET), see, for instance [16]. For ocean areas having the rock equivalent height $H^{*}$ is given by

$$H^{*}=H\left(1-{\rho }_w/\rho \right)$$

(6)

where ${\rho }_w=1030 \mathrm{kg}/{\mathrm{m}}^3$ is the ocean water mass density and $\rho =2670 \mathrm{kg}/{\mathrm{m}}^3$ is the standard topographic rock density. Now Equation (5) can be written as

$$d{N}_P=\frac{G\rho }{\gamma }\iint \frac{H_Q}{{\ell }_{PQ}}d{s}_Q$$

(7)

Thus, when integration point $Q$ is on the ocean $H_Q=H_Q^{*}$.

2.2. Mass Discretization

With respect to mass discretization, the tesseroid approach [39] is employed in this study, which is more efficient than the methods based on the classical prism technique. The method is not rigorous. However, since we are dealing with relatively small quantities, such as the contribution of topographic masses to geoid undulations, the errors are negligible. Moreover, since tesseroids are bodies bounded by surfaces of constant height and geographical grid lines or curves of parallels and meridians over the terrestrial sphere or the ellipsoid, they are easily implemented with fast Fourier techniques, such as the 1D-FFT method implemented for 2D integration over the sphere [50]. This method, which combines 1D-FFT along parallels with direct numerical integration along meridians, is as rigorous as direct numerical integration but can handle a lot more data, being also more efficient. On the other hand, is not as efficient as the 2D FFT method but is not affected by the meridian convergence error. Equation (7), can be written as the spherical convolution integral.

$$g_P=\iint k\left({\psi }_{PQ}\right)F_{Q}d{s}_Q$$

(8)

According to [50], all the values of $g$ along the parallel of latitude ${\varphi }_P$ are given by

$$g_{\varphi P}\left({\lambda }_P\right)=\mathsf{\Delta }\varphi \mathsf{\Delta }\lambda {\hat{F}}_1^{-1}\left({\displaystyle \sum }_{{\varphi }_q={\varphi }_1}^{{\varphi }_n}{\hat{F}}_1\left(k\left(\mathsf{\Delta }{\lambda }_{PQ}\right)\right){\hat{F}}_1\left(f_{\varphi Q}\left({\lambda }_Q\right)cos{\varphi }_Q\right)\right)$$

(9)

where $\mathsf{\Delta }\varphi$, $\mathsf{\Delta }\lambda$ are the sampling intervals, ${\hat{F}}_1$ and ${\hat{F}}_1^{-1}$ represent the 1-D Fourier transform operator and its inverse, respectively, $g_{\varphi P}\left({\lambda }_P\right)$ is the $1\times m$ vector along the parallel with latitude ${\varphi }_P$, $f_{\varphi Q}$ is the $1\times m$ vector of $f$ along the parallel with latitude ${\varphi }_Q$, $k_{{\varphi }_Q}\left(\mathsf{\Delta }{\lambda }_{PQ}\right)$ is the $1\times m$ vector of the kernel function $k$ relating points of parallels ${\varphi }_P$ and ${\varphi }_Q$, n is the number of parallels, and m is the number of meridians. After the discretization of the integral in Equation (8), with a constant sampling interval along parallels, Equation (9) yields the integral evaluation for all points of parallel with latitude ${\varphi }_P$.

2.3. Corrector Surface

In practice, the evaluation and/or validation of GGM geoid undulations ($N^{GGM})$ are made through a set of GNSS/leveling stations that produces independent geoid undulation $N^{GL}$ values. The associated differences $\mathsf{\Delta }N$ at those stations will be affected by datum inconsistencies, GNSS/leveling errors, DEM errors, and GGM commission errors as well (assuming the topographic mass contribution has been applied to them). These errors effects can be significantly reduced, in some cases to a centimeter level, by regression or parametric models. For instance, a polynomial regression model of order n, with coefficients $a_i,b_i$ and in terms of spherical coordinates ($\varphi ,\lambda$) can be expressed.

$$\begin{array}{cc}\mathsf{\Delta }N\left(\varphi ,\lambda \right)\hfill & =N^{GGM}\left(\varphi ,\lambda \right)-N^{GL}\left(\varphi ,\lambda \right)\hfill \\ & =a_0+a_1\varphi +b_1\lambda +a_2{\varphi }^2+b_2{\lambda }^2+a_3{\varphi }^3+b_3{\lambda }^3\dots +a_n{\varphi }^n+b_n{\lambda }^n\hfill \end{array}$$

(10)

See for instance [51], the coefficient in Equation (10) $a_0$ represent a bias. While $a_1,b_1$ are associated to the model’s bias and tilts. In the present study, only the first and second order regression models were analyzed. Since our focus is for local application, for which geoid surfaces tend to be fairly smooth, solutions of higher models become unfeasible. The Kriging method was also applied to generate the corrector surface only for comparison and interpolation purposes. For larger areas other models can be applied, like the one presented in [52].

Kriging is an estimator approach that make use of a priori variance–covariance values of a function, predicting its values at unsampled points by the computing a weighted average of the known values of the function in the neighborhood of the point [53,54,55]. Ordinary kriging is one of the methods more commonly applied in geosciences. This approach assumes the functions to be stationary. It produces minimum error variance, being a BLUE, best linear unbiased estimator. We have that, at an unsampled point, the unknown true value Z is determined by

$$\hat{z}={\displaystyle \sum }_{i=1}^n{w}_i{z}_i$$

(11)

where $w_i$ is the weight of the corresponding known function value $z_i$. The weights are obtained such as the prediction error variance is minimized, imposing the constrain $\sum _{i=1}^n{w}_i=1$. The set of weights that minimize the error variance under this constraint satisfies the n + 1 system of equations:

$$\begin{gathered} {\displaystyle \sum }_{i=1}^n{w}_i{\tilde{C}}_{ij}+\lambda ={\tilde{C}}_{j0 } \forall j=1,\dots , n \\ {\displaystyle \sum }_{i=1}^n{w}_i=1 \end{gathered}$$

(12)

Being ${\tilde{C}}_{ij}$ the covariance between the function at points i and j. ${\tilde{C}}_{i0}$ is the covariance between point i and the estimation point $\lambda$ is the Lagrange multiplier.

The system can be written in matrix format

$$\begin{gathered} CW=D \\ \left[\begin{array}{cc}{\tilde{C}}_{11}& \begin{array}{ccc}\cdots & {\tilde{C}}_{1n}& 1\end{array}\\ \begin{array}{}⋮& {\tilde{C}}_{n1}& 1\end{array}& \begin{array}{}\begin{array}{ccc}\ddots & ⋮& ⋮\end{array}& \begin{array}{ccc}\dots & {\tilde{C}}_{nn}& 1\end{array}& \begin{array}{ccc}\dots & 1& 0\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{w}_1\\ ⋮\\ \begin{array}{c}{w}_n\\ \lambda \end{array}\end{array}\right]=\left[\begin{array}{c}{\tilde{C}}_{10}\\ ⋮\\ \begin{array}{c}{\tilde{C}}_{n0}\\ 1\end{array}\end{array}\right] \end{gathered}$$

(13)

Matrix $C$ is of size (n + 1)(n + 1), while vectors $W$ and $D$ have n + 1 elements.

Then

$$W=C^{-1}D$$

(14)

The minimized error variance ${\hat{\sigma }}_R^2$ or mean square error MSE is given by

$$MSE={\hat{\sigma }}_R^2={\hat{\sigma }}^2-\sum w_i{\tilde{C}}_{i0}-\lambda ={\hat{\sigma }}^2-W^{T}D$$

(15)

This method assumes that the random variables have the same mean and variance. These assumptions allow the use of the variogram to obtain the variance and covariance values. Using the relation

$$\gamma \left(d_{ij}\right)={\hat{\sigma }}^2-{\tilde{C}}_{ij}$$

(16)

$\gamma \left(d_{ij}\right)$ is the variogram function, which depends on the distance $d_{ij}$ between points$i,j$, for $d_{ij}=0$, $\gamma \left(0\right)={\hat{\sigma }}^2$.

3. Study Area and Data Application

In Mexico, the National Bureau of Statistics Geography and History INEGI (Instituto Nacional de Geografia e Historia) for the last decade, has been launching campaigns to stablish GNSS/leveling stations all over the country. The main purpose of this effort is to increase the access of geodetic references for survey engineering projects, either for vertical and horizontal geo-referencing. The data set employed on this study, specifically the geodetic coordinates and orthometric heights of the GNSS/leveling stations that belong to the Mexican National Geodetic network, which is referenced to the Mexican National Geodetic Active Network. This network consists of 32 Continuously Operating Reference Stations (CORS). The precise geodetic positioning of the GNSS/leveling stations was made through measuring sessions involving data from the CORS stations ensuring 5 cm or better geodetic positioning accuracy. With respect to the orthometric heights, the leveling stations that are part of the Mexican national vertical control network were determined to employ first order second class geodetic leveling network standards. According to them, the accumulation error of the spirit leveling should be of 4 mm$\sqrt{\mathrm{km}}$ for 30–50 km circuits, maintaining vertical position precision of 2–3 cm (one $\sigma$). On some urban areas, the geodetic stations generate a local geodetic network. We have identified so far five cities with these networks: Durango, Hermosillo, Merida, Mexico City, and Monterrey (Figure 1).

The networks were of different number of points as well as area coverage. The one with smallest number of stations was Monterrey with 166, having a coverage area of about 13,000 square kms (1.2 square degrees), whereas Hermosillo with similar extension has a network of 337 stations. Durango and Merida networks have the smallest area coverage with 3000 and 1800 square kms (0.3 and 0.2 square degrees) respectively, but with corresponding 250 and 182 geodetic points. The Mexico City, or central Mexico area was the largest one, 57,000 square kms corresponding to about 5 square degrees, and with the largest number of stations (811) as well. The GGMs employed in this work are EGM2008 and XGM2019e_2159. They are high resolution geopotential models with spherical harmonic coefficient expansion of 2190 with respect to degree and order. According to ICGEM (http://icgem.gfz-potsdam.de) and from a comparison with 4898 GNSS/leveling points in Mexico, they have a precision of 21 and 17 cm, respectively. To get the corresponding RTM, which were going to be utilized to compute topographic masses contribution to the geoid, the STRM V4.1 with 1 arc-second (30 m) resolution was used as a DEM and as a reference field of topographic masses the TGFM models DTM2006.0 [22] and EARTH2014 [28] with normalized spherical harmonic coefficients were employed for EGM2008 and XGM2019e_2159, respectively. For comparison purposes, the most recent gravimetric geoid for Mexico GGM10 [48] was also analyzed. This geoid was generated from gravity data through the Stokes–Helmert remove restore technique. On its elaboration, the EIGEN-GRACE02S [56] satellite-only global geopotential model was employed for the long wavelength content of the field. The GGM10 estimated precision is of 20 cm.

4. Results and Analysis

As an initial analysis, an assessment was made about accuracy and consistency of the geoid models. For every local network, a comparison of geoid undulation (N) values, as produced by the geoid models EGM2008, XGM2019e_2159, and GGM10 with respect to the corresponding values N obtained from GNSS/leveling, was made. Some of those differences are shown on Figure 2 while

Table 1. Statistics of Geoid undulation differences between the geoid models (EGM2008, XGM2019e_2159, and GGM10), and GNSS/leveling stations for five Mexican urban areas. Units: Meters.

	GNSS/LEV-EGM2008				GNSS/LEV-XGM2019_2159				GNSS/LEV-GGM10
AREA	MAX.	MIN.	MEAN	STANDARD DEVIATION	MAX.	MIN.	MEAN	STANDARD DEVIATION	MAX.	MIN.	MEAN	STANDARD DEVIATION
DURANGO	0.752	0.303	0.525	0.063	0.708	0.123	0.440	0.091	0.813	0.288	0.585	0.064
HERMOSILLO	−0.004	−0.460	−0.118	0.076	0.048	−0.455	−0.085	0.080	0.425	−0.247	0.263	0.094
MERIDA	−0.472	−0.765	−0.598	0.046	−0.444	−0.789	−0.615	0.048	−0.332	−0.632	−0.468	0.046
MONTERREY	0.180	−0.189	−0.017	0.066	0.198	−0.155	0.039	0.076	0.738	0.257	0.499	0.104
MEX CITY	0.800	0.100	0.487	0.169	0.858	−0.141	0.453	0.199	1.409	0.238	0.898	0.229

shows the corresponding statistics. From this table, by observing the mean values of the geoid differences, we can appreciate that all the models have local biases ranging from −60 cm of EGM2008 at Merida region to 90 cm of GGM10 at Mexico City. In this City, the other two models had a mean of about 50 cm. Moreover, over this city, the three geoid models had large maximums of the differences with values of 80, 86, and 141 cm for EGM2008, XGM2019e_2159, and GGM10, respectively. With respect to the minimums of geoid differences, the three models had the largest value at Merida City, ranging from −63 to −79 cm. With respect to the standard deviation, EGM2008 models performed relatively better than the other two models, having values of 17 cm at Mexico City and 5–8 cm for the rest of the working areas. XGM2019e_2159 had a standard deviation of 20 cm at Mexico City and 5–9 cm for the others cities, and finally the gravimetric geoid GGM10 had 23 cm for Mexico City and 5–10 cm for the rest of the cities. The smallest standard deviation values were obtained at Merida City with the three models, 5 cm.

In the next step, estimation of the topographic masses’ contributions to N for both geopotential models EGM2008, and XGM2019e_2159 were computed for all the urban areas. The Residual Terrain model was formulated from the EDM STRM V4.1 model and the TGFM models DTM2006.0 (for EGM2008) and EARTH2014 (for XGM2019e_2159). A FORTRAN90 program was developed to compute the topographic masses contributions using the spherical tesseroids discretization method, and considering the REC method to account for mass densities. The program performs the numerical integration using 1D-FFT along parallels with integration along meridians, according to Equation (9). The geometrical configuration of the data employed consists of 7.5 seconds intervals and two square degree area coverage, with solution at one central square degree. Figure 3 exposes the surface plots while

Table 2. Statistics of topographic mass contributions to the geoid for 5 Mexican urban areas. Discretization method: spherical tesseroids. Integration method 2D spherical integration with 1D-FFT. Units: Meters.

STAT	DUR	HILLO	MER	MTY	MEX CITY
MIN	−0.056	−0.030	−0.003	−0.117	−0.088
MAX	0.058	0.060	0.003	0.147	0.167
MEAN	0.001	−0.001	0.000	0.001	0.000
ST DEV	0.013	0.007	0.001	0.030	0.024

shows the corresponding statistics. From Figure 3, one can appreciate that Mexico City area has the most rugose surface followed by Monterrey City. This is probably due to the fact that they are in or near mountainous regions. It can also be seen from Table 2 that the contributions are relatively higher for these two cities with maximum absolute values around 15 cm. For the other three cities, the topographic masses contributions to the geoid where generally smaller than a centimeter with maximum absolute values of about 6 cm: Similarly, for the standard deviations, in Table 2, Monterrey and Mexico Cities show the larger values of about 3 cm, while the other three cities have values of 13 mm or less.

For the final numerical analysis of this work, corrector surfaces were generated from geoid heights differences $\mathsf{\Delta }N$ obtained at the geodetic networks between N values from GNSS/leveling and the ones from the geoid models EGM2008, XGM2019e_2159, and GGM10. For the case of the two GGMs employed, the topographic mass contributions were previously applied to the corresponding $N^{GGM}$ values. The corrector surfaces were estimated applying polynomial regression models, according to Equation (10), and the ordinary kriging interpolation method, according to Equations (11)–(16), to the $\mathsf{\Delta }N$ values. With respect to the polynomial regression method, on all the study areas only the first and second expansion were able to be estimated (n = 1, n = 2). This is most likely due to the geoid smoothness at a local/regional scale. All the residuals obtained were statistically similar for the three geoid models. The two polynomial regression models yielded very similar results, with the n = 2 polynomial slightly better with differences up to 1 cm on the standard deviation. These two models produced residuals with standard deviations of about 15 cm for Mexico City and 4–8 cm for the other areas (

Table 3. Statistics of the corrector surfaces obtained from applying polynomial regressions of order 1, 2, and kriging interpolation to the geoid undulation differences between EGM2008 and GNSS/leveling stations for five Mexican urban areas. Units: Meters.

DURANGO
	EGM2008			GGM10			XGM2019
	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING
MAX	0.214	0.215	0.182	0.228	0.213	0.182	0.296	0.224	0.178
MIN	−0.239	−0.213	−0.230	−0.263	−0.239	−0.195	−0.263	−0.228	−0.192
MEAN	0.000	0.000	−0.060	0.000	0.000	0.000	0.000	0.000	−0.001
ST DEV	0.061	0.059	0.058	0.061	0.060	0.050	0.071	0.063	0.049
HERMOSILLO
	EGM2008			GGM10			XGM2019
	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING
MAX	0.123	0.200	0.209	0.122	0.189	0.205	0.124	0.201	0.204
MIN	−0.313	−0.296	−0.220	−0.340	−0.301	−0.225	−0.332	−0.282	−0.217
MEAN	0.000	0.000	0.001	0.000	0.000	0.001	0.000	0.000	0.001
ST DEV	0.075	0.066	0.048	0.079	0.066	0.049	0.078	0.068	0.048
MERIDA
	EGM2008			GGM10			XGM2019
	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING
MAX	0.156	0.178	0.034	0.164	0.184	0.171	0.187	0.188	0.167
MIN	−0.156	−0.154	−0.271	−0.160	−0.157	−0.144	−0.173	−0.160	−0.141
MEAN	0.000	0.000	−0.130	0.000	0.000	0.000	0.000	0.000	0.000
ST DEV	0.044	0.042	0.036	0.044	0.042	0.035	0.047	0.046	0.035
MONTERREY
	EGM2008			GGM10			XGM2019
	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING
MAX	0.203	0.186	0.172	0.189	0.165	0.151	0.181	0.184	0.169
MIN	−0.155	−0.170	−0.163	−0.149	−0.164	−0.175	−0.176	−0.174	−0.174
MEAN	0.000	0.000	0.002	0.000	0.000	0.001	0.000	0.000	0.000
ST DEV	0.065	0.056	0.054	0.070	0.062	0.053	0.057	0.056	0.054
MEXICO CITY
	EGM2008			GGM10			XGM2019
	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING	POLY-1	POLY-2	KRIGING
MAX	0.339	0.318	0.308	0.469	0.453	0.315	0.457	0.413	0.311
MIN	−0.539	−0.452	−0.413	−0.497	−0.491	−0.396	−0.551	−0.476	−0.397
MEAN	0.000	0.000	0.001	0.000	0.000	0.000	0.000	0.000	0.000
ST DEV	0.140	0.137	0.086	0.165	0.162	0.083	0.150	0.141	0.085

). For all the geodetic networks, the kriging method gave the optimum solution, producing standard deviations for the residuals of 8 cm for Mexico City and 4–6 cm for the other study areas. With the formulation of the required variogram, the spherical model was employed. For every data set, this variance model was adjusted in order to be consistent with the data variance. For every city it was found that the variance of all the data sets were very similar. The data variances of the cities ranged from 0.0035 ${\mathrm{m}}^2$ at Durango City to 0.018 ${\mathrm{m}}^2.$ For the case of the variogram nugget effect and range parameter, they went respectively from 0.001 and 0.2 at Merida City to 0.01 and 1.0 in Mexico City.

5. Discussion and Conclusions

An analysis was made for the accuracy improvement of geoid undulation values as produced by a high expansion Global Geopotential Model such as EGM2008 and XGM2019e_2159. The applied procedure involves the use of a Digital Elevation Model, a Topographic Reference Field, and a set of GNSS/leveling stations. The scope of the study was local or regional and the areas of application were five urban areas of Mexico. As expected, the masses’ contributions to the geoid signal were found to be relatively small, generally smaller, or equal, to a centimeter, except with the Mexico and Monterrey Cities cases that showed standard deviations of about 2 cm and maximum absolute values of 17 and 15 cm, respectively. This was probably produced by the proximity to the Popocatepetl volcano area which is very mountainous. On the comparison of geoid undulation values over the GNSS/leveling stations between the GGM and the ones of the GNSS/leveling stations, local biases with the range −60 to 90 cm were observed. This was also the case with the Mexican gravimetric geoid GGM10, which was used for comparison purposes. Disregarding these biases and according to the standard deviations of geoidal differences, EGM2008 could be considered to have 17 cm precision on Mexico City area and 5–8 cm precision over the other urban areas, showing better performance than the one reported at ICGEM website. At the same time XGM2019e_2159 and GGM10 would have, respectively: 20 and 23 cm precision for Mexico City and 5–9 and 5–10 cm precision for the other areas.

The corrector surfaces determined from the geoid heights differences between the ones from GSSN/leveling and the ones from the three geoid models, produced very similar values. When generating the corrector surfaces with polynomial regressions of order one and two from the geoid differences, the improvement of N produced by the GGM was about 60% in most of the cases, the kriging method being the one with the optimum results for all the cases. From the results obtained on the present work, one can say the geoid produced by the applied procedure is of 8 cm for the Mexico City area and of 4–6 cm for the rest of the analyzed urban areas. Therefore, a sub-decimeter geoid is possible applying the present methodology. On the other hand, we cannot draw any conclusions about the relationship of the performance or quality of the corrector surfaces and the GNSS-leveling networks configuration (extension, distribution, and number of stations), since all of them treated on this study are very dissimilar, having differences about topographic roughness as well. Possible future research could focus on broader areas or at a national scope as the Mexican government keeps producing more geodetic data, specifically establishing more GNSS/leveling stations all over the country.