Gravimetric geoid modeling from the combination of satellite gravity model, terrestrial and airborne gravity data: a case study in the mountainous area, Colorado

Constructing a high-precision and high-resolution gravimetric geoid model in the mountainous area is a quite challenging task because of the lack of terrestrial gravity observations, rough topography and the geological complexity. One way out is to use high-quality and well-distributed satellite and airborne gravity data to fill the gravity data gaps; thus, the proper combination of heterogeneous gravity datasets is critical. In a rough topographic area in Colorado, we computed a set of gravimetric geoid models based on different combination modes of satellite gravity models, terrestrial and airborne gravity data using the spectral combination method. The gravimetric geoid model obtained from the combination of satellite gravity model GOCO06S and terrestrial gravity data agrees with the GPS leveling measured geoid heights at 194 benchmarks in 5.8 cm in terms of the standard deviation of discrepancies, and the standard deviation reduces to 5.3 cm after including the GRAV-D airborne gravity data collected at ~ 6.2 km altitude into the data combination. The contributions of airborne gravity data to the signal and accuracy improvements of the geoid models were quantified for different spatial distribution and density of terrestrial gravity data. The results demonstrate that, although the airborne gravity survey was flown at a high altitude, the additions of airborne gravity data improved the accuracies of geoid models by 13.4%–19.8% in the mountainous area (elevations > 2000 m) and 12.7%–21% (elevations < 2000 m) in the moderate area in the cases of terrestrial gravity data spacings are larger than 15 km.


Introduction
In 2017, The Joint Working Group (JWG) 0.1.2 (Strategy for the Realization of the International Height Reference System (IHRS)) and JWG 2.2.2 (the 1 cm geoid experiment) of the International Association of Geodesy (IAG) jointly launched the Colorado geoid experiment. The goal of this experiment is to assess the repeatability of gravity potential values as IHRS coordinates using different geoid modeling methods, and to compare and evaluate the corresponding gravimetric geoid models. In this frame, the National Geodetic Survey (NGS) of the United States of America (USA) provided the geodesy community with terrestrial, airborne gravity and GPS (global positioning system) leveling data as well as digital elevation model (DEM) for a mountainous area of about 400 thousand km 2 in Colorado, which allowed the comparison of different methods and softwares for geoid computation using the same input dataset in this challenging area.
The experiment area is covered by both terrestrial and airborne gravity data, which are not only different in spatial distribution and spectral contents, but also in their error characteristics. Nowadays, satellite gravity models are routinely used to provide accurate long wavelength gravity field information for regional geoid modeling, while terrestrial and airborne gravity data contribute the medium and short wavelengths of gravity field. The major challenge for geoid modeling in this area is the proper combination of satellite gravity model, terrestrial and airborne gravity data.
The commonly used approaches for the combination of Earth gravity model, terrestrial and airborne gravity data can be divided into three categories. In the first approach, the airborne gravity data are firstly downward continued from the flight altitude to the ground or geoid to merge with the terrestrial gravity data, then the merged gravity data are used to compute the geoid by Stokes' integral (Novak and Heck 2002;Bayoud and Sideris 2003;Hwang et al. 2007;Forsberg et al. 2000Forsberg et al. , 2012Jekeli et al. 2013) or least squares collocation (Hwang et al. 2007;Scheinert et al. 2008;Shih et al. 2015). However, this approach involves two-step processing of airborne gravity data, introducing the edge effects twice. Additionally, the combination of different datasets is not weighted according to their spectral contents. The second approach is the least squares collocation (LSC) which combines all the input data and compute the geoid in one step (Hwang et al. 2007;Forsberg and Olesen 2010). The main advantage of LSC is that it can accommodate inhomogeneous data of different types and spatial resolutions. However, for application on massive datasets in large areas, the computation efforts involved are excessive to afford. The third approach is the spectral combination proposed by Wenzel (1982) and Sjöberg (1981), which can lead to a geoid solution with minimum least squares error if the spectral weights of each dataset are determined properly. The combination of terrestrial and airborne gravity data can be done in one step by means of spectral weights in surface integrals. It allows the integration of a great amount of gravity data in an efficient way on the contrary to LSC, and provides a flexible control on the data combination by spectral weights (Jiang and Wang 2016). In addition to these three approaches, radial base functions and wavelets can also be used for the combination of different types of gravity data, interested readers are referred to (Schmidt et al 2007;Klees et al 2008;Wittwer 2009;Panet et al 2010).
As a result, we used the spectral combination method for geoid determination in Colorado as our contribution to the geoid experiment. This paper summarizes the method, procedure, data, results and analyses of our Colorado geoid modeling experiment. Several gravimetric geoid models in Colorado were computed from different combination modes of satellite gravity models, terrestrial and airborne gravity data. The derived geoid models were then validated and compared using the high-precision GPS leveling measured geoid heights. Moreover, the contributions of airborne gravity data to geoid modeling were quantitatively evaluated for different data combination modes and terrestrial gravity data conditions. Gravimetric geoid modeling method is presented in Sect. 2. Section 3 describes the Colorado geoid experiment and the data used in this experiment. The geoid computations, results and analyses are explicated in detail in Sects. 4 and 5. Conclusions are given in Sect. 6.

