DEVELOPMENT OF A GRAVIMETRIC GEOID MODEL AND A COMPARATIVE STUDY

Site specific geoid model is prerequisite for accurate determination of orthometric heights. No geoid model has been developed so far for India or any of its part. So, development of a geoid model for India or its part is of utmost need to make use of GNSS data towards determination of orthometric heights. In this research work, an attempt has been made to develop geoid undulation models by gravimetric method using Molodensky’s concept. Component parameters in line with the Remove – Compute – Restore (RCR) technique have been used recursively. Models have been developed for two study areas: one of these lies in and around Dehradun (30° 19′ N, 75° 04′E) in Uttarakhand state, India in lower Himalayan region having highly undulating topography and the other near Hyderabad (17° 30′N, 78°30′E) in Telengana state of India having gentle topography. The model has been tested for 7 stations in the first study area and accuracy has been found to be 17.5 cm; whereas, for the second area accuracy has been found to be 7.0 cm for 24 test stations. Further, the performances of the developed models have been evaluated with those from three global geoid models namely EIGEN6C4, EIGEN6C3stat and EGM2008; and have been found to be similar or better in case of first study and for second study area far more superior. Thus, local/regional geoid undulation model requiring accuracy better than 20 cm for any study area may be developed adopting the method. However, the optimality in the number and density of gravity stations may be considered as a future scope of work.


Introduction
Development of a local, regional or global geoid model is one of the forefront research works conducted by geodesists for decades. One of the primary uses of geoid model is to compute geoid undulation or geoid height (Seeber 2003), N which gets subsequently used to determine orthometric height. Thus, in the pursuit for accurate determination of orthometric heights, researchers work towards development of more and more accurate geoid models. As a result, there are many geoid models available in the world, at present. Based on the methods of development, the geoid models are broadly divided into four types (Featherstone et al. 1998). But, gravimetric models have been found to be most widely prevalent.
Gravimetric methods make use of geo-potential model in conjunction with earth's surface gravity of the performance of the developed model with reference to GPS BM data at various regions of USA shows that standard deviation of errors of GEOID09 varies from 6mm to 2 cm whereas USGG2009 provided errors from 2 cm to 7.1 cm (Roman et al. 2010). IRG04, a geoid model for Iran, was developed using Stoke's formula, GRACE GGM02 Global geo potential model and SRTM high resolution Digital Elevation Model (DEM) along with GPS/levelling data. The results from the geoid model were compared with those from GPS and field level of BMs. The RMS error of the model has been found to be 0.09 m which is nearly 4 times more accurate than those obtained from pure gravimetric geoid model (Kiamehr 2007). Carrion et al. (2009) and Srinivas et al. (2012) have reported their attempts towards development of geoid models for the southern part of India using derived gravity data. The results obtained have been compared with those from field levelling (of undefined quality). But, neither known geoid undulation for any station nor the effect of variation of mass density above geoid was considered in these studies. Thus, there are serious deficiencies in the methodology as well as in the validation of the developed geoid models.
Thus, it has been found that, so far, most of the geoid models (except GEOID09) made use of EGM96 as the gravitational potential function which is significantly less accurate (Krynski 2009) than EGM2008. The terrain model used has been derived from satellite data which is less accurate as may be derived from field levelling data. Quality and density of gravity data particularly observed terrestrial gravity data also influence the quality of a model (Hong et al. 2009), which is often not followed properly. Therefore, overcoming of these limitations is expected to lead more accurate geoid models. Thus, scope lies in development of better geoid models using improved geo potential model, accurate digital terrain model and improved (quality as well as denser) terrestrial gravity data. Geoid model also depends on terrain as well as levelling (height) system which varies from region to region. Since, terrestrial gravity as well as levelling data of India has not been shared with any development of global geoid models, so development of a proper gravimetric geoid model for India or its part is of utmost need. In this research work, an attempt has been made to develop geoid models for terrains in India.
Thus the objective of this study is to develop undulation models by gravimetric method for two study areas in India and to carry out a comparative study of the models with some state-of-art global geoid models in the pursuit to validate the developed models for India. The outcome of this study is being expected to provide expressions for geoid undulation, N as a function of latitude, longitude, ellipsoidal height (f, l, h) valid for considered terrains in India.

