Quasi Geoid and Geoid Modeling with the Use of Terrestrial and Airborne Gravity Data by the GGI Method—A Case Study in the Mountainous Area of Colorado

This article concerns the development of gravimetric quasigeoid and geoid models using the geophysical gravity data inversion technique (the GGI method). This research work was carried out on the basis of the data used in the Colorado geoid experiment, and the mean quasigeoid (ζm) and mean geoid (Nm) heights, determined by the approaches used in the Colorado geoid experiment, were used as a reference. Three versions of the quasigeoid GGI models depending on gravity data were analyzed: terrestrial-only, airborne-only, and combined (using airborne and terrestrial datasets). For the combined version, which was the most accurate, a model in the form of a 1′×1′ grid was calculated in the same area as the models determined in the Colorado geoid experiment. For the same grid, the geoid–quasigeoid separation was determined, which was used to build the geoid model. The agreement (in terms of the standard deviation of the differences) of the determined models, with ζm and Nm values for the GSVS17 profile points, was ±0.9 cm for the quasigeoid and ±1.2 cm for the geoid model. The analogous values, determined on the basis of all 1′×1′ grid points, were ±2.3 cm and ±2.6 cm for the quasigeoid and geoid models, respectively.


Introduction
Regional gravimetric geoid and quasigeoid models are mainly built on the basis of Stokes and Molodensky integrals and the least squares collocation (LSC) method. A general basis for these approaches can be found in physical geodesy textbooks, e.g., [1][2][3]. It should be emphasized that the above-mentioned general solutions have many modifications that significantly improve the quality of the models determined and enable the use of a global geopotential model (GGM) and digital elevation model (DEM). These additional data are most often used in the remove-compute-restore (RCR) procedure, in which the calculation runs in three steps. In the first step, a global model part (long-wavelength part) and a part based on the DEM (short-wavelength part) are removed from the original gravity data, determining the residual gravity values. In the second step, the residual gravity data are transformed into residual geoid or quasigeoid heights by the mentioned techniques. Finally, in step three, the GGM and DEM components of geoid undulations or height anomalies are determined and added to their residual values, giving the final results. Among the ways that the RCR technique is implemented in relation to the use of topography in steps two and three, the Helmert method of condensation, e.g., [4][5][6][7] and the residual terrain model method, e.g., [8][9][10][11] are the most common. A modification to this approach is the window remove-restore technique [12].
Step two is often carried out with a fast Fourier technique, e.g., [13][14][15][16][17]. An alternative to the RCR method is the method developed at the