Gravimetric geoid modeling method
The spectral combination of satellite gravity model, terrestrial and airborne gravity data is performed in Molodensky's theoretic frame. Molodensky's harmonic continuation method is employed to solve the geodetic boundary value problem (Heiskanen and Moritz 1967, p. 312;Hofmann-Wellenhof and Moritz 2005, p. 303-308;Wang et al. 2012); thus, the quasigeoid heights (height anomalies) are firstly computed and then need to be transformed into the geoid heights. Remove-compute-restore procedure with a high-degree Earth gravity model, e.g., the Earth Gravitational Model of 2008 (EGM08, Pavlis et al 2012Pavlis et al , 2013, is applied to account for the contribution outside local gravity data coverage. Residual terrain model (RTM) is used to represent the short wavelength components of gravity field generated by the high-frequency part of topography (Frosberg 1984).
For the spectral combination of satellite gravity model, terrestrial and airborne gravity data, the gravimetric geoid height can be decomposed into three components contributed from each dataset: where N is the gravimetric geoid height, ζ Sat , ζ Ter and ζ Air are the height anomaly contribution of the satellite gravity model, terrestrial and airborne gravity data, respectively, ζ 0 is the zero-degree term of height anomaly, and is the geoid-quasigeoid separation term.
Applying the classical remove-compute-restore procedure with a high-degree reference gravity model and representing the frequency gravity effects by the residual terrain model (Forsberg 1984), ζ Ter and ζ Air can be further decomposed and Eq. (1) is written as where ζ Res Ter is the residual height anomaly derived from terrestrial gravity data, and ζ Ref Ter is the reference height anomaly corresponding to terrestrial gravity. ζ Res Air and (1) Air are the corresponding terms of airborne gravity data, respectively. ζ RTM is the RTM effect of height anomaly, restoring the topographic effects which have been removed when computing the residual terrestrial gravity anomaly and the residual airborne gravity disturbance. ζ Sat can be computed from potential coefficients of the satellite gravity model using spectral weights W Sat (n) . ζ Res Ter is computed from residual terrestrial gravity anomalies using the Stokes' integral with spectral weights W Ter (n) . ζ Res Air can be computed directly from residual airborne gravity disturbances at flight altitude in one step using the generalized Hotine's integral with spectral weights W Air (n) . ζ the reference gravity model by spherical harmonic synthesis using spectral weights W Ter (n) and W Air (n) , respectively. ζ RTM is computed using the formula of rectangular prism for the RTM effects of height anomaly (Frosberg 1984;Nagy et al. 2000). Detailed formulas for the computation of above terms are referred to Eqs. (1-7) in Jiang and Wang (2016). Note that the gravity disturbance is used for airborne gravity data, because geodetic heights of the aircraft can be accurately known from the onboard GPS kinematic positioning.
The key problem for a decent geoid solution from the combination of satellite gravity model, terrestrial and airborne gravity data is to determine the proper spectral weights of each dataset, W Sat (n) , W Ter (n) and W Air (n) . The spectral weights of each dataset can be derived from the corresponding error degree variances using the condition of least squares residuals (Wenzel 1982;Sjöberg 1981). For Colorado geoid experiment, we selected the KTH (Royal Institute of Technology, Sweden) error degree variance estimation method, which applies the white and colored noise models to estimate the error degree variances of terrestrial and airborne gravity data (Sjöberg 1986;Ågren 2004). Formulas for spectral weight determination and KTH error degree variance estimation are referred to Eqs. (10-13 and 24-28) in Jiang and Wang (2016).
As the Colorado geoid experiment is standardized to be consistent with the IHRS definition, and geodetic reference system 1980 (GRS80, Moritz 2000) is adopted as the reference ellipsoid, and the zero-degree term of height anomaly is computed by where GM and GM GRS80 are the geocentric gravitational constants adopted by the IHRS and GRS80, W 0 is the IHRS reference gravity potential and U 0 is the normal gravity potential on GRS80 ellipsoid, r P is the geocentric radial distance of the computation point, and γ Q is the normal gravity of the corresponding point on the telluroid. The geoid--quasigeoid separation term is computed by (Flury and Rummel 2009) where g BO is the refined Bouguer gravity anomaly, H is the orthometric height, − γ is the mean normal gravity, and V TOP P and V TOP P 0 are the gravitational potentials of the topographic masses evaluated at the computation point and its projection point on the geoid. G is the constant of gravitation, ρ 0 is the average density of topographic masses, and g TC P is the terrain correction evaluated at the computation point. Equation (4) is an extension of the well-known approximation of the geoid-quasigeoid separation term in Heiskanen and Moritz (1967, Eqs. (8-102)) As a summary, the procedure of gravimetric geoid modeling by spectrally combining satellite gravity model, terrestrial and airborne gravity data is shown in Fig. 1.

The Colorado geoid experiment and the data
The reason for selecting Colorado as the experiment area was that it is the test area of the geoid slope validation survey 2017 (GSVS17) of NGS, where aboundant terrestrial, airborne gravity, GPS leveling and other terrestrial survey data are available. GSVS17 is the third survey after GSVS11 (Smith et al. 2013) and GSVS14 (Wang et al. 2017). The purpose of GSVSs is to evaluate the reachable accuracy of gravimetric geoid models and quantify the contribution of airborne gravity data of the 'Gravity for the Redefinition of the American Vertical Datum' (GRAV-D) project to the improvement of geoid models. While GSVS11 was performed over a low and flat topographic area in Texas and GSVS14 took place in Iowa in an area with moderate topography but significant gravity variation, GSVS17 selected a rough topographic area in Colorado which made it the most challenging case for geoid modeling among the three GSVS tests. Figure 2 shows the topography of the experiment area based on the shuttle radar topography mission (SRTM) DEM data (Farr et al. 2007) and the GSVS17 terrestrial survey traverse. The average, minimum and maximum topographic elevation of this area is 1733 m, 314 m and 4385 m, respectively. (4) Terrestrial gravity observations along with orthometric heights at 59,303 points in the area bounded by 35 • ≤ ϕ ≤ 40 • and 250 • ≤ ≤ 258 • were extracted from the gravity database of NGS. Over the area bounded by 34.5 • ≤ ϕ ≤ 38.8 • and 250.8 • ≤ ≤ 258.6 • , the GRAV-D airborne gravity data were collected using the Micro-g LaCoste TAGS (turn-key airborne gravimetry system). The aircraft flown in the nominal ground speed of 460 km/h at the average altitude of 6186 m (geodetic height). For the differential kinematic positioning of the aircraft, two ground static GPS stations were set up within the area. GPS data were collected and then processed based on the reference frame of GRS80 and ITRF2008 (international terrestrial reference frame 2008). The average precision of horizontal and vertical position of the aircraft is 5.2 cm and 8.7 cm (95% confidence interval), respectively. There are 49 west-east going data lines consisting of 283,716 observations and 7 north-south going cross lines, forming 269 crossover points. The designed spacing of adjacent data lines is 10 km and that of adjacent cross lines is 80 km. The total RMS (root mean square) error of the crossover discrepancies is estimated to be 2.2 mGal. An along track timedomain Gaussian filter of 120 s was applied to reduce the high-frequency noises in the airborne gravity data (GRAV-D Science Team 2017a, 2017b, 2017c).
To validate the gravimetric geoid models, NGS provided historic GPS leveling data at 194 benchmarks for the area bounded by 36 • ≤ ϕ ≤ 39 • and 251 • ≤ ≤ 257 • . The error budget of these GPS leveling measured geoid heights is estimated to be around 3 cm (Wang et al. 2020), which makes them the reliable independent data for the preliminary validation and analysis of the gravimetric geoid models. Note that the  (2009) is used, the Helmert orthometric heights should be corrected to the rigorous orthometric heights using consistent correction. From Hofmann-Wellenhof and Moritz (2005, Eqs. 4-31), Santos et al. (2006, Eq. 15) and Flury and Rummel (2009, Eq. 18-19), we obtained the correction ε H as where − g is the integral-mean value of gravity along the plumb line between the geoid and the Earth surface, − g H is the Helmert approximation of the mean gravity along the plumb line, H is the Helmert orthometric height, H O is the rigorous orthometric height, and g P is the gravity at the surface point. In the Colorado case, the gravity values at the 194 benchmarks were interpolated from the terrestrial gravity data.
The spatial distribution of terrestrial gravity points, airborne gravity survey lines, and historic GPS leveling benchmarks are plotted in Fig. 3. The statistics of the terrestrial gravity anomalies and the GRAV-D airborne gravity disturbances in this area are listed in Table 1.

General parameters and procedures
For the computation and analysis of gravimetric geoid models, the following general parameters and procedures were adopted: Average density of topographic masses ( ρ 0 ): 2670 kg m −3 . 2. Geographic limits of the computed geoid models: 35.5 • ≤ ϕ ≤ 39.5 • and 250.5 • ≤ ≤ 257.5 • . Grid resolution: 1′ × 1′. 3. Geoid computations were performed in the tidefree system. GRS80 was used as the reference ellipsoid, and the conventional constants are referred to Moritz (2000). 4. Atmospheric corrections were applied on the terrestrial and airborne gravity data using the formula (Dimitrijevich 1987, p. 4): where the orthometric height H in km. 5. After removing reference gravity values and RTM effects on gravity, the residual terrestrial gravity anomalies and residual airborne gravity disturbances were gridded into 1′ × 1′ grid using the program GEOGRID (Tscherning et al. 1991), separately. The maximum degree of the terrestrial gravity contribution was set to be 10,800 corresponding to the 1′ grid spacing of terrestrial data. 6. The radius of spherical cap was empirically chosen as 1 • for Stokes' and Hotine's integration. Geoid models based on the integration radius of 0.5 • , 1 • , 1.5 • and 2 • were tested against the GPS leveling data, among which the 1 • radius yielded the best agreement.
The topography of the experiment area in Colorado based on the SRTM DEM. The black line represents the GSVS17 terrestrial survey traverse 7. The RTM effects on gravity were computed using the 3′' × 3′' SRTM v4.1 DEM data and a mean DEM by the program TC (Forsberg 1984) with the integration radius of 100 km, and the resolution of mean DEM depends on the truncation degree of reference degree model. 8. For the KTH error degree variance estimation, the error degree variances of satellite gravity models were computed using their error coefficients. Considering that terrestrial gravity data errors are unknown and the RMS of the crossover discrepancies of airborne gravity data is 2.2 mGal (GRAV-D Science Team 2017c), we empirically assigned σ w = 3mGal, σ c = 1mGal for the terrestrial gravity anomalies and σ w = 1.5mGal, σ c = 0.5mGal for the airborne gravity disturbances, where σ w and σ c are the standard derivation (SD) of the white and colored noise, respectively. The other criterion for choosing σ w and σ c was that they yielded reasonable and realistic spectral weights for satellite gravity model, terrestrial and airborne gravity data (Figs. 4 and 5). The Nyquist frequencies of terrestrial and airborne gravity data were set to be N Ter Q = 10800 and N Air Q = 2000 , respectively. 9. Following the procedure in Fig. 1, ζ Sat was derived from the spectral weighted potential coefficients of satellite gravity models, ζ Ter was computed from the terrestrial gravity grid using the degree weighted Stokes' integral and ζ Air was computed directly from the airborne gravity grid at flight altitude using the degree weighted Hotine's integral. In this way, the satellite gravity model, terrestrial and airborne gravity data were spectrally combined to obtain the geoid models.

Combination of satellite gravity model and terrestrial gravity data
As the first part of the Colorado geoid modeling experiment, gravimetric geoid models are computed from the combination of satellite gravity model and terrestrial gravity data, which is the classical mode for data combination in areas where airborne gravity data are not available. Thus the terms for airborne gravity data in the   Figure 4 shows the spectral weights of terrestrial gravity data and each satellite gravity model, respectively. The contributions of terrestrial gravity data start at degree 126, 115, 99 and 125 for the combination with GOCO06S, ITSG-Grace2018s, TIM6 and DIR6. The spectral weights of terrestrial gravity data are equal to those of each satellite gravity model at degree 171, 149, 170 and 181. At these degrees, the satellite models and terrestrial gravity data have equal contribution in the data combination. After these degrees, the spectral weights of four satellite gravity models quickly decrease to zero at degree 231, 182, 231 and 253, while those of terrestrial gravity data increase to 1 symmetrically.
To compare the performance of the four satellite gravity models in the data combination, four gravimetric geoid models were computed by combining terrestrial gravity data with GOCO06S, ITSG-Grace2018s, TIM6 and DIR6, respectively. EGM2008 model up to degree and order 2190 was used as the reference gravity model. Statistics of the differences between the four geoid models and the GPS leveling measured geoid heights are shown in Table 2. It turns out that the biases and standard deviations of differences between the four geoid models and the GPS leveling measured geoid heights agree with each other in millimeter level, though the four satellite gravity models are derived from different sources of dataset. It seems that each satellite model is suitable for the data combination. Considering the multiple source of satellite Fig. 4 Spectral weights of satellite gravity models and terrestrial gravity data. The satellite gravity models are: a GOCO06S. b ITSG-Grace2018s. c TIM6. d DIR6 data and the time span and amount of observation data used in the modeling of GOCO06S, it was selected as the satellite gravity model to be combined with terrestrial and airborne gravity data for geoid modeling.

Choice of reference gravity model
Due to the limited coverage of terrestrial gravity data, reference gravity model is applied in a remove-compute-restore fashion to account for the contribution outside the terrestrial data domain. To select the proper reference gravity model and its truncation degree, three high-degree gravity field models up to degree and order 2190, EGM2008, EIGEN-6C4 (Förste et al. 2014) and XGM2019 (Zingerle et al. 2019) were compared. Table 3 shows the statistics of the differences between the gravimetric geoid models based on the three reference gravity models and the GPS leveling measured geoid heights; each reference model was truncated to degree 2190, 1080 and 720, respectively. The geoid models were computed by combining the GOCO06S model and terrestrial gravity data. The comparison results suggest that the three reference gravity models yield almost identical geoid model accuracy for the same truncation degree. The optimal degree of truncation for each reference gravity model is 2190, which results in the best accuracy of 5.8 cm in terms of the standard deviation of the differences. Lower truncation degree causes larger truncation errors of the gravity field, which will degrade the accuracy of geoid solution. The closeness of the results in Tables 2 and 3 demonstrates the reliability and stability of the spectral combination method for geoid computation; no matter which satellite gravity model or reference gravity model was used in the data combination, the results agree in millimeter level.
In consideration of the two facts: (1) EIGEN-6C4 used EGM2008-derived gravity anomalies over continents (Förste et al. 2014) and XGM2019 used a global 15′ × 15′ gravity anomaly data grid provided from the database of NGA (National Geospatial-Intelligence Agency) of the USA (Zingerle et al. 2019), which is nearly the same terrestrial dataset for the EGM2008 model.
(2) EGM2008 model adopted much less satellite observations compared with EIGEN-6C4 and XGM2019, and no GOCE data were used. It has the minimum data overlap with the satellite gravity model GOCO06S, which is an advantage for the explicit combination of GOCO06S with terrestrial and airborne gravity data in this case. We selected EGM2008 up to degree and order 2190 as the reference gravity model for geoid computation.

Combination of satellite gravity model, terrestrial and airborne gravity data
Because of the thriving of airborne gravity campaigns in many countries or regions, the combination of satellite gravity model, terrestrial and airborne gravity data are expected to be the main data combination mode for regional geoid determination. In this section, the geoid models are computed by spectrally combining GOCO06S model, terrestrial and GRAV-D airborne gravity data and then validated using the GPS leveling measured geoid heights.

Spectral weights of satellite gravity model, terrestrial and airborne gravity data
The spectral weights of satellite gravity model GOCO06S, terrestrial and GRAV-D airborne gravity data are plotted in Fig. 5. GOCO06S model takes nearly full weight in the data combination below degree 114, and after that the weight rapidly reduces to zero around degree 217. Starting from the degree of 139, the spectral weight of airborne gravity data quickly rises and reaches the maximum value of 0.63 until degree 204, then slowly reduces  to zero around degree 1355. The main contributions of airborne gravity data concentrate at the spectral band from degree 148 to 512 ( W Air (n) ≤ 0.2 ), and the airborne data weight more than the terrestrial data in the data combination below degree 308. After the degree of 1000, the contributions of terrestrial gravity data are dominant and those of airborne data can be neglected. Note that the along track time-domain Gaussian filter of 120 s could not remove all the airborne data noises. The characteristic of decreasing spectral weights of airborne gravity data at medium and high degrees is crucial for suppressing the high-frequency data noises, so that the one-step geoid computation directly from the airborne gravity data at flight altitude via degree weighted Hotine's integral can be stabilized.

Gravimetric geoid validation
The statistics of the differences between gravimetric geoid models based on different data combination modes and GPS leveling measured geoid heights are shown in Table 4, and the geoid height differences are plotted in Fig. 6. The differences for the geoid heights computed from EGM2008, EIGEN-6C4 and XGM2019 are also included for comparison. EGM2008 performs better than EIGEN-6C4 and XGM2019 in the experiment area. In comparison with EGM2008 model, the accuracy of the geoid model derived from the combination of GOCO06S model and terrestrial gravity data is slightly improved by 3 mm in terms of the standard deviation, and the two models show a 1.7 cm discrepancy in the bias. Since both geoid models are based on the exactly same terrestrial gravity dataset in this area, the differences in their performance can be attributed to the better satellite gravity data used in the latter and the difference between the local and global modeling method for gravity field refinement. After the addition of GRAV-D airborne gravity data, the geoid model accuracy is improved from 5.8 cm to 5.3 cm in terms of the standard deviation. Considering the rough topography in this area and the error budget of the historic GPS leveling data, this accuracy level is promising. Figure 7 shows the gravimetric geoid model based on the combination of GOCO06S model, terrestrial and GRAV-D airborne gravity data. The 86.3 cm bias between this gravimetric geoid model and the GPS leveling measured geoid heights is caused by the different potential values (W 0 ) adopted by the IHRS and the NAVD88.

Geoid-quasigeoid separation
For geoid determination at centimeter accuracy level in the Colorado experiment area with high, rugged and geologically complex topography, the transformation from height anomalies to geoid heights needs to be dealt with carefully. We used the rigorous formula (Eq. 4) derived by Flury and Rummel (2009) to compute the  geoid-quasigeoid separation term. This approach rigorously accounts for the contribution of the attraction of topographic masses to the mean gravity along the plumb line; it differs from the approximation (Eq. 6) of Heiskanen and Moritz (1967) in the term of V TOP  Table 5, and V TOP P 0 − V TOP P / − γ at the computation grids are shown in Fig. 8. Two geoid models were computed from the combination of GOCO06S model, terrestrial and airborne gravity data, one with the rigorous geoid-quasigeoid separation FR and the other with the approximation HM . Statistics of the differences between the geoid models and the GPS leveling measured geoid heights are given in Table 6. Comparison results show that the accuracy of geoid model based on FR is 8.6% better than that of geoid model based on HM , demonstrating the necessity of applying the rigorous modeling of geoid-quasigeoid separation rather than using the approximation in this mountainous area.

Contribution of airborne gravity data
To quantify the contribution of airborne gravity data to geoid modeling in the combination with terrestrial gravity data of different spatial distribution and density, terrestrial data points were resampled at the spacing of 5 km, 10 km, 15 km, 20 km, 25 km, 30 km, 35 km and 40 km on the basis of the original data. Two groups of gravimetric geoid models were computed; group A consists of nine geoid models based on the combination of GOCO06S model, airborne and terrestrial gravity data with original spacing and the spacing of 5 km, 10 km, 15 km, 20 km, 25 km, 30 km, 35 km and 40 km, respectively, while group B includes the nine counterparts derived from the combination of the GOCO06S model and these terrestrial data.
The changes of geoid heights contributed from airborne gravity data are shown in Fig. 9, which are the differences between the corresponding geoid models with and without the airborne data (group A-group B), and the statistics of the geoid height changes are listed in Table 7. In the cases of terrestrial gravity, the point spacings are less than or equal to 15 km, and the magnitude and characteristics of the geoid height changes are very close with each other, exhibiting short wavelength features. In the cases of terrestrial point spacing of 20 km, 25 km, 30 km, 35 km and 40 km, the magnitudes of geoid height differences increase with the widening of the terrestrial data spacing, and the medium wavelength dominant characteristics of the changes in geoid heights behave quite differently from one to another.
To evaluate the contribution of airborne gravity data to the improvement of geoid model accuracy in different topographic areas, we divided the 194 GPS leveling benchmarks into two groups, with 90 benchmarks at elevations > 2000 m (mountainous area) and 104 benchmarks at elevations < 2000 m (moderate area). Table 8 presents the statistics of the differences between the gravimetric geoid models and the GPS leveling measured geoid heights at elevations > 2000 m and elevations < 2000 m, respectively. With the increasing of the spacing between terrestrial gravity points, the number of resampled terrestrial data points decreases rapidly. For the geoid models in group B, the standard deviations of the differences rise quickly with the terrestrial gravity point spacing above 10 km in the mountainous area and 20 km in the moderate area, which is not the case if terrestrial gravity point spacings are less than or equal to 10 km in the mountainous area and 20 km in the moderate area, respectively.
Comparing the performance of gravimetric geoid models in group A and group B (Table 8), the accuracy improvements contributed by airborne gravity data are relevant to the spacing of the used terrestrial gravity data. In the cases of terrestrial data with original spacing and the spacing of 5 km, 10 km and 15 km, the inclusions of airborne gravity data improved the accuracies of the geoid models by 6.3%-11.9% in the mountainous area and 9.4%-12.3% in the moderate area. In the cases of terrestrial data with the spacing of 20 km, 25 km, 30 km, 35 km and 40 km, the additions of airborne data improved the geoid models in the accuracy by 13.4%-19.8% in the mountainous area and 12.7%-21% in the moderate area. The results demonstrate that: (1) airborne gravity data can only slightly improve the accuracy of geoid model if the used terrestrial gravity data are densely distributed with the spacing less than or equal to 15 km; (2) airborne gravity data are capable of effectively filling the data gaps of terrestrial gravity and obviously improving the geoid model accuracy when combined with sparsely distributed terrestrial data with the spacing larger than 15 km.
Overall, the improvement rates in the geoid model accuracies contributed by airborne data in the mountainous area are not superior to those in the moderate area. On the contrary, slight better improvement rates in the geoid model accuracies are observed in the moderate area than in the mountainous area. This contradicts to our expectation that airborne gravity data will lead to greater accuracy improvement for the geoid models in the mountainous area than those in the moderate area. The possible reason may be that, in the Colorado experiment area, the spatial distribution of terrestrial gravity data in the mountainous area is as dense as that in the moderate area, and even denser in some particular areas (Fig. 3), which is not the usual case for mountainous areas.
For each selection of terrestrial gravity point spacing in Table 8, the combinations of GOCO06S model, terrestrial and airborne gravity data lead to the gravimetric geoid models (group A) with the accuracy ranging from 4.7 cm to 5.2 cm in terms of the standard deviation in the moderate area. In the mountainous areas, the gravimetric  . 9 Changes of geoid heights contributed from airborne gravity data when combined with terrestrial gravity data of different spatial distribution and density. Terrestrial gravity data point spacings are: a original spacing. b 5 km. c 10 km. d 15 km. e 20 km. f 25 km. g 30 km. h 35 km. i 40 km geoid models of group A agree with the GPS leveling measured geoid heights in the standard deviation of the differences better than 7.2 cm, except for the case of terrestrial gravity data spaced at 35 km. This level of geoid model accuracy demonstrates the good quality of the GRAV-D airborne gravity observations in this area and the correctness and reliability of the spectral combination based geoid modeling method and procedure.

Conclusions
We have described the computations of 1′ × 1′ gravimetric geoid models as our contribution to the IAG Colorado geoid experiment. A series of gravimetric geoid models in the varied topography area of Colorado were derived based on the different data combination modes of satellite gravity model, terrestrial and GRAV-D airborne gravity data using the spectral combination method, and then validated against the historic GPS leveling measured geoid heights at 194 benchmarks provided by the NGS. The gravimetric geoid model, based on the combination of GOCO06S model and terrestrial gravity data, agrees with the GPS leveling measured geoid heights in 5.8 cm in terms of the standard deviation of the discrepancies, and this agreement reduces to 5.3 cm after the inclusion of airborne gravity data into the combination.
Based on the comparisons and analyses using the GPS leveling data, the accuracies of geoid solutions based on four satellite gravity models and three high-degree reference gravity models truncated to the same degree agree with each other in millimeter level. Additionally, the rigorous modeling of geoid-quasigeoid separation (Flury and Rummel 2009) improved the geoid model accuracy by 8.6% on the basis of the traditional approximation formula (Heiskanen and Moritz 1967), confirming the necessity of using the rigorous formula