EVALUATION OF RECENT GLOBAL GEOPOTENTIAL MODELS BY GNSS/LEVELLING DATA: INTERNAL AEGEAN REGION

Global geopotential models of spherical harmonic coefficients are used to determine the external gravitational field of the Earth. These coefficients are derived from satellite orbit perturbations, terrestrial gravity anomalies and altimeter data. Hundreds of thousands of coefficients and standard deviation values for these coefficients are estimated from millions of observation. Measurement amount, homogenous distribution of the measurements of global scale, different measurement types reflecting the different frequencies of the gravity signal and measuring-assessment techniques affect the model accuracy directly. Starting from 1960’s and lasts to the present day and also gaining new acceleration with the satellite gravity field missions, every outcome of the studies related to the determination of the global Geopotential model is experienced by a series of validation tests. Accuracy of the model can either be determined from the estimated error degree variances concerning the coefficients (interior validation) or comparison of geoid heights, gravity anomalies, gravity disturbances and components of vertical deflection calculated from the model with terrestrial measurements directly (outer validation). In this paper, recent global geopotential models are primarily explained. Global geopotential models are compared with GNSS/levelling data of the study area. The objective of this comparison is to determine the best fit global geopotential model which will contribute to the study of Turkish geoid determination.


INTRODUCTION
The geoid surface serves as a reference for most applications that require a datum for determining topographic heights or ocean depths. The improvements derived from recent satellite gravity missions have significantly improved earth gravity field knowledge, such that global geopotential models (GGMs) representing the Earth's gravity field have acquired greater importance to the geosciences.
The technological and scientific developments in satellite techniques and computation algorithms provide significant improvements in the determination of the global gravity field models. Since the launch of the CHAllenging Minisatellite Payload (CHAMP), Gravity Recovery And Climate Experiment (GRACE), and Gravity field and steady-state Ocean Circulation Explorer (GOCE) missions (2000, 2002, and 2009 respectively), numerous GGMs have become available to the scientific community through the public domain (http://icgem.gfzpotsdam.de/ICGEM). Especially, the releases of the Earth Gravitational Model 2008 (EGM2008) by the US National Geospatial Intelligence Agency (Pavlis et al., 2008) and European Improved Gravity model of the Earth by New techniques (EIGEN-5C) by the GFZ-GRGS cooperation (Förste et al., 2008) are significant achievements in the determination of the Earth's mean gravity field. These high-degree models lead to significant improvement of our knowledge of the long wavelength part of the Earth's static gravitational field, and thereby of the long wavelengths of the geoid. Therefore, corresponding improvements are expected for precise regional geoid model determination because regional geoid models typically include a GGM as underlying geopotential representation (Erol et al., 2009).
The geodesy community engaged in comprehensive efforts for the comparison and validation of GGMs using several techniques and independent data sets that were not used for the development and evaluation of these GGMs. To improve local geoid models, it is essential to select the best GGM for the studied area. In the selection of a GGM for geoid determination, published error estimates for GGMs are frequently not used to judge which GGM is best for a certain region. This is because the published quality estimates may be too optimistic and/or presented as global averages and thus not necessarily representative of the performance of the GGM in a particular region. Hence, the user of a GGM should perform his own accuracy and precision verifications (Kiamehr and Sjöberg, 2005).
Continuous developments in the acquisition, modelling and processing of GPS data have provided geodesists with highly reliable and precise external control to evaluate global and regional models for the Earth's gravity field (Kotsakis, 2008). The main objective of this study is comparing EGM2008, EIGEN-6C4, GOCE and EGM2008 COmbined model (GECO), The Combined Gravity Models (GGM05C and GOCO05C). Geoid heights determined from GNSS/Levelling over the Internal Aegean Region study area were used to quantify the GGMs' accuracy in order to find the geopotential model that best fits the study area for further geoid determination at regional and national scales.

GNSS/Levelling
Global Navigation Satellite System (GNSS)-derived ellipsoidal heights refer to a reference ellipsoid, while orthometric heights refer to an equipotential reference surface determined through levelling. When these heights are collocated at the same benchmark, their difference can be used to determine geoid height through a geometrical approach. The GNSS/Levelling geoid heights are computed by the following equation (Heiskanen and Moritz, 1967): where N is the geoid height, h is the ellipsoidal height computed from GNSS and H is the orthometric height computed from levelling ( Fig. 1). Geoid heights have been computed based on the known ellipsoidal and orthometric heights (Banarjee et al., 1999). Eq. (1) is not exact due to the ignorance of the deflection of the vertical (). Nevertheless, it is accurate enough for most practical applications, because  has a negligible influence (sub mm-order) on the orthometric height (Tenzer et al., 2005).

Global Geopotential Model
For a better determination of orbits and height systems in science and engineering, it is necessary to significantly improve our knowledge of the gravity field of the Earth, both in terms of accuracy and spatial resolution (Rummel et al., 2002). The GGM is used to determine the long wavelength part of the earth's gravity field and comprises a set of fully-normalized, spherical harmonic coefficients that are obtained from geopotential solutions (Mainville et al., 1992). These coefficients are determined from the incorporation of satellite observations, land and ship-track gravity data, marine gravity anomalies derived from satellite radar altimetry and airborne gravity data (Rapp, 1997).
The geoid height (N) can be represented by a set of spherical harmonic coefficients (in spherical approximation) with the following equation (Heiskanen and Moritz, 1967): where (θ, λ) co-latitude and longitude of the computation point, R is the mean radius of the Earth, m P is the associated Legendre polynomials, m C and m S  are the spherical harmonic coefficients for degree l and order m, respectively.

Study Area and Source Data
The area is located in the internal Aegean region of Turkey within the geographical boundaries: 37 0 .3083 N ≤ φ ≤ 40 0 .4417 N; 28 0 .4833 E ≤ λ ≤ 32 0 .7167 E defining a total area of  133000 km 2 (350 km x 380 km) with a rough topography (Fig.  2). All our GGM evaluation tests based on geoid height refer to the 87 points that belong to Turkish National Fundamental GPS Network (TNFGN) (Fig. 3).  Ellipsoidal heights at 87 points have been determined using dual-frequency GNSS receivers and antennas with respect to TNFGN (aligned to ITRF96) (reference epoch 2005.00) and orthometric heights at these points have been determined through spirit levelling with respect to the Turkish National Vertical Control Network (fixed to local mean sea level of the Antalya tide gauge). Geoid heights at 87 TNFGN points have been computed according to the Eq. (1) based on the known ellipsoidal and orthometric heights above.

Earth Gravitational Model 2008
EGM2008 is a spherical harmonic model of the Earth's gravitational potential, developed by a least squares combination of the ITG-GRACE03S gravitational model and its associated error covariance matrix, with the gravitational information obtained from a global set of area-mean free-air gravity anomalies defined on a 5 arc-minute equiangular grid. This grid was formed by merging terrestrial, altimetry-derived, and airborne gravity data. Over areas where only lower resolution gravity data were available, their spectral content was supplemented with gravitational information implied by the topography. EGM2008 is complete to degree and order 2159, and contains additional coefficients up to degree 2190 and order 2159 (Pavlis et al., 2012). The national geoid model for Turkish territory, Turkish Hybrid Geoid 2009 (THG-09) (Kilicoglu et al., 2011) was computed depending on EGM2008

The Latest Combined Global Gravity Field Model Including GOCE Data up to Degree and Order 2190
EIGEN-6C4 is a static global combined gravity field model up to degree and order 2190. It has been elaborated jointly by GFZ Potsdam and GRGS Toulouse and contains the following satellite and ground data: -LAGEOS (degree 2 -30): 1985 -2010 -GRACE RL03 GRGS (degree 2 -130): ten years 2003 -2012 -GOCE-SGG data.
-DTU12 ocean geoid data and an EGM2008 geoid height grid for the continents (max degree 370). The combination of these different satellite and surface data sets has been done by a band-limited combination of normal equations (to max degree 370), which are generated from observation equations for the spherical harmonic coefficients.
The resulted solution to degree and order 370 has been extended to degree and order 2190 by a block diagonal solution using the DTU10 global gravity anomaly data grid (Förste et al., 2015).

The Global Gravity Model by Locally Combining GOCE Data and EGM2008
GECO is a global gravity model, computed by incorporating the GOCE-only TIM R5 solution into EGM2008. The input data of GECO: -EGM2008 spherical harmonic coefficients and corresponding error standard deviations -EGM2008 global grid of geoid error standard deviations (5' x 5' resolution) -GOCE TIM R5 spherical harmonic coefficients -GOCE TIM R5 block-diagonal coefficient error covariance matrix. EGM2008 geoid undulations are computed on a global spherical grid of resolution 0.5° x 0.5° by making a synthesis from EGM2008 coefficients up to degree 359. The GOCE geoid on the same grid are computed by making a synthesis from the TIM R5 coefficients up to degree 250. Two geoid grids are merged by least-squares adjustment. Finally, the GECO spherical harmonic coefficients are computed by making an analysis of the combined global geoid grid. The analysis is performed up to degree 359 (consistently with the 0.5° x 0.5° resolution). From degree 360 to degree 2190 the GECO coefficients are the same of EGM2008. The GECO coefficient errors are computed as a weighted average of the coefficient errors of EGM2008 and the TIM R5 solution (Gilardoni et al., 2016).

The Combined Gravity Model GGM05C
GGM05C was estimated to spherical harmonic degree and order 360 from a combination of GRACE and GOCE gravity information (based on GGM05G) and surface gravity anomalies from DTU13. The 2 minute resolution anomalies were used, assuming that they were classical gravity anomalies (i.e., defined on the ellipsoid). The first step was a low pass filter applied to the DTU13 global anomaly field. This was followed by a spherical harmonic analysis of the gravity anomaly set on the ellipsoid, where the coefficients were analytically transformed to degree 540, but only the coefficients up to degree 360 were used. Rather than reprocess the surface gravity data, the full covariance from GGM03C was adopted as apriori. The covariance was then modified so that, below degree 240, the terrestrial information was severely downweighted in order to preserve the accuracy of the GRACE and GOCE gravity contribution. This artificial covariance was used to combine the surface gravity information with GGM05G to obtain the GGM05C solution (Ries et al., 2016).

The Combined Gravity Model GOCO05C
GOCO05C is a static global combined gravity field model up to degree and order 720 based on full normal equation systems (more than 500000 parameters). It has been elaborated by the Gravity Observation Combination (GOCO) Group. GOCO05C is a combination model based on the satellite-only gravity field model GOCO05S and several gravity anomaly datasets (Arctic, Australia, Canada, Europe, Oceans, South America, USA), constituting a global 15'x15' data grid. For the remaining land areas (Central America, Asia, Africa, Antarctica) fill-in datasets (NIMA96, GOCO05S, RWI_TOIS2012) were used (Fecher et al., 2016).
GGMs that are compared over the study area are listed in Table  1 with model characteristics. In GGM evaluation, geoid heights based on GNSS-derived ellipsoidal heights and spirit levelled orthometric heights at discrete points provide an estimated accuracy of the GGM's. The usual and accepted practice is to adopt for a reference model that GGM that is a best fit to the geoid height point estimates determined from the GNSS/levelling. The evaluation of GGMs focuses on the correspondent geoid height differences between the GGMs and GNSS/levelling using the equation below:

Model
where ∆N is the geoid height residual, NGNSS/Lev is the geoid height estimated from GNSS/levelling , and NGGM is the geoid height estimated from GGMs. For the statistical analysis of geoid height differences, minimum and maximum values of ∆N are determined and the overall performance of GGMs is assessed through RMSE accuracy measure defined by: where n is the number of the points used for the accuracy verification and k refers to the residual sequence.

CASE STUDY
For the evaluation process, the geoid heights based on GGMs are interpolated from the closest grid points using software obtained from International Centre for Global Earth Models (ICGEM) web page <http://icgem.gfz-potsdam.de/ICGEM> using the Kriging interpolation method and refer to the reference system World Geodetic System 1984 (WGS84).
The differences between GNSS/levelling based geoid heights and GGM-based geoid heights may be affected by datum inconsistencies. In order to minimize these offsets (i.e. bias and tilt) a 4-parameter transformation is used. The geoid heights obtained from GGMs are compared with discrete geoid heights based on GNSS/levelling data after fitting the tilt. The statistical values of the height data sets that were used for GGM evaluation are given in Table 2 The graphical representations have been adopted for the comparative evaluation of GGMs by producing a residual map for each GGM (Fig. 4-8) that indicates the occurrence and magnitude of geoid height differences. The residual maps are produced by the Surfer  13 software before fitting the tilt.

COMPARATIVE RESULTS AND CONCLUSIONS
The visual analysis of the geoid height residual maps shows that EGM2008 has a better terrain approximation than the other GGMs. It is visible from Fig.4 that the deviation of EGM2008 based geoid heights from GNSS/levelling based geoid heights is reduced for most parts of the study area ( -0.8 m. before fitting the tilt).
Global statistics of geoid height residuals based on GGMs are presented in Table 3. When the statistics summarized in Table 2 are evaluated, the following conclusions can be drawn based on this study: (i) EGM2008, EIGEN6C4, and GECO has better results because of their higher frequency content. (ii) EGM2008 provides more accurate results than other GGMs.
From the statistical values of NGNSS/Lev -NGGM, RMSEs were used to infer the best fit of the GGMs to the GNSS/levelling data for model evaluating because any gravimetric determination of the geoid is deficient in the zero and first-degree terms. Obviously, EGM2008 fit the GPS/levelling data better than other GGMs over the study area. The results of GGM evaluation in this study have indicated the outstanding of EGM2008 to the other GGMs. EGM2008 better statistics than the other GGMs and fits best to the THG-09 at ± 0.2803 m. agreement despite the coefficient errors and GNSS/levelling dataset that can not be considered as an entirely errorless. From our GGM evaluation results we can conclude that EGM2008 can be used as a reference earth geopotential model for further geoid determinations at regional and national scales.

Model
Due to advancements and improvements in instrumentation, software, processes, applications, and understanding, high resolution GGMs (e.g. up to degree and order 2190) are major steps to represent the gravity field of the Earth with a high accuracy. Nowadays global gravity field models, mainly derived from satellite measurements, become more and more detailed and accurate. These gravity field models should be combined with terrestrial gravity anomalies) and GNSS/levelling-derived or altimetry-derived geoid heights. Furthermore, an important task of geodesy is to make the gravity field functionals available to other geosciences. For all these purposes, it is necessary to calculate the corresponding functionals as accurately as possible or, at least, with a well-defined accuracy from a given global gravity field model. Therefore, in order to achieve major improvements for the future high-accuracy gravimetric geoid models in Turkey, further and future analysis of high resolution GGMs (e.g. GOCE-based GGMs) will be needed.