Background theory
The geoid undulation at any station or point is defined as the difference in elevation between its ellipsoidal height (h) and orthometric height (H) and is given by where N is called geoid undulation or geoid height (Seeber 2003). It depends on anomalous gravity potential (T) at the considered point and may be computed from gravity anomaly (Δg) (Heiskanen, Moritz 1987). However, to make the computation of N independent of density of earth, concept of height anomaly (ζ), as proposed by Molodensky, has also been adopted in this study. According to Molodensky, height anomaly (ζ) is related to the ellipsoidal height (h) and normal height (H*) by the relation (Heiskanen, Moritz 1987): where H * is the normal height. Height anomaly, ζ may be obtained from T γ (Heiskanen, Moritz 1987) where T is the disturbing potential at ground level, γ is the normal gravity at telluroid (Heiskanen, Moritz 1987). The orthometric height of a station or point, H and its normal height, H * may be obtained from where C is geopotential number, g is the mean gravity along plumb line between geoid and ground and γ is the mean normal gravity along the ellipsoidal normal between reference ellipsoid and telluroid. Thus, geoid undulation, N can be reduced to the relation (Heiskanen, Moritz 1987), As both geoid undulation, N and height anomaly, ζ depends on anomalous gravity potential, these are also related to anomalous gravity. Gravity anomaly (Δg) at geoid corresponding to any point/station on the surface of earth can be determined from where g obs is the gravity value observed at a point on surface of earth in gal, γ is normal gravity mathematically computed at ellipsoid by using international Gravity Formula 2 2 1 0.00193185138639 sin 9.7803267714 1 0.00669437999013sin gal and H is orthometric height i.e. height of the point above Indian mean sea level (Indian geoid) in meter. The gravity anomaly (Δg) is being considered to be consisted of three components and is expressed as where egm g ∆ , part of gravity anomaly represented by earth geoid model; rtm g ∆ represents the component due to residual topographical mass; and res g ∆ denotes the residual component of gravity anomaly which varies randomly. Further, each component of gravity anomaly has been considered to provide a corresponding component both for height anomaly, ζ and geoid undulation, N. For height anomaly, components are considered as ζ egm , ζ rtm and ζ res and that of geoid undulation N egm , N rtm and N res respectively. Thus, Height anomaly, ζ = ζ egm + ζ rtm + ζ res ; Geoid undulation, N = N egm + N rtm + N res . (9) All these parameters and their components have been computed from their corresponding gravity anomalies. The value of height anomaly (ζ ) and correction (N-ζ), for stations whose geoid height (N) and orthometric height (H) are known, has been extended to other stations using local covariance function. res g ∆ and ' res g ∆ i.e., residual gravity anomalies at two adjoining points, separated by a distance say "s", are correlated and being a random quantity can be treated statistically. Variance/covariance function, C(s) is defined by C(s) Cov ) { } M( res res res where M has been considered as average, extended over all pairs of points P and P′ for which PP′ = s = constant. The fitting of the empirical covariance functions to local covariance function (analytical models) has been performed using the formula (Arabelos, Tscherning 2003): where n is the degree of expansion of the geopotential model for reduction of gravity anomalies, r, r′ are the distances of the point X, Y from the earth's centre, ˆi σ error anomaly degree variances associated with geopotential model coefficients, R B the radius of Bjerhammar sphere, P i the Legendre polynomial of degree i, y XY the spherical distance between X and Y while A is the gravity anomaly variance. Error degree variances of a geopotential model reflect the behaviour of the model globally, in regional or local computation it is multiplied by a scale factor l which is defined through fitting procedure.

Methodology
In this study, gravity anomaly is being primarily used to as the fundamental parameter to develop the geoid model, using height anomaly as an intermediate parameter. The methodology for developemnt of the proposed geoid model consists of recursive use of component parameters in line with the Remove -Compute -Restore (RCR) technique (Featherstone et al. 2004). First, Δg egm and Δg rtm are to be computed and then remove from Δg (6) to obtain Δg res , at the gravity (observation) stations. These Δg res are to be used to determine the local covariance function and analytical covariance parameters i.e., depth of Berjahamar sphere, (R E -R B ); scale factor, l and variance at zero height (A) of the study area by fitting procedure with global empirical covariance function. Next, control stations (known horizontal location [φ, λ], ellipsoidal height[h], geoid undulation [N] and orthometric height [H]), within the study area, are to be considered. At these points, predicted Δg res are to be computed adopting the least square collocation method using local covariance function with associated parameters of the study area. Further, for these points, Δg egm and Δg rtm are to be computed by linear interpolation. These three components are to be added to find the gravity anomaly (Δg CP ) and subsequently it is to be converted to the height anomaly (ζ) and the correction, (N -ζ). Thus gravimetric geoid obtained is drapped on control points in conformity to the local vertical datum. ζ egm and ζ rtm respectively are to be computed from their corresponding components of gravity anomaly and to be subtracted from ζ to obtain ζ res at the Control Stations. In this way, the framework for geoid undulation model, consisting of model covariance function, covariance parameters and drapped geoid has to be developed for the study area. And, model covariance function, height anomaly and correction alogwith their components on prediction points are also to be made available for determination of geoid undulation, N of any station or point (φ, λ) within the study area. The location of the study areas and those of station positions are as shown in Figure 1. The stations used for development of geoid model have been marked in green while test stations have been marked in yellow colour.

