Centimeter Precision Geoid Model for Jeddah Region (Saudi Arabia)

: In 2014, the Jeddah Municipality made a call for an estimate of a centimetric precision geoid model to be used for engineering and surveying applications, because the regional geoid model available at that time did not reach a su ﬃ cient precision. A project was set up to this end and dedicated sets of gravity and Global Positioning System (GPS) / levelling data were acquired in the framework of this project. In this paper, a thorough analysis of these newly acquired data and of the last available Global Gravity Field Models (GGMs) has been done in order to obtain a geoid undulation estimate with the prescribed precision. In the framework of the Remove–Compute–Restore (RCR) approach, the collocation method was used to obtain the height anomaly estimation that was then converted to geoid undulation. The remove and restore steps of the RCR approach were based on GGMs, derived from the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) and Gravity Recovery and Climate Experiment (GRACE) dedicated gravity satellite missions, which were used to improve the long wavelength components of the Earth’s gravity ﬁeld. Furthermore, two di ﬀ erent quasi-geoid collocation estimates were computed, based on gravity data only and on gravity plus GPS / levelling data (the so-called hybrid estimate). The best solutions were obtained with the hybrid geoid estimate. This was tested by comparison with an independent set of GPS / levelling geoid undulations that were not included in the computed solutions. By these tests, the precision of the hybrid geoid is estimated to be 3.7 cm. This precision proved to be better, by a factor of two, than the corresponding one estimated from the pure gravimetric geoid. This project has been also useful to verify the importance and reliability of GGMs developed from the last satellite gravity missions (GOCE and GRACE) that have signiﬁcantly improved our knowledge of the long wavelength components of the Earth’s gravity ﬁeld, especially in areas with poor coverage of terrestrial gravity data. In fact, the geoid models based on satellite-only GGMs proved to have a better performance, despite the lower spatial resolution with respect to high-resolution models (i.e., Earth Gravitational Model 2008 (EGM2008)).