The Utilized Solution: Overview
In the GGI approach, the disturbing potential in the elaboration area is represented by a three component model [34,43].
The first component covers the potential of local topographic masses included in a volume Ω, usually denoted as T Ω . The second component (T κ ) is the potential of disturbing masses located between the geoid and the compensation surface, included in the volume κ, corresponding in its horizontal ranges to the volume Ω ( Figure 1). The first component covers the potential of local topographic masses included in a volume , usually denoted as . The second component ( ) is the potential of disturbing masses located between the geoid and the compensation surface, included in the volume , corresponding in its horizontal ranges to the volume ( Figure 1). The components and are given by Newton's integrals, where and are density distribution functions in defined volumes and , respectively, r is the distance between the attracting masses and the attracted point P located on the Earth's surface, G is Newton's gravitational constant, and and are elements of volumes.
Both volumes ( and ) are horizontally limited, slightly beyond the study area; thus, they do not cover all the disturbing masses with a local disturbing potential. This makes it necessary to introduce an additional potential , which represents the potential of the disturbing masses not covered by the and components. The potential also covers the systematic and long-wavelength errors of the used data. This assumes its low variability, so it can be represented in the form of low degree harmonic polynomials: where , , are the coordinates of the point P and , … , are polynomials coefficients.
In general, the disturbing potential at a point on the terrain surface will be written as: In the model presented above, the unknown parameters are the 3D density distribution functions and , as well as the coefficients of the polynomials defining the external potential. In the original concept of the method, these parameters are determined by the least squares method, on the basis of a dense network of gravity points with known gravity anomalies or gravity disturbances and a sparse network of GNSS/levelling points with known height anomalies, converted into disturbing potential values. The gravimetric solution requires the model parameters to be determined without the use of GNSS/levelling data. The proposed process for their omission will be described at the end of this section.
The determination of the density functions and is carried out by the procedure of the linear inversion of gravity data [44], which requires their discretization. Thus, the The components T Ω and T κ are given by Newton's integrals, where ρ and δ are density distribution functions in defined volumes Ω and κ, respectively, r is the distance between the attracting masses and the attracted point P located on the Earth's surface, G is Newton's gravitational constant, and dV Ω and dV κ are elements of volumes. Both volumes (Ω and κ) are horizontally limited, slightly beyond the study area; thus, they do not cover all the disturbing masses with a local disturbing potential. This makes it necessary to introduce an additional potential T E , which represents the potential of the disturbing masses not covered by the T Ω and T κ components. The T E potential also covers the systematic and long-wavelength errors of the used data. This assumes its low variability, so it can be represented in the form of low degree harmonic polynomials: where X P , Y P, H P are the coordinates of the point P and a 1 , . . . a 5 , are polynomials coefficients. In general, the disturbing potential at a point on the terrain surface will be written as: In the model presented above, the unknown parameters are the 3D density distribution functions ρ and δ, as well as the coefficients of the polynomials defining the external potential. In the original concept of the method, these parameters are determined by the least squares method, on the basis of a dense network of gravity points with known gravity anomalies or gravity disturbances and a sparse network of GNSS/levelling points with known height anomalies, converted into disturbing potential values. The gravimetric solution requires the model parameters to be determined without the use of GNSS/levelling data. The proposed process for their omission will be described at the end of this section.
The determination of the density functions ρ and δ is carried out by the procedure of the linear inversion of gravity data [44], which requires their discretization. Thus, the Ω and κ volumes are divided into finite volume blocks, and a constant density is assigned to each of the blocks. In previous and current studies, a very simple division was used, taking into account only the lateral density variation in two layers: one for the volume Ω, and one for the volume κ. Consider that the volume Ω is defined by blocks of the digital elevation model (DEM), which are grouped into zones of constant, determined densities. The κ volume consists of constant density blocks, extending from the geoid to the compensation surface. The horizontal size of these blocks corresponds to constant density zones of the volume Ω.
Due to the local nature of the GGI models built so far, the calculations were carried out in the local Cartesian coordinate system. The Z-axis of the coordinate system is directed towards the geodetic zenith, at the origin point located in the middle of the study area. The X and Y axes lie on the plane of the horizon and are directed toward the north and east, respectively. In such a defined coordinate system, the Ω and κ volumes are represented in the form of rectangular prisms, for which the solutions of Newton's integrals for gravity potential (Equations (1) and (2)) and their derivatives are well-known, e.g., [45,46].
With the Ω and κ volumes defined in this way, Equations (1) and (2) can be written as follows [41]: where n is the number of constant density zones of the DEM; ρ k is the searched constant density of zone k; m k is the number of rectangular prisms of the DEM in zone k; s is the number of rectangular prisms defining the κ volume; and δ j is the searched density of the rectangular prism j. The K i and K j coefficients are solutions of Newton's integrals for the rectangular prisms: i of the DEM and j of the κ volume: y j1 x j2 x j1 1 r j dx j dy j dz j (8) where x i1 , x i2 , y i1 , y i2 , z i1 , z i2 are the coordinates defining the rectangular prism i of the x j1 , x j2 , y j1 , y j2 , z j1 , z j2 are the coordinates defining the rectangular prism j of the κ volume; and It should be noted that the implementation of Equations (7) and (8) takes into account the sphericity of the Earth. Therefore, for each attracted point P we assume Z P = H P , while the height coordinates defining DEM blocks and prisms of the κ volume take into account where d is the horizontal distance of the point P and the DEM block center, and R is the mean of the Earth's radius.
In calculations, a certain reference density model is usually used, wherein a constant value is assumed for the volume Ω, denoted as ρ 0 . Based on this, the reference density model for the volume κ (δ 0 ) is adopted as the density, which balances topographical masses of the volume Ω. Assuming that the constant density zone i of the Ω volume is located exactly above the constant density block j of the volume κ, the reference density of this block δ 0 j is defined by the following equation: where H i is the mean height of zone i of the Ω volume and h i is the height of block j of the volume κ. Assuming that the determined density model would be written as: And a reference density model is written in the form: The vector of unknowns can be written as: where dτ = τ − τ 0 is the determined vector of residual density and a = [a 1 , . . . , a 5 ] T . Assuming the vector of determined parameters defined by Equation (12), the observation equations for the disturbing potential (T) and gravity disturbance (δg) at point P will have the following form [43]: where f and f z are the vectors of known parameters resulting from the equation of the disturbing potential model (4) and its Z P derivative respectively, and v T and v δg are adjustment errors.
The approximate observation quantities are determined based on the vector δg 0 = −f T z x 0 For the series of observations, a system of Equations (13) is built, which we can write as: where v T = v T . . . , v δg . . . is the vector of adjustment errors, A is the design matrix of known coefficients, and L T = [T − T 0 , . . . , δg − δg 0 , . . . ] is the vector of observations. The use of more than one layer of the determined density distribution leads to ambiguity in the solution of the gravity inversion problem. This ambiguity is solved by introducing a deep weighting function and condition proposed by [47].
For all assigned unknowns the condition can be written as: where W x = W a 0 0 W τ , W a is the zero weighting matrix assigned to the vector a and W τ is the density model weighting matrix, defined based on [47].
A detailed description of the definition and determination of the deep weighting function parameters that were used can be found in [43].
With regard to condition (16), with a defined weight observation matrix P, the least squares objective function is written as: This condition leads to the solution of equation system 15 in the following form: Note that the use of a deep weighting function (matrix W τ ), not only solves the problem of solution ambiguity, but also allows for some control over the depth of the location of the designated densities. The second factor influencing density location is related to the utilized observations and can be controlled, to some extent, by the weights given to individual groups of them. Note that in the model under discussion, the gravity data are more sensitive to density variation located closer to the terrain surface, and the potential values are more sensitive to density variation in deeper locations. Hence, it is important to properly design the weight matrix P. In previous studies, it has been assumed that the weights for the disturbing potential values are closely related to their accuracy (the accuracy of the GNSS/levelling data). Regarding the gravity values, it has been assumed that their weights are related, not to the accuracy of their determination, but more to the accuracy of the GGI model in terms of gravity. For the calculations, a value 2-3 times greater than the prediction accuracy has been assumed as the basis for estimating the weights of gravity values.
Calculations can also be performed with the use of global geopotential models. In this case, the RCR procedure is used. If we calculate gravity potential values (W) for GNSS/levelling points, the RCR procedure may be carried out according to the following scheme:

1.
Removal from the original data of a global model part (long-wavelength part) and determination of the residual disturbing potential values δT r and the residual gravity δg r .
where T GM , W GM , and g GM are the disturbing potential, gravity potential, and gravity from the global model, while g is the measured gravity.

2.
On the basis of the residual data, the residual model of the disturbing gravity potential is built using the GGI method.
where the components δT E , δT Ω , and δT κ are residual parts of the componenets T E , T Ω , and T κ (after removing from them parts contained in the global geopotential model).

3.
For new points, the residual disturbing potential values from the GGI model and the gravity potential values from the global model are determined. Consequently, the values of gravity potential (W GGI ) and gravity (g GGI ) can be determined.
where δg r GGI is the residual gravity from the GGI model.
On the basis of the gravity potential values, it is possible to determine, for example, normal heights H N GGI or height anomalies (ζ GGI ) using known equations (Heiskanen and Moritz 1967): where W o is the gravity potential on the geoid, γ is the mean value of the normal gravity between the ellipsoid and telluroid, and γ Q is the normal gravity on the telluroid. The observation Equations (13) for residual values can be written as The values T 0 and δg 0 , calculated on the basis of the reference density model (14), remain unchanged.
Since determination of the gravimetric geoid and quasigeoid models excludes the use of GNSS/levelling data, these data have to be replaced in the calculation procedure. A natural solution to the problem is to replace this data with values determined directly from the global model. In this case, the entire solution remains unchanged. We assume equality, W = W GM , so the Equation (25) will take the following form: The advantage of this solution is that we can include in the calculations a regular grid of points with known disturbing potential values, the size and density of which can be freely defined. On the other hand, these data are mostly less accurate than the GNSS/levelling data, which may have a negative impact on the accuracy of the model.
Notably, as a result of the solution discussed, the gravimetric GGI model will be a model locally fitted to the global model used.
Geoid to quasigeoid separation can be determined, for example, based on an equation, e.g., in [2]: where H and H n are the orthometric and normal heights, respectively, g is the mean gravity between the geoid and the terrain surface, and γ is the mean normal gravity between the ellipsoid and telluroid. This value can be used for geoid model determination using the above-described solution. Therefore, we finally determine the geoid model using the following relation: The presented solutions for geoid and quasigeoid determination require supplementation resulting from the specific conditions applied for the analyzed area: The determined quasigeoid model should be consistent with the IHRS reference level W o = 62, 636, 853.400 m 2 /s 2 [48], which is different from the normal potential on the reference ellipsoid surface GRS80, U o = 62, 636, 860.850 m 2 /s 2 . Hence, determining the height anomaly based on the gravity potential value (24), the difference ∆W o = −7.45 m 2 /s 2 has been taken into account.