Gravity data
In this study, gravity observations have been taken with relative gravimeters having 1 micro gal resolution, having capability to take 6 readings in 1 second, capable of maintaining standard deviation (s.d.) of 360 readings in 1 minute within 30 µgal and in stable condition, 5 readings taken at a station within 6 minutes agree within 10 µgal. Gravity observations were done at every 3-4 km in plain and at every 1-2 km in hilly terrain. Gravity observations have been taken during the period 2007 to 2014 in International Gravity Standardization Net 1971 (IGSN71). During observation the drift of the instrument has been maintained within 200 µgal. Gravity observation lines were closed at known gravity stations whose gravity values have been determined by absolute gravimeter. In this study, gravity observations were taken at 45 stations in study area 1 and 58 number of gravity data in study area 2 have been used for geoid modeling.

GPS observations and computations
A set of dual frequency geodetic GPS receivers along with geodetic antenna have been used to collect GPS data. In the first study area, a network of 18 stations has been considered for GPS observations in which two stations, namely EE and WE, as reference stations. In order to establish the reference stations, GPS data had been collected in static surveying mode for 24 hours simultaneously on these two stations and precise ephemerides from International GNSS Service (IGS) (ftb://iqscb.jpl.nasa.gov/iqscb/product/) had been used to process the data through a scientific s/w. The maximum RMSE in coordinates (f, l, h) was found to be less than 0.0042 m. For other stations, GPS observations were taken for a few hours in relative mode. Data had been processed in commercial post processing s/w with broadcast ephemerides and elevation mask 10° for baseline processing followed by network adjustment. Adjustment was conducted with 2 reference stations held 'fixed' and setting desired accuracies for position for all. Adjusted coordinates were taken better than 0.011 m in latitude, 0.009 m in longitude and 0.031 m in ellipsoidal height in the entire network at 95% precision confidence level. For the other study area, GPS observations were taken at eight reference stations situated at four corners of each of the two areas (DT06,19,26,31 and HP15,19,24,25) for 7 hours at 15 seconds epoch interval having GDOP less than 5. These GPS data (c) Study Area 2 -Telangana near Hyderabad were processed using scientific software with precise ephemerides with respect to IGS station at IISc, Banglore. The maximum RMSE in coordinates (f, l, h) was taken less than 0.0031 m. For other 54 GPS stations, observed were taken each of 45 minutes session at 15 second epoch interval. The observed data were processed through commercial software considering the reference station 'fixed' . Networks of 32 stations of one part and 26 stations of the other part were adjusted separately with maximum error 0.0696 m in coordinates at 95% confidence interval.

Level data
Orthometric heights of eighteen stations in first study area and fifty eight stations in the second study area were taken from SOI [Survey of India] archives. The accuracy of measurement has been controlled by the permissible discordance between forward and backward levelling between two consecutive BMs and also of a levelling line to be less than 4√k millimetres where k is the direct distance traversed in kilometres. Further the observed geometric height difference of a circuit and individual lines have been statistically analysed by Lallemand formula (Avers 1935). The probable accidental error and systematic error have not exceeded ±1.0 mm/km and ±0.2 mm/km respectively. The mean accidental and systematic error have not exceeded ±1.5 mm/km and ±0.3 mm/km respectively. Geometric height data in conjunction to gravity data have been used to derive geopotential numbers which have been used to obtain Helmert orthometric height (Heiskanen, Moritz 1987).

Digital Terrain Models
For study area 1 i.e., around Dehradun, DTM has been prepared from the Survey of India contour maps of 1:50,000 scale and for study area 2 i.e., around Hyderabad, shuttle radar terrain DTM (3") called as SRTM3" has been used.

Development of geoid model
In this study, GRAVSOFT software (Tscherning et al. 1992;Tscherning, Forsberg 2008) has been used to develop the geoid model. The Remove/Compute/ Restore operations have been implemented repeatedly to operate on different modules of the software. Gravity observations and known geoid undulations of 45 and 11 stations respectively in the study area 1; and that of 58 and 34 stations respectively in study area 2 have been used during development of these models. Further, the Digital Terrain Models (DTM) of the study areas (4.4) and Global Geoid Model, EGM2008 have been used as primary inputs to the GRAVSOFT. EGM2008 which uses ellipsoidal harmonics (of degree 2190 and order 2159) has been used as base model. GEOEGM and TC modules of Gravsoft has been used to interpolate free air gravity anomaly represented by ellipsoidal harmonics of Global Geoid Model (GGM) (Δg EGM ) and terrain effect (Δg RTM ) respectively at observed gravity stations. These two have been removed from free air gravity values (6) to obtain residual free air gravity (Δg res ). The digital terrain data at 15′′ and 60′′ grid have been calculated as reference and coarse data by TC. Covariance function of residual gravity of the region has been developed by running EMPCOV with residual gravity anomaly (∆g-egm-tc.dat) as input. The output of this execution has been used in COVFIT program. Values of depth of Bjerhammar sphere i.e. RE-RB, scale factor, l and gravity variance, A, to be used in GEOCOL module, have been noted through command prompt. These values for Dehradun area have been found to be 0.0759, -424.22 and 150.81 respectively for scale factor, RE-RB and gravity anomaly variance and for Hyderabad area 0.0176, -200.55 and 1.51. SELECT module with GPS data as input was then used to find residual gravity anomaly at GPS points within the study area having known geoid undulation (N). Residual free air gravity anomaly have been computed at the GPS stations by running GEOCOL program with Δg res and GPS data as input files. With GPS data as input file, by running GEOEGM, free air gravity anomaly component represented by Global Geoid Model (GGM) has been found. Similarly by running TC program with GPS data as input, free air anomaly component due to terrain effect has been found. Both are added to residual free air anomaly at GPS points by running FCOMP program to find free air gravity anomaly at GPS points. By running N2ZETA program using free gravity anomaly at GPS points with known N, height anomaly (ζ) and N-ζ at GPS stations have been obtained. With operations so far, the derived geoid has been drapped in local orthometric height terms. A flow chart showing the different steps, inputs, modules for remove/compute/restore operations is shown in Figure 2. With these, the geoid undulation model gets developed. The geoid undulation, N of any station or point (j, l) within the study area may now be computed using model covariance function and remove/compute/add/restore procedure.

Terrestrial gravity anomaly
To study the terrestrial gravity anomaly (∆g) and its variation as well as variation of its components in the study area 1, statistics of the observed terrestrial gravity anomaly such as maximum value, Minimum value, Mean and Standard Deviation was computed. Same were computed after removing contribution from GGM (∆g -∆g egm08 ) and also for residual gravity anomaly (∆g res ) . The result has been summarised in Table 1.

Covariance functions
In order to study the nature of the covariance function developed for the study area 1, empirical covariances against distances have been plotted, as shown in Figure 3. Also, for the same distances, the analytical model covariance values have also been plotted. It has been found that there are large variations between the empirical and model covariance functions. Next, empirical co-variances against distances have been plotted, as shown in Figure 4 and also those by the analytical covariance function of the study area 2. In this case, the two functions have been found to be well fitted.
It is seen that value of residual gravity is very significant in this area. Variation of gravity anomaly is also high which in in conformity to the highly undulating terrain.
Similarly statistics of ∆g, ∆g -∆g egm08 and ∆g res for study area 2 were also computed and summarised in Table 2 below. As expected, being a flat terrain, variation of gravity anomaly is small in study area 2. Residual gravity anomaly is also very small in this area.
It is evident from Table 1 and Table 2 that EGM2008 represents the major portion of gravity anomaly but it is unable to represent it fully due to limitation of degree and order of spherical harmonic coefficients. This limitation is more pronounced in undulating terrain.

Geoid undulations or Geoid height, N
The geoid heights, N of the test stations determined by using Equation (1) ellipsoidal height (h) (from GPS observation) and precise levelling height (H) have been considered as true geoid undulations and thus, used for comparative evaluation.
In order to test the geoid models developed in this study, first the geoid undulations (N GRAV ) of 7 stations in the first study area have been computed using the developed model. To determine the geoid undulation of a station point, first the components of gravity anomaly at the desired points have been computed by method of interpolation and least square collocation. These are then converted to their corresponding components of height anomaly (ζ) as well as the correction factor, (N-ζ) for height anomaly. Finally, these are converted to components of geoid undulations i.e., N res , N egm , N rtm and are subsequently added to find the geoid undulation (N GRAV ) at desired location. The computed geoid undulations are then evaluated with respect to the true geoid undulations. It may be noted that the mean error in model based geoid undulation (N GRAV -N) has been found to be 0.175 meter and that RMSE is 0.190 meter. Further, these geoid heights have been compared with those from some standard models: EIGEN6C4 (N EIGEN6C4 )(www.gfz-potsdam. de), EIGEN6C3stat (N EIGEN6c3 ) (www.gfz-potsdam. de) and EGM2008 (N EGM2008 ) (www.icgem.gfz-potsdam.de/ICGEM/modelstab.html). The mean errors of N EIGEN6C4 , N EIGEN6c3 and N EGM2008 have been found to be 0. 175 meter,0.187 meter,0.197 meter respectively;and those RMSEs are 0.190meter,0.197 meter,0.208 meter respectively. The deviations of these geoid heights have been listed under columns 2, 3 and 4 i.e. columns headings (N EIGEN6C4 -N), (N EIGEN6c3 -N) and (N EGM2008 -N). The entire result for study area 1 has been summarised in Table 3.
For the study area 2 i.e., the study area around Hyderabad, 24 stations have been considered for testing purpose. Geoid undulations (N GRAV ) of these 24 stations have been computed using the developed model of the second area and these values have been   N). The entire result for study area 2 has been summarised in Table 4.

Discussion
The empirical covariance function of the study area 1 has been found to be deviated much from its analytical covariance model. As the physiography of the area is much undulating and area belong to lower Himalayas i.e., the area is under the influence of huge mass accumulation around, there is much variations in gravity anomalies at different stations and thus, its empirical covariance function. On the other hand, the empirical covariance function of the study area 2 has been found to be fitted with its analytical model, This is commensurate with the physiography of the Hyderabad area which is quite flat and associated with no accumulated mass around. Thus, the deviations or fitting of the derived covariance functions with those from their analytical models are as expected to their field conditions.
The accuracy of the model for the first study area has been found to be 17.5 cm. However, for different test stations variations are quite prominent with errors varies from 26.6 cm to 2.3 cm; these are attributed to the great variations in orthometric heights of the stations. Further, the presence of Himalayas around the stations deteriorate the field condition for gravity anomaly and thus, in the adopted method. Further, errors in geoid undulations determined by global geoid models namely EIGEN6C4, EIGEN6C3stat and EGM2008 have been found to be same or more than what has been found by geoid model. For the second study area, the accuracy of the model has been found to be 7.0 cm i.e., far superior to first model. Moreover, errors in geoid undulations determined by EIGEN6C4, EIGEN6C3stat and EGM2008 have been found to be 45 cm, 50 cm and 42 cm respectively which are far more than what have been found by developed model. Thus, it can be found that the gravimetric geoid model developed is as accurate as EIGEN6C4 for the study area 1 but for study area 2 developed model is five times more accurate than EIGEN6C4. The superiority may be attributed to the nature of terrain and/or consideration of large number of stations for development of the model. The gravimetric models developed in this research works have been found to be better than all the three best global geoid model for both study areas. This is due to consideration of Molodensky's concept for the determination of N as well as use of true N of some stations in the area for drapping of geoid in local vertical terms. However, by increasing the density and number of gravity stations as well as quality of gravity data still better covariance model may be developed, of any study area. So there is a scope for further improvement of the developed models.

Conclusions
In this study, geoid undulation models have been developed by gravimetric method for two different types of topography. The developed models have also been validated with some global models and have been found to perform better than the best global models. The method consumes less time, involves less efforts and more economic than geometric method for development of geoid model. It appears that adopted gravity method has the potential to provide most economical local or regional geoid model. Thus, local/regional geoid undulation model having accuracy better than