Introduction
Since the available regional model for Saudi Arabia did not achieve the required centimetric precision [1], in 2014 the Municipality of Jeddah (Saudi Arabia) set up a project for the estimation of a geoid undulation model that is appropriate for promoting the use of Global Navigation Satellite System (GNSS) techniques in positioning and surveying applications. This project was managed by the authors, who estimated a centimetric geoid model for the Jeddah Municipality. This estimate has been revised, based on a deeper investigation of the new available Global Gravity Models (GGMs).
In recent years, the geodetic community has had the opportunity of reprocessing the entire data span of the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) mission [2] and combining these new satellite data with other types of data, like the satellite gravity data provided by the Gravity Recovery and Climate Experiment (GRACE, [3,4]), the Satellite Laser Ranging information from different satellites and also terrestrial gravity observations. Thus, a new generation of GGMs is now available and is routinely accessible through the International Center for Global Earth Models (ICGEM) website [5].
In light of these considerations, as mentioned before, the geoid estimation in the Jeddah Region was revisited in order to see if an improved solution could be computed based on these last new GGMs. Thus, the entire geoid estimate procedure has been repeated using the existing data and the most recent releases of the GOCE-based high-resolution EIGEN-6C4, GOCO05c and XGM2016 models ( [2,[6][7][8], respectively).
The procedure adopted for this estimate is Remove-Compute-Restore with the least squares collocation method and the geoid undulation estimate was based on gravimetry data, as well as on and Global Positioning System (GPS) and levelling measurements collected in this region in the framework of the aforementioned project.
Two solutions were computed, one based on gravity data only and one based on the integration of gravity and GPS/leveling data. In this second estimate, gravity data have been integrated with geoid undulation values derived by GPS and leveling measures on the same benchmarks to obtain a hybrid model, according to the approach in ( [9,10]).
These new local geoid undulation models computed for the Jeddah Region improved the previous estimate, thus matching the required precision. This result was also achieved due to the use of the satellite-only global geopotential models that proved to be effective in recovering/modeling the long wavelength of the Earth's gravity field.
In this paper, this revised geoid undulation computation procedure is described. In Section 2, the basics of the adopted methodology are described, while in Section 3 the GPS/levelling and the gravity data are documented. In Section 3, GGMs that were used in the computation are also described, while the Digital Terrain Model (DTM) used for evaluating the gravity effect of topography is also illustrated. Section 4 is devoted to the geoid estimate: the pure gravimetric geoid undulation estimate and the hybrid solution are described, together with their validation.
Finally, in Section 5, a discussion on the obtained results is given.

Methodology
The computation of the gravimetric geoid for Jeddah Municipality has been performed using the new terrestrial gravity data collected in 2014 and the Remove-Compute-Restore (RCR) procedure ( [11,12]), where the compute step, leading to the residual geoid from residual gravity data, has been accomplished by applying the least squares collocation solution operator.
Since the 90s, the RCR has been a well-tested technique to compute local geoid/quasigeoid undulation from gravimetric observations in the context of the Molodensky theory [13][14][15] and/or the Stokes theory. In this approach, the free-air gravity anomalies are modelled and/or processed according to their wavelengths. The long-wavelength components of gravity and height anomalies are modelled using a GGM, while the high-frequency components are evaluated by using a detailed DTM. In the so-called remove step, the long-wavelength components of the Earth's gravity field are removed from the observed free-air gravity anomalies (∆g FA ) using the GGM signal (∆g GGM ). In this long-wavelength component, the low-frequency gravity signal due to the topographic masses is also taken into account. The higher frequency of the topographic effect is then modelled and removed from Remote Sens. 2020, 12, 2066 3 of 18 the data by computing the Residual Terrain Correction (∆g RTC ) ( [11,12]). At the end of the remove step, the so-called residual gravity anomalies ∆g res are obtained according to Equation (1) In the compute step, the residual gravity ∆g res is then used to get the residual height anomalies (ζ res ) by applying, e.g., least squares collocation [13].
The long-wavelength components (ζ GGM ) coming from the GGM and the topographic high-frequency components (ζ RTC ) of the height anomaly are then computed and added to the residual height anomaly (ζ), to obtain the gravimetric model of the height anomaly ζ in the restore phase.
Finally, the geoid undulation N can be obtained from the height anomaly ζ using the following relationship [16]: where ζ is the height anomaly and N the geoid undulation expressed in meters, ∆g B is the Bouguer anomaly in gal and H is the orthometric height in km. This is an approximate formula, the validity of which has been discussed in [17]. However, given that the area under investigation is quite flat, formula (2) should be accurate enough for the purpose. The procedure for the gravimetric geoid estimate described previously is summarized in the scheme of Table 1. Table 1. Scheme of the Remove-Compute-Restore (RCR) procedure for the estimate of the geoid undulation, using gravity anomalies ∆g and the Molodensky theory. Free air (FA); derived by a Global Gravity Field Model (GGM); residual terrain correction (RTC); residual component (res). N is the geoid undulation; ζ is the height anomaly.

Remove step:
∆g FA − ∆g GGM − ∆g RTC = ∆g res Compute step: ζ res = Collocation(∆g res ) Restore step: The project that was developed in the Jeddah Region also provided a new set of points where the ellipsoidal and the orthometric heights were measured by GPS technique and spirit levelling, respectively. Thus, a set of about 435 geoid undulations is available from these points.
These so called GPS/levelling geoid undulations N are completely independent from the values obtained by the gravimetric geoid modeling. This independent database of geoid undulations allowed us to quantify the precision of the gravimetric estimation. These undulations were then also used in the computation of a hybrid geoid undulation. This hybrid geoid has, again, been obtained by applying the RCR procedure in which both the gravimetric data and the GPS/levelling undulations are combined using the least squares collocation method.
As described above, the data used in a geoid estimation procedure are the terrestrial gravity data, GGMs, the DTM for the computation of the residual terrain correction and GPS/levelling data. In the following sections, all steps contributing to the estimation of the pure gravimetric and the hybrid geoid are discussed in detail.

GPS/Leveling Data
The knowledge of ellipsoidal h and orthometric H heights of the same benchmarks is important for evaluating the precision of the gravimetric geoid undulation, because the difference between the ellipsoidal height and the orthometric height gives an independent set of geoid undulations, Remote Sens. 2020, 12, 2066 4 of 18 N GPS/lev . Moreover, this set of data can be used in the so-called hybrid geoid undulation estimate, where the gravity and the GPS/levelling data are used together to give a combined estimate of the geoid undulation.
For this reason, a network of 435 benchmarks were monumented in the Jeddah Region, following the main and the secondary roads of the Region, and GPS and spirit levelling measurements were performed. The spatial distribution of the data can be seen in Figure 5.
Spirit levelling surveys in closed loops were performed using digital levels and invar rods. The least squares adjusted orthometric heights have a precision that is less than 2 mm.
Concerning the GPS measurements to obtain the ellipsoidal heights h, static GPS measurements were performed for all the points of the network, using dual-frequency GPS receivers. This survey required 165 measurement sessions. Each point was measured twice, using different GPS receivers to minimize the random errors. The least squares adjustment of the network was done using Trimble Business Center software with a minimum constraint fixing the Jeddah Continuously Operating Reference Stations (CORS). The precision of the GPS-derived ellipsoidal heights was less than 3 cm. Thus, by error propagation, we can assume that the precision of the GPS/levelling geoid undulation is around 3 cm.

The Gravity Data Set
In order get a reliable geoid estimate in the Jeddah region, dedicated gravimetric campaigns were performed and the Micro-g LaCoste Company collected the measurements. A set of 20 values of absolute gravimetry points, plus a set of 1648 relative values, have been measured using an A10X gravimeter by Micro-g La Coste Inc and an Autograv gravimeter, respectively. The distribution of the gravimetry points was in the range of 2-5 km over the entire Jeddah region.
All absolute gravity measurements were taken on pre-located and installed concrete pillars throughout the municipality of Jeddah. For the sites located in populated areas, the pillar location was typically installed inside government-controlled facilities in order to take advantage of presumed site longevity. The pillars were large concrete blocks with a (roughly) 60 cm × 60 cm surface area and all exhibited excellent stability with respect to instrument noise.
All the gravity absolute measures were within the expected instrumental repeat values, that is the repeat standard deviation was under 5 µGal at quiet sites.
The relative gravity measurements of 1648 points were performed using Autograv CG5 gravimeter and the coordinates of these gravimetric benchmarks were determined using STOP&GO GPS observation techniques. The precision of these relative gravimetric acquired data is consistent with the precision of the absolute gravity measures.
The measured gravimetry values g obs have been reduced to free-air gravity anomalies, computing the normal gravity γ according to the Geodetic Reference System 1980 (GRS80) [18] using the following formula ∆g FA (P) = g obs (P) -γ(Q) where P represents a point on the actual topographic surface and Q is the corresponding point that has an ellipsoidal height h(Q) equal to the orthometric height H(P).
The observed gravity values have been checked for outliers using the procedure described in [19]. This procedure requires the use of reduced free-air gravity anomalies. To this end, the GGM Earth Gravitational Model 2008 (EGM2008) to a degree and order (d/o) of 2190 has been used [20]. Based on the correlation length of the empirical covariance function of these residuals (∆g res ), which is about 7.7 km (0.07 • ), subsets of the gravity data were selected and the related statistics were computed, block by block, as the standard deviation σ. The measured gravity values that were out of the 2σ interval centered in the block mean have been marked as possible outliers (∆g res *). To confirm that these selected values were outliers, Student's t-inference test at a significance level of 5% was performed. The hypothesis H 0 : ∆ĝ res = ∆g * res was tested, where ∆ĝ res represents the predicted value obtained using a least squares collocation interpolator based on the data included in a 0.07 • cap centered in each Remote Sens. 2020, 12, 2066 5 of 18 investigated gravity point. If the test fails, the selected gravity value is confirmed as an outlier. In the collocation procedure, the following covariance function (Eq. 4) has been used: where ψ is the spherical distance between gravimetric points, C 0 is the variance of the data included in the current cap and α is the correlation length.
Following this procedure, 55 gravity values (3.3%) were selected as possible outliers, but only 21 of them were confirmed as such (1.2% of the total data set, with maximum and minimum difference values of 110.796 mgal and −32.590 mgal, respectively). In Figure 1, free-air gravity anomaly values are plotted in a color scale. The black stars are the suspected outliers pointed out by the 2σ rule, whereas the red ones are the outliers confirmed by Student's t-inference test. The new terrestrial gravity dataset was then combined by the sea gravity data provided by Bureau Gravimétrique International (BGI) [21] to improve the spatial coverage of the data and reduce the edge effects, and thus the gravity free-air database was formed with 4411 points.

The Global Geopotential Model
Nowadays, different Global Geopotential Models (GGM) are available thanks to the new satellite gravity missions, especially the Gravity Recovery and Climate Experiment (GRACE) [3,4] and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) [2]. Different GGMs were taken into account and tested, looking for the one giving the best fit to the new terrestrial gravity dataset used in this study. Among the available GGMs investigated, satellite-only GGMs have also been considered, even if their maximum degree and order (300) in the spherical harmonic expansion is significantly lower than the one of combined GGMs (satellite and terrestrial combined solutions) that are complete to a higher degree and order (d/o) (e.g., 2190).
This has been done in order to check how effective new satellite-only GGMs, based on the GRACE and GOCE satellite missions, are in representing the low-frequency components of the gravity field, while also comparing them to the high-degree GGMs that are based on the combination of satellite and ground data. What we wanted to investigate with this comparison is the effectiveness of the high-degree GGMs in areas where only a small amount of terrestrial data were available, as is the case for the Arabic Peninsula. As is well known, in the EGM2008 computation [20], the 5 arc minute area mean gravity anomalies used for this geographic area are labelled as "fill-in" data. This means that, for this region, the 5 arc minute gravity anomalies were derived by 15 arc minute area mean value classified data. Thus, in this region, the quality of the data is quite poor. Furthermore, these "fill-in" data are also present in the other high-degree GGMs since they have been estimated using the same terrestrial database of EGM2008.
In this paper, GGMs labelled as timewise (GO_TIM), direct (GO_DIR) and spacewise (GO_SPW) are the satellite-only GGMs that have been considered [2]. The R4, R5 and R6 labels indicate the release of the GOCE-based GGMs. For instance, the R5 and R6 releases use GOCE data, spanning the complete mission lifetime, including the low-orbit data collected in the re-entry phase of the satellite in November 2013. The three solutions (GO_TIM, GO_DIR and GO_SPW) differ for the processing strategies of the same GOCE observations. GO_TIM and GO_DIR are based on the so-called timewise (GO_TIM) and direct (GO_DIR) approaches, respectively, whereas the latter is based on the so-called spacewise (GO_SPW) method [2]. The R6 release in each processing strategy (direct, timewise and spacewise) differs from the R5 release for an improved re-calibration of the GOCE satellite gradiometer components [22].
In Table 2, the global potential models are tested and their main characteristics are reported. The models highlighted in a grey color are the ones already tested and used in the 2014 geoid undulation estimate, whereas the others are the most recent models used in the new validation procedure for the geoid estimate. Table 2. Maximum degree/order (d/o) harmonic coefficients of the Global Gravity Models (GGMs) used in the computation and tested in relations to the spherical window size of the moving average of the detailed Digital Terrain Model (DTM). For each GGM, the fourth and fifth columns represent the chosen combination parameters (d/o and window size), giving the best statistics for the gravity residuals. In the last column, S means satellite, A altimetry and T gravity terrestrial data. The spherical models highlighted in grey are GGMs already tested in 2014 when the project was set up. Type of Data EGM2008 [20] 2190 360, 720, 1080, 1440, 1800, 2190 1800 5' S-A-T EIGEN-6c2 [5] 1949 360, 720, 1080, 1440, 1800, 2190 1800 5' S-A-T EIGEN-6c4 [6] 2190 360, 720, 1080, 1440, 1800, 2190 1800 5' S-A-T GOCO05c [7] 720 360, 720 360 35' S-A-T GO_SPW_R4 [5] 280