2.
Since the GM used is in the zero-tide system and the models developed in the Colorado geoid experiment were developed in the non-tidal system, for the determined height anomalies, a small correction dζ zn was added [49].
where ϕ is the point latitude, k ≈ 0.3 is the Love number, and the result is in meters. Finally, the height anomalies were determined by the following equation: 3.
It should be added that the mean gravity g (used in Equation (27)) should be consistent with the orthometric height type used. As in the analyzed area, Helmert orthometric heights are used, the mean gravity g in our analyses were calculated in accordance with these heights, according to the following equation [2]: which assumes a density of the topography ρ = 2670 kg/m 3 .

The Used Data
The data provided by the US National Geodetic Survey included several datasets and were the same as used in the Colorado geoid computation experiment [28]. The data covered the area between the parallels 35 • and 40 • North (latitude) and the meridians 110 • and 102 • East (longitude) and included the following: A detailed description of these datasets can be found in [28].
The GSVS17 data set consists of points forming a profile with a length of over 350 km, in which, among others, Helmert geoid undulations, the deflection of the vertical components, and gravity values were determined. These points were used to compare the individual geoid and quasigeoid models developed in the Colorado geoid experiment, wherein the average value of the 13 developed models was used as a reference for the evaluation of the individual models [28]. These averaged values were also used in our analyses, to evaluate and compare the developed models using the GGI method.
In the calculations for the primary data, we used 55,161 terrestrial gravity points and 24,030 airborne gravity values. The location of these points is marked in Figure 2a. The Helmert orthometric heights of terrestrial gravity points were converted to the normal heights used in the calculations, while the ellipsoidal heights of these points were determined using the height anomalies from the global model XGM2016 [50], which were also used in the calculations as a reference GM in the RCR procedure. The height anomalies from the XGM2016 model were also used to determine the normal heights of airborne gravity points. 3. It should be added that the mean gravity ̅ (used in Equation (27)) should be consistent with the orthometric height type used. As in the analyzed area, Helmert orthometric heights are used, the mean gravity ̅ in our analyses were calculated in accordance with these heights, according to the following equation [2]: which assumes a density of the topography = 2670 ⁄ .

The Used Data
The data provided by the US National Geodetic Survey included several datasets and were the same as used in the Colorado geoid computation experiment [28]. The data covered the area between the parallels 35° and 40° North (latitude) and the meridians 110° and 102° East (longitude) and included the following: GSVS17 GPS/levelling data set (223 points) A detailed description of these datasets can be found in [28]. The GSVS17 data set consists of points forming a profile with a length of over 350 , in which, among others, Helmert geoid undulations, the deflection of the vertical components, and gravity values were determined. These points were used to compare the individual geoid and quasigeoid models developed in the Colorado geoid experiment, wherein the average value of the 13 developed models was used as a reference for the evaluation of the individual models [28]. These averaged values were also used in our analyses, to evaluate and compare the developed models using the GGI method.
In the calculations for the primary data, we used 55,161 terrestrial gravity points and 24,030 airborne gravity values. The location of these points is marked in Figure 2a. The Helmert orthometric heights of terrestrial gravity points were converted to the normal heights used in the calculations, while the ellipsoidal heights of these points were determined using the height anomalies from the global model XGM2016 [50], which were also used in the calculations as a reference GM in the RCR procedure. The height anomalies from the XGM2016 model were also used to determine the normal heights of airborne gravity points.   the analyzed area is presented in Figure 2b. The DEM shown in Figure 2b also defines the horizontal range of the Ω and κ volumes used. In accordance with the adopted strategy of developing a gravimetric quasigeoid model, 910 points with known disturbing potential values from the GM were accepted into the calculations. These points evenly cover the entire elaboration area and are marked in Figure 2b with white dots. The red lines in Figure 2 show the location of the GSVS17 profile points and the black rectangles define the area for which the geoid and quasigeoid models with a 1 × 1 grid were determined.

Results
The calculation process included two stages, which led to the development of the final quasigeoid and geoid models. The first, preliminary stage, covered the analyses indicating parameters important for the modeling results, which were used in the second stage, in which the final models were developed.