Digital Terrain Model
The availability of a detailed Digital Terrain Model (DTM) is a key point in the computation of the residual terrain correction. Although a high resolution DTM (5 m × 5 m) is available for the Jeddah Municipality, the 3" × 3" DTM from the Shuttle Radar Topography Mission (SRTM3 [28] was preferred, because the gravity measurements are affected by the signal of the topography in an area larger than the computation area. Thus, a digital terrain model in an area wider than the one containing the gravity data is needed. For this reason, a homogeneous terrain model was preferred with respect to the highest resolution Jeddah DTM, which was available only over the Jeddah Municipality region. As the computation zone includes a portion of the Red Sea, the DTM model must be integrated with a bathymetric model. To this end, the bathymetric General Bathymetric Chart of the Oceans (GEBCO, [29]) grid has been merged with the SRTM3 model. The GEBCO grid has been interpolated by means of a bilinear interpolation to have the same resolution as the SRTM3 grid. Furthermore, the GEBCO grid, with a resolution of 30" × 30", has also data on land. These values have been used to fill the gaps in the SRTM3 grid. Figure 2 illustrates the assembled detailed DTM.
As is well known, Residual Terrain Correction (RTC) consists of the computation of the gravimetric effect of the terrain bounded by two surfaces. We computed this effect using the TC software of the GRAVSOFT package [11], selecting the standard density values 2670 kg m −3 and 1030 kg m −3 for the land and for the water, respectively. The first surface is the topographic surface, represented by the detailed DTM in Figure 2. The reference DTM is smoother than the first ones and should be correlated to the topography effect that is already taken into account by the used GGM (see Table 2). This surface is not known a priori and it depends on the maximum degree of the spherical harmonic expansion of the applied GGM. This smoother surface can be statistically identified by testing different smoothed versions of the available detailed DTM, computed by a moving average procedure based on different window sizes (in this work, from 5' to 50'). The statistics of the residual gravity anomalies (∆g res ) are then used in the choice of the best window size for the DTM averaging procedure. The smoothed reference DTM that gives the lowest values of the mean and standard deviation of the residual gravity anomalies is the proper reference DTM to be considered in the RTC computation. As the average DTM is correlated to the geopotential model used to compute the residuals, for each considered GGM, an analysis of the size of the averaging window has been performed. By using this  As is well known, Residual Terrain Correction (RTC) consists of the computation of the gravimetric effect of the terrain bounded by two surfaces. We computed this effect using the TC software of the GRAVSOFT package [11], selecting the standard density values 2670 kg m −3 and 1030 kg m −3 for the land and for the water, respectively. The first surface is the topographic surface, represented by the detailed DTM in Figure 2. The reference DTM is smoother than the first ones and should be correlated to the topography effect that is already taken into account by the used GGM (see Table 2). This surface is not known a priori and it depends on the maximum degree of the spherical harmonic expansion of the applied GGM. This smoother surface can be statistically identified by testing different smoothed versions of the available detailed DTM, computed by a moving average procedure based on different window sizes (in this work, from 5' to 50'). The statistics of the residual gravity anomalies (Δgres) are then used in the choice of the best window size for the DTM averaging procedure. The smoothed reference DTM that gives the lowest values of the mean and standard deviation of the residual gravity anomalies is the proper reference DTM to be considered in the RTC computation. As the average DTM is correlated to the geopotential model used to compute the residuals, for each considered GGM, an analysis of the size of the averaging window has been performed. By using this approach, we optimized both the d/o of each considered GGM and the reference DTM to be used in the RCR procedure (Table 2).

The Remove-Compute-Restore Results
In the remove step of the RCR procedure, gravity data are reduced for the low-and high-frequency components. In Table 3, the statistics of the remove step are given for each used GGM. For each solution based on different GGMs, the statistics of the ∆g res show the different features of the gravity residuals in relation to the d/o of the used GGM. This can also be seen in the structure of the empirical auto-covariance functions that are described herein.
The prediction on a 1 × 1 mesh grid of the residual height anomaly estimation ζ res was performed using the least squares collocation (LSC) method ( [13,14]). The LSC method is based on the estimate of the empirical auto-covariance function (ACF) of the residual gravity anomalies and its interpolation with a positive definite function, called the covariance function model ( [30][31][32]).
The magnitude of the variance of the signal and the correlation length of the ACFs are different, due to different GGMs and RTC computations used in the remove step. In Figures 3 and 4, empirical auto-covariance functions and best-fit covariance models for the reduced values obtained with EGM2008 and GO_DIR_R5 are plotted, respectively. As expected, the residual free-air gravity anomalies obtained with GO_DIR_R5 have a longer correlation length than the corresponding ones obtained using a higher order degree model, i.e., EGM2008, even though the difference in the correlation length is not as sharp.   Table 2). Auto-covariance function (ACF).  Table 2). Auto-covariance function (ACF).   28.0 18.9 −47.0 47.0 The Compute step of the RCR procedure was then followed by the restore step. The long-wavelength components ( GGM) coming from the different GGMs considered in this study and the topographic short-wavelength components ( RTC) of the height anomaly were then computed and added to the residual height anomaly ( ), obtained by LSC.  The Compute step of the RCR procedure was then followed by the restore step. The long-wavelength components (ζ GGM ) coming from the different GGMs considered in this study and the topographic short-wavelength components (ζ RTC ) of the height anomaly were then computed and added to the residual height anomaly (ζ), obtained by LSC.

Validation of the Gravimetric Geoid Model
The precision of the estimated gravimetric geoid models has been assessed by comparing their geoid undulation values with those obtained by GPS/levelling observations collected in the framework of this project for a set of points distributed in the investigated area (see Figure 5).

The Hybrid Estimate of the Height Anomaly and Geoid Undulation.
The availability of a dense set of GPS/levelling points-that is, points with known orthometric and ellipsoidal heights, allows for the estimation of a hybrid geoid model for the Jeddah Municipality ( [9,10]). The gravimetric observations and the height anomalies coming from the difference in the GPS and spirit levelling observations can be integrated together, again, by using the LSC estimator.
The integration procedure is based on a revised "Remove-Compute-Restore" approach: the residual gravity anomalies (Δgres) are integrated into the LSC procedure with the residual GPS/levelling height anomalies ( GPS/lev_res obtained by GPS/levelling geoid undulations and applying formula (2)), found by subtracting the long-wavelength component ( GGM), using a global geopotential model and the topographic short-wavelength component ( RTC). In the LSC estimator, the vector of the observations is now given by = (Δ , / _ ) and the covariance matrix C is a block matrix, composed of the gravity covariance matrix CΔgΔg, the covariance matrix Cζζ of the As we have adopted the Molodensky theory, the estimated height anomalies have to be transformed in geoid undulations using Equation (2) to be coherently compared with the geoid undulations obtained by the orthometric H and the GPS ellipsoidal heights h GPS , that is N GPS/lev = h GPS -H. In Table 5, statistics of the differences between the gravimetric geoid undulations and N GPS/lev are reported. Independently of the GGM used, the differences present significant discrepancies. This is mainly due to the different reference systems in which the gravimetric and the GPS/levelling-derived geoid undulations are given. Thus, in order to enable a proper comparison, a datum transformation must be estimated. According to Section 2 [33], a three-parameter transformation (dx, dy, dz) using formula (5) can be used: N grav = N GPS/lev + ∆N(θ,λ) = N GPS/lev + dx sinθ cosλ+ dy sinθ sinλ + dz cosθ (5) where θ is the co-latitude, that is θ = 90 • − φ, and λ φ are the longitude and the latitude of the computation point. The statistics of the least squares differences represent an estimate of the precision of the gravimetric geoid undulation (Table 6). An outlier rejection is also performed on the least squares differences of the datum transformation, under the hypothesis of normally distributed values at a significance level of 1%. As can be observed, the number of points reported in Table 5 is different from one solution to another since different geoid solutions were used and the outlier rejection gave us different results. Table 6. Statistics of the differences between the gravimetric geoid undulations (labelled with the name of the GGM used) and the GPS/levelling values after the datum transformation (cm).