Preliminary Phase Analyses
In the preliminary phase, the accuracy of the XGM2016 model, in terms of height anomalies, and the accuracy of the gravity prediction by the GGI method were estimated. These values are necessary for the correct estimation of the observation weights. To estimate the accuracy of gravity prediction, initial modeling was performed in a smaller, central part of the elaboration area for various ρ 0 values, assuming that the gravity errors adopted were at the level of ±3 mGal (1 mGal = 10 µm/s 2 ). Based on this, the accuracy of the gravity prediction was estimated to be at the level of ±2.2 mGal, which, as stated above, indicates that the gravity errors used to estimate the gravity weights should be approximately in the range ±(4 − 6) mGal. Therefore, the value of ±5 mGal was adopted for further analyses. The height anomalies determined from the XGM2016 model were compared along the GSVS17 profile; with the mean height anomalies obtained from the 13 models participating in the Colorado geoid experiment hereinafter marked as ζ Re f [28]. The obtained standard deviation of the differences, at the level of ±10 cm, was adopted for the weight estimation used in further calculations. So far, the GGI method has been applied to land plains and low mountains. For such areas, the significance of the adopted reference density model for the accuracy of the GGI model in terms of gravity values was demonstrated [41]. The currently analyzed area covers high mountains (average height over 2200 m and with the highest point above the height of 4000 m) and the disturbing potential values used have lower weights, so their impact on the estimated densities will be smaller. Therefore, we can also expect a significant impact of the reference density model on the accuracy of the model, in terms of the height anomalies. To analyze this issue, a series of test calculations were performed, covering a wide range of reference densities ρ 0 (from 0 to 2670 kg/m 3 in increments between 100 and 500 kg/m 3 ). We assumed that the size of the constant density zones was 6 × 6 km. The calculations were performed for three versions of the gravity data used: 1.
combined (with the use of terrestrial and airborne gravity data).
Height anomalies determined from the GGI model (ζ GGI ) for each version of the gravity data and for different ρ 0 values were compared to the reference values (ζ m ) for the GSVS17 points. The standard deviations of the differences, ∆ζ = ζ GGI − ζ m , were next determined as the model accuracy parameters. Relations between these standard deviations and the ρ 0 values are presented in Figure 3 by solid lines. For each analyzed version of the reference density model, the mean density of the Ω volume (ρ mean ) from the GGI model was also determined. This allowed us to state the value of the mean density changes ∆ρ 0 = ρ mean − ρ 0 for each modeling cycle. These values are represented in Figure 3 by dotted lines. Individual versions of the gravity data used are marked with different colors: airborne-only = green, terrestrial-only = black, and combined = red.
The results presented in Figure 3 show a clear relationship between the reference constant density model and the accuracy of the height anomalies determined by the GGI method. The least accurate are the results for = 0 ⁄ . There is a noticeable increase in the accuracy with the increase in the reference density values up to an optimal range of , for which the accuracy is the highest. This optimal range is the narrowest for the terrestrial version and significantly wider for the two other versions The optimal ranges of values for individual versions are presented in Table 1. The ranges were adopted as the ranges, giving an accuracy no worse than 0.2 from the best results. The shown strong dependence of the GGI model accuracy on the values causes a certain problem for the estimation of the optimal density ( ) giving the highest quasigeoid accuracy. The best solution of this problem would be the use of a set of accurate reference GNSS/levelling data for the estimation of the value. In the absence of such data, it is necessary to rely on the experience in various implementations of the method. Taking into account the wide range of optimal values for the airborne-only and combined versions, it can be expected that for these versions, it will not be very difficult to estimate the value based on experience. For the terrestrial-only version, however, it will be much more difficult. In this case, the range of optimal values of is relatively narrow. Moreover, beyond this range, the decrease in accuracy is significant (compared to the versions using airborne gravity data). Hence, some additional parameters supporting the determination of the value should be indicated. In this regard, the Δ values defined above may turn out to be helpful. The relationships between The results presented in Figure 3 show a clear relationship between the reference constant density model and the accuracy of the height anomalies determined by the GGI method. The least accurate are the results for ρ 0 = 0 kg/m 3 . There is a noticeable increase in the accuracy with the increase in the reference density values up to an optimal range of ρ 0 , for which the accuracy is the highest. This optimal range is the narrowest for the terrestrial version and significantly wider for the two other versions The optimal ranges of ρ 0 values for individual versions are presented in Table 1. The ranges were adopted as the ranges, giving an accuracy no worse than 0.2 cm from the best results. The shown strong dependence of the GGI model accuracy on the ρ 0 values causes a certain problem for the estimation of the optimal density (ρ 0 optimal ) giving the highest quasigeoid accuracy. The best solution of this problem would be the use of a set of accurate reference GNSS/levelling data for the estimation of the ρ 0 optimal value. In the absence of such data, it is necessary to rely on the experience in various implementations of the method. Taking into account the wide range of optimal ρ 0 values for the airborne-only and combined versions, it can be expected that for these versions, it will not be very difficult to estimate the ρ 0 optimal value based on experience. For the terrestrial-only version, however, it will be much more difficult. In this case, the range of optimal values of ρ 0 is relatively narrow. Moreover, beyond this range, the decrease in accuracy is significant (compared to the versions using airborne gravity data). Hence, some additional parameters supporting the determination of the ρ 0 optimal value should be indicated. In this regard, the ∆ρ 0 values defined above may turn out to be helpful. The relationships between ∆ρ 0 and ρ 0 values are represented in Figure 3 by dotted lines. Note that these relationships are linear for all versions. Importantly, for versions using terrestrial gravity data (terrestrial-only and combined), the values of ∆ρ 0 achieve 0 for a ρ 0 value inside or close to the optimal range. These values of reference density (ρ 0 zero ), for individual versions, are marked in Figure 3 by arrows and can be used to estimate the ρ 0 optimal value. For the combined version, the ρ 0 zero value is in the middle of the optimal range, so we can assume ρ 0 optimal = ρ 0 zero . For the terrestrial-only version, there is a slight but visible shift of the ρ 0 zero value from ρ 0 optimal . This shift shows that the use of the ρ 0 zero value to estimate the ρ 0 optimal value requires a further, broader and more detailed analysis. However, the presented results give hope that in the terrestrial-only or combined versions, this will be an effective method.
The above analyses were performed by comparing the modeling results with the reference data on the GSVS17 profile points. The best results were obtained for the combined version; hence, the final quasigeoid model was determined for this data version. The ρ 0 optimal value, used in further calculations, for this version was computed based on the ρ 0 zero value, assuming the equality of both quantities. As the changes of ∆ρ 0 are linear, in order to determine the ρ 0 zero value, it is sufficient to find the value of ∆ρ 0 for two ρ 0 , e.g., for minimal ρ 0 min = 0 kg/m 3 and maximal ρ 0 max = 2700 kg/m 3 (∆ρ 0 max ) values, and then to determine the ρ 0 zero value by a simple relation:

Determination of GGI Geoid and Quasi Geoid Models for Colorado Geoid Experiment Area
Geoid and quasigeoid models were determined according to the scheme presented in Figure 4.  On the basis of the preliminary phase analyses, the gravity quasigeoid model was determined for the area used in the methods that participated in the Colorado geoid experiment. The data set presented in Figure 2 was used in the calculations. In this figure, the area for which the 1′ × 1′ grid of the quasigeoid model was determined is marked by the black line. In the calculations, zones of constant density with a size of 6 × 6 were assumed, and the = value was determined using Equation (32). A map of On the basis of the preliminary phase analyses, the gravity quasigeoid model was determined for the area used in the methods that participated in the Colorado geoid experiment. The data set presented in Figure 2 was used in the calculations. In this figure, the area for which the 1 × 1 grid of the quasigeoid model was determined is marked by the black line. In the calculations, zones of constant density with a size of 6 × 6 km were assumed, and the ρ 0 optimal = ρ 0 zero value was determined using Equation (32). A map of the determined quasigeoid model is presented in Figure 5. On the basis of the preliminary phase analyses, the gravity quasigeoid model was determined for the area used in the methods that participated in the Colorado geoid experiment. The data set presented in Figure 2 was used in the calculations. In this figure, the area for which the 1′ × 1′ grid of the quasigeoid model was determined is marked by the black line. In the calculations, zones of constant density with a size of 6 × 6 were assumed, and the = value was determined using Equation (32). A map of the determined quasigeoid model is presented in Figure 5.  The geoid model was determined on the basis of the quasigeoid model, based on relation (28). The g values necessary to calculate the geoid to quasigeoid separation (27) were calculated according to formula (31). The gravity values g at each 1 × 1 grid point were determined by the interpolation of the complete Bouguer anomalies (CBA) using the kriging method. The CBA values were calculated using the same DEM as in the quasigeoid modeling, taking into account topographic masses within a radius of 167 km and assuming their constant density ρ = 2670 kg/m 3 . Figure 6 shows the CBA values used for gravity prediction and Figure 7 presents the geoid to quasigeoid separations used for geoid model determination.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 The geoid model was determined on the basis of the quasigeoid model, based on relation (28). The ̅ values necessary to calculate the geoid to quasigeoid separation (27) were calculated according to formula (31). The gravity values at each 1′ × 1′ grid point were determined by the interpolation of the complete Bouguer anomalies (CBA) using the kriging method. The CBA values were calculated using the same DEM as in the quasigeoid modeling, taking into account topographic masses within a radius of 167 and assuming their constant density = 2670 ⁄ . Figure 6 shows the CBA values used for gravity prediction and Figure 7 presents the geoid to quasigeoid separations used for geoid model determination.       The map of the determined geoid model is presented in Figure 8.

Discussion
In order to evaluate the developed models, height anomalies (ζ GGI ) and geoid undulations (N GGI ) determined by the GGI method in the 1 × 1 grids were compared to the average values obtained from the models analyzed in the Colorado geoid experiment (13 for the quasigeoid and 14 for the geoid) at the GSVS17 points. These values for height anomalies are marked in the same way as before (ζ m ), and for geoid undulations they are hereafter referred as N m . Subsequently, the differences ∆ζ = ζ GGI − ζ m and ∆N = N GGI − N m were determined, the basic statistics of which are presented in Table 2. Table 2. Basic statistics of the differences ∆ζ = ζ GGI − ζ m , ∆N = N GGI − N m , as well as ∆ζ GNSS/lev = ζ GGI − ζ GNSS/lev and ∆N GNSS/lev = N GGI − N GNSS/lev , based on the GSVS17 profile points (in centimeters). In brackets are ranges of corresponding values obtained for methods participating in the Colorado geoid experiment (based on [28] Ellipsoidal and Helmert orthometric heights measured at the GSVS17 profile points were used for the independent evaluation of the developed models. Based on the data provided, we determined the measured geoid undulations (N GNSS/lev ). The gravity values measured at the GSVS17 profile points were used to convert the Helmert orthometric heights into the normal heights (according to Equation (27)), and consequently we calculate the GNSS/levelling height anomalies (ζ GNSS/lev ) by the difference between the ellipsoidal and normal heights. Those values were compared to ζ GGI and N GGI values. Hence, the differences ∆ζ GNSS/lev = ζ GGI − ζ GNSS/lev and ∆∆ GNSS/lev = N GGI − N GNSS/lev were calculated; the basic statistics of which are also presented in Table 2.
All of the differences, the statistics of which are presented in Table 2, are also presented in Figure 9. From ∆ζ GNSS/lev and ∆N GNSS/lev , offsets of 90 cm (associated with the longwavelength errors of the US vertical datum [51]) were removed. The differences ∆ζ = ζ GGI − ζ m , ∆N = N GGI − N m , ∆ζ GNSS/lev = ζ GGI − ζ GNSS/lev , ∆N GNSS/lev = N GGI − N GNSS/lev , and ellipsoidal heights at GSVS17 points. From the ∆ζ GNSS/lev and ∆N GNSS/lev , offsets of 90 cm were removed.
The determined geoid and quasigeoid models were also compared with the N m and ζ m values for all 1 × 1 grids points. The most important statistics of the differences ∆ζ = ζ GGI − ζ m and ∆N = N GGI − N m are presented in Table 3. Table 3. Basic statistics of the differences ∆ζ = ζ GGI − ζ m and ∆N = N GGI − N m based on all 1 × 1 grids points (in centimeters). In brackets are the ranges of corresponding values obtained for methods participating in the Colorado geoid experiment (based on [28] When assessing the determined geoid and quasigeoid models, it is useful to refer to the results obtained by the methods participating in the Colorado geoid experiment. On the basis of [28], we summarized the ranges of the corresponding values obtained for all methods taking part in the experiment. These values are presented in Table 2; Table 3 in brackets, next to the values determined for the GGI method. Accepting the standard deviation (STD) values of the differences between the quantities determined by the GGI method and their reference values as the basic parameters of the accuracy assessment, we can determine the relation of the GGI method to the methods participating in the Colorado geoid experiment. The STD values determined in reference to the mean geoid and quasigeoid models, calculated based on the GSVS17 profile points (Table 2), were ±0.9 cm for the quasigeoid and ±1.1 cm for the geoid model. Using the results of [28], we can notice that only in one quasigeoid model and one geoid model from the Colorado geoid experiment were lower STD values obtained. The analogous STD values calculated on the basis of all 1 × 1 grid points (Table 3) are worse; they are ±2.3 cm and ±2.6 cm for the quasigeoid and geoid, respectively. Additionally, the relationships with the methods participating in the Colorado geoid experiment are worse. Slightly lower STD values were obtained for seven quasigeoid models and four geoid models. An independent evaluation of the developed models is their comparison with the GNSS/levelling height anomalies and geoid undulations. The STD values determined in reference to these values at the GSVS17 profile points (Table 2) for the quasigeoid and geoid models were ±2.6 cm and ±2.5 cm, respectively. Among the Colorado geoid experiment models, lower STD values were obtained for three quasigeoid and three geoid models. It should also be noticed that there was a relatively small decrease or no decrease in the accuracy parameters of the determined geoid model in relation to the quasigeoid model. Taking into account the technique used for the determination of the geoid model, the accuracy of the obtained geoid-quasigeid separation should be assessed as high.
Note that the differences of both ∆ζ and ∆N, represented in Figure 9 by red solid and dotted lines, respectively, are within very narrow ranges and are close to each other. The graphs of the differences ∆ζ GNSS/lev and ∆N GNSS/lev (blue solid and dotted lines) are also interesting. There is a clear change in their values in their final sections (approximately 80 endpoints of the profile). The shift of the mean values for the points in the start and end sections is at the level of −5 cm and is also visible for other geoid and quasigeoid modeling methods participating in the Colorado geoid experiment [28]. This shift is the main reason for the significant increase in the STD values of the ∆ζ GNSS/lev and ∆N GNSS/lev differences.

Conclusions
The conducted analyses lead to several important GGI method conclusions. Primarily, the proposed use of the disturbing potential values determined from the GGM, instead of the GNSS/levelling data, turned out to be very effective. This allows for the determination of a local gravimetric quasigeoid model, fitted to the GGM used.
Test calculations were based on terrestrial and airborne gravity data. Such a data set made it possible to carry out analyses on the three versions of the gravity data used: terrestrial-only, airborne-only, and combined, which used the airborne and terrestrial datasets. The GGI models developed for each version turned out to be dependent on the reference density model adopted for the calculations. By comparing the determined quasigeoid models with the reference values obtained in the Colorado geoid experiment, it was possible to estimate the optimal density ranges giving the most accurate results. The narrowest and the most unfavorable range of optimal density values was obtained for the terrestrial-only version. Significantly wider ranges of the optimal values of ρ 0 (more favorable) were obtained for the version using airborne gravity data (airborne-only and combined). Moreover, the highest accuracy was achieved for the combined version (Table 1). Therefore, the supplementing of the data set with airborne gravity data leads to an improvement of the quality of geoid and quasigeoid models.
The demonstrated relationship between the accuracy of the GGI model and the ρ 0 value causes problems for determining an accurate quasigeoid model and indicates the need to determine or estimate the optimal ρ 0 as one of the most important parameters in the modeling process. At the present stage, the estimation of the optimal reference density model should be based on the GNSS/levelling control dataset. In the absence of such data, one could follow the results of previous applications of the method and the optimal ρ 0 estimation, the first of which are the presented results obtained by the Colorado geoid experiment. Analyses of changes in the mean density of topographic masses (∆ρ 0 = ρ mean − ρ 0 ) can also be helpful in searching for the optimal ρ 0 values. For the terrestrial-only and combined versions, these values achieve 0 for a ρ 0 value inside or close to the optimal ranges. Therefore, they can be used for the estimation of the optimal densities for these versions. However, it should be clearly emphasized that the study of these relationships requires further analysis, mainly on the use of GGMs with different resolutions and for other types of terrain.
In the analyzed case, these relations allowed determining the optimal density for the combined version, for which the model in the form of a 1 × 1 grid was calculated in the same area used for the models determined in the Colorado geoid experiment. For the same grid, the geoid-quasigeoid separation was determined, and this was used to build the geoid model. The standard deviations of the differences between the geoid and quasigeoid models determined by the GGI method and the mean models determined for the GSVS17 profile points were ±0.9 cm for the quasigeoid and ±1.2 cm for the geoid model. Analogous values, determined on the basis of all 1 × 1 grid points, were ±2.3 cm and ±2.6 cm for the quasigeoid and geoid models, respectively. The presented results indicate the high-quality of the determined quasigeoid and geoid models. A relatively small difference between the accuracy parameters of the geoid and the quasigeoid also proves that this was a very accurate determination of the geoid-quasigeoid separation.