Nres
Number As the standard deviation of the differences is about 7.5-8.0 cm, the geoid estimates obtained using the satellite-only GGMs and the XGM2016 GGM are those that are in better agreement with the GPS/levelling geoid undulation values with respect the model based on the other GGMs. These statistics confirm our hypothesis about the reliability of the satellite-only GGMs in areas poorly covered by terrestrial gravity data and show the important contribution of the GOCE mission.

The Hybrid Estimate of the Height Anomaly and Geoid Undulation
The availability of a dense set of GPS/levelling points-that is, points with known orthometric and ellipsoidal heights, allows for the estimation of a hybrid geoid model for the Jeddah Municipality ( [9,10]). The gravimetric observations and the height anomalies coming from the difference in the GPS and spirit levelling observations can be integrated together, again, by using the LSC estimator.
The integration procedure is based on a revised "Remove-Compute-Restore" approach: the residual gravity anomalies (∆g res ) are integrated into the LSC procedure with the residual GPS/levelling height anomalies (ζ GPS/lev_res obtained by GPS/levelling geoid undulations and applying formula (2)), found by subtracting the long-wavelength component (ζ GGM ), using a global geopotential model and the topographic short-wavelength component (ζ RTC ). In the LSC estimator, the vector of the observations is now given by y = (∆g res , ζ GPS/lev_res ) and the covariance matrix C is a block matrix, composed of the gravity covariance matrix C ∆g∆g , the covariance matrix C ζζ of the height anomalies and the cross-covariance matrix C ζ∆g . The covariance matrix of the gravity data is the same used in the pure gravimetric solution (Figure 1), whereas the covariance matrix C ζζ and the cross-covariance C ζ∆g are obtained via covariance propagation from C ∆g∆g . On the main diagonal of the C ζζ matrix, the height anomaly variance σ 2 n (ζ) is added. This value is chosen by assuming that the standard deviation of the height anomalies obtained from GPS/levelling data is about 3 cm, which represents a realistic value, taking into account the precision in the GPS and the levelling observations given in Section 3.
The result of the LSC computation is the set of hybrid residual height anomalies ζ H res , as shown in Equation (5), to which the long-wavelength term (ζ GGM ) and the topographic high-frequency component (ζ RTC ) will be added in order to complete the restore step.
As we have presented in the previous section, the available GPS/levelling data set consisted of 443 points. Some of these values have been detected as outliers in the validation procedure of the pure gravimetric geoid model, so these points have been removed from the data set considered in the hybrid estimate. In the hybrid geoid estimate, only the ζ GPS/lev_res values in the range from −10 to 10 cm, with respect to the pure gravimetric geoid undulations, have been taken into account. These values, about 80% of the total database, are selected according a 3σ criterion, based on the precision of the GPS-derived heights (see Section 3). Furthermore, this reduced data set was randomly divided into two subsets, A and B, where about the 2/3 of the entire dataset is in set A. Set A was integrated with the gravimetric observations to obtain the hybrid quasigeoid model, whereas set B was used as an independent dataset to validate the results.
At the end of this estimate, the height anomaly values were converted into geoid undulation, according to the approximate formula of Equation (2).

Validation of the Hybrid Geoid Model
The first validation of the data was performed on set A of the GPS/levelling geoid undulations. This is only a consistency control to check whether the relative weights between gravity and geoid undulation data are correctly set up. In Table 7, the statistics of the differences between the hybrid geoid models and GPS/levelling geoid undulations, after the datum shift reduction, are reported. The second independent validation was then done for set B. The hybrid geoid model was interpolated onto the points in set B with a bilinear interpolation and the differences between the GPS/levelling geoid undulations were computed (statistics in Table 7).
The integration procedure has significantly improved the precision of the geoid estimate in the Jeddah Municipality. The standard deviations of the best pure gravimetric solutions are of the order of 7.5-7.8 cm, while the standard deviation of the hybrid solution decreases to 3.7 cm. The different hybrid solutions obtained by using different satellite-only GGMs have similar statistics, better than the ones of the combined GGMs, such as EGM2008, EIGEN-6c2, EIGEN-6c4 and GOCO05c. Furthermore, the recent XGM2016 model gives results of the same order of precision of the solutions based on the satellite-only GGMs.
The comparison of different geoid undulation solutions shows good agreement in the area covered by terrestrial gravity data and GPS/levelling points, as shown in Figure 6, where the differences between the hybrid geoid model determined by the solution GO_SPW_R5 and the ones developed by the solution GO_DIR_R5 are reported. Outside the area covered by ground data, the differences are larger, in the range from −10 to 10 cm when we compare solutions based on satellite-only models. On the other hand, when comparing, for instance, the solution based on GO_SPW_R5 and the solution based on XGM2016, one has higher differences in the sea, even where gravity data are present (Figure 7). This behavior seems to prove that, between the two geoid models, there is a tilt effect that has been mitigated only in the central part of the computation area where the GPS/levelling data were collected.

Conclusions
In 2014, the Jeddah Municipality made a call for an estimate of a centimetric precision geoid model to be used for engineering and surveying applications.. A project was set up to this end and dedicated sets of gravity and GPS/levelling data were acquired in the framework of this project. In this paper, a thorough analysis of these newly acquired data, of the DTM and of the available GGMs has been exploited in order to obtain a geoid undulation estimate with the prescribed precision. In the framework of the RCR approach, the least squares collocation method was used to estimate the height anomaly estimation, which was then converted to geoid undulation. This procedure has been applied based on different GGMs with different d/o resolutions. Furthermore, two different quasi-geoid collocation estimates were computed, based on gravity data only and on gravity plus GPS/levelling data (the so-called hybrid estimate). The best solutions were obtained in the hybrid geoid estimate. This was tested by comparison with an independent set of GPS/levelling geoid undulations that were not included in the computed solutions. Through these tests, the precision of the hybrid geoid, i.e., 3.7 cm, proved to be better, by a factor of two, compared to the precision of the pure gravimetric geoid, i.e., 7.5 cm.
As a byproduct, the estimate of the geoid undulation for the Jeddah Municipality using the newly collected data also gave us the opportunity to test the newly developed GGMs, estimated with the data from the GRACE and the GOCE satellite gravity missions.
All the solutions based on satellite-only GGMs had higher precision in terms of the geoid undulation estimate. Moreover, the old releases of the satellite-only GGMs, as well as the R4 release

Conclusions
In 2014, the Jeddah Municipality made a call for an estimate of a centimetric precision geoid model to be used for engineering and surveying applications. A project was set up to this end and dedicated sets of gravity and GPS/levelling data were acquired in the framework of this project. In this paper, a thorough analysis of these newly acquired data, of the DTM and of the available GGMs has been exploited in order to obtain a geoid undulation estimate with the prescribed precision. In the framework of the RCR approach, the least squares collocation method was used to estimate the height anomaly estimation, which was then converted to geoid undulation. This procedure has been applied based on different GGMs with different d/o resolutions. Furthermore, two different quasi-geoid collocation estimates were computed, based on gravity data only and on gravity plus GPS/levelling data (the so-called hybrid estimate). The best solutions were obtained in the hybrid geoid estimate. This was tested by comparison with an independent set of GPS/levelling geoid undulations that were not included in the computed solutions. Through these tests, the precision of the hybrid geoid, i.e., 3.7 cm, proved to be better, by a factor of two, compared to the precision of the pure gravimetric geoid, i.e., 7.5 cm.
As a byproduct, the estimate of the geoid undulation for the Jeddah Municipality using the newly collected data also gave us the opportunity to test the newly developed GGMs, estimated with the data from the GRACE and the GOCE satellite gravity missions.
All the solutions based on satellite-only GGMs had higher precision in terms of the geoid undulation estimate. Moreover, the old releases of the satellite-only GGMs, as well as the R4 release of the GO_DIR approach, show a good performance, somehow even better than the newer solutions.
Thus, for this reason, the geoid undulation model computed in 2014 for the Jeddah Municipality can be still considered reliable.
Finally, it must be remarked that the new XGM2016 GGM is also a promising model in view of the computation of the new generation of combined GGMs. In fact this model can be considered the precursor step for the combined Earth Gravitational Model 2020 (EGM2020).
Through the definition of the Jeddah2014 Vertical Reference Frame, different research groups have improved and validated geoid undulation models all over the Kingdom of Saudi Arabia (KSA), i.e., KSA-GEOID17 and SAGE013 [34][35][36], so we believe that this work, although focused on a small region of the KSA, could give important information about the GGMs that can be used and about their expected precision.