An experimental Indian gravimetric geoid model using Curtin University’s approach

University’s approach R. Goyal, W.E. Featherstone, S.J. Claessens, O. Dikshit, N. Balasubramanian 1) Department of Civil Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India; rupeshg@iitk.ac.in, onkar@iitk.ac.in, nagaraj@iitk.ac.in 2) School of Earth and Planetary Sciences, Curtin University of Technology, GPO Box U1987, Perth WA 6845, Australia; w.featherstone@curtin.edu.au, s.claessens@curtin.edu.au


INTRODUCTION
Ideally, Stokes's (1849) integral should be implemented over the entire Earth with continuous gravity anomalies on the geoid and with the condition that there must not be any gravitating masses above it.
However, in practice, the availability of gravity observations is limited to a specific area, so the integration domain has to be truncated. Also, the gravity anomalies usually exist discontinuously on or above the Earth's surface so various types of downward continuation and regularisation have been proposed. The gaps between theoretical and practical aspects induce several kinds of errors, which geodesists have tried to reduce, but usually requiring assumptions and approximations.
Based on various ideas, philosophies and numerical approaches, what we consider the four most commonly used approaches/techniques are adopted for geoid computation experiments in India.
For India, the first geoid map was developed more than five decades ago. It was based on astrogeodetic observations (Fischer 1961) and with respect to the Everest 1956 ellipsoid (cf. Singh andSrivastava 2018). No more information is available on this geoid, apart from distorted hardcopy contour maps that are difficult to digitise reliably. The levelled height information presently available < A c c e p t e d M a n u s c r i p t > Page 4 of 38 in India is more than a century old. When these heights were observed, neither the concept of foresight and backsight levelling nor the use of invar staves were considered. Observed gravity values were not available as this was before the development of the low-cost portable relative gravimeter. The Indian vertical datum defined in 1909 was based on constraining the levelling to nine tide-gauges along the Indian coast to zero height (Burrard 1910). We will show later that this approach may have caused a north-south tilt (cf. Fischer 1975;1977), most probably due to the ocean's time-mean dynamic topography (cf. Featherstone and Filmer 2012).
Frequent seismic activity in various parts of the Indian sub-continent and so-caused crustal movement also necessitate the introduction of a new height system, probably to be based on geopotential numbers and Helmert's orthometric heights (or 'rigorous' orthometric heights as formulated by Santos et al. 2006). The Survey of India (SoI) carried out a re-levelling program [2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017] with gravity observations at fundamental benchmarks to provide a densified network of Helmert's orthometric heights as a part of the Redefined Indian Vertical Datum 2009 (Singh 2018;G&RB 2018). However, these data are not yet in the public domain, so we are unable to use them to validate our geoid and quasigeoid models. In addition, the national geodetic agencies have proposed to compute a precise national geoid model to serve as the new vertical datum for the country. This can be viewed as following the suit of New Zealand (LINZ 2016), Canada (Véronneau and Huang 2016), and the USA (NGS 2017;. Such an approach is being considered in many other countries too. Researchers and government organisations have made some efforts to develop local gravimetric geoid models for regions in India (Singh 2007;Carrion et al. 2009;Srinivas et al. 2012;Mishra andGhosh 2016, Singh andSrivastava 2018), but only using the GRAVSOFT package with residual terrain modelling (Forsberg 1985). Despite these efforts, a state-of-the-art national gravimetric geoid model for the whole India remains elusive (Goyal et al. 2017). Therefore, in this study, we present the first-ever nationwide geoid and quasigeoid computation results over India with < A c c e p t e d M a n u s c r i p t > Page 5 of 38 the available data sets using the CUT method implemented using our own computation package developed in MATLAB ™.

Terrestrial Gravity
Pointwise observed gravity data is confidential in India. Therefore, with this predicament, we obtained a grid of Indian terrestrial gravity data from GETECH (https://getech.com/) that is claimed to come from the Gravity Map Series of India (GMSI), a joint project of five Indian organisations, viz., SoI, Geological Survey of India (GSI), Oil and Natural Gas Corporation (ONGC), National Geophysical Research Institute (NGRI) and Oil India Limited (cf. Tiwari et al. 2014). The GETECH gravity data comprises a 0.02˚× 0.02˚ grid of simple Bouguer gravity anomalies over India (except a for few regions in northern India), with an overall estimated precision of ± 1.5 mGal (GETECH, 2006). According to the GETECH manual for Indian gravity data, they used i) the normal gravity formula from WGS84 (Macomber 1984) (Sandwell et al. 2019). The SCRIPPS data is also accompanied with an error grid that we have shown, for our study area, in Figure 1. The data contains a 1'x1' grid that also covers the land, but we used the SCRIPPS data only for the oceanic region because the land data, in the SCRIPPS dataset, is from EGM2008 to avoid aliasing (Gibbs fringing) at the coasts.

PLACE FIGURE 1 NEAR HERE
We do not have gravity data from the countries neighbouring India and a well distributed As discussed next, we merged these three datasets to get a complete free-air gravity anomaly grid of 0.02°x0.02°, avoiding aliasing or the contamination of land data (both GETECH and EGM2008 individually) with the marine data or vice-versa.
There exist numerous sophisticated space-domain and frequency-domain methods for merging heterogenous gravity anomaly datasets (e.g., Strykowski and Forsberg 1998;Olesen et al. 2002;Catalao 2006;McCubbine et al. 2017). However, we chose to work with the comparatively

< A c c e p t e d M a n u s c r i p t >
Page 8 of 38 straightforward CUT space-domain method (cf. Featherstone et al. 2011;. This choice is somewhat arbitrary because we are working with the land gravity of unknown quality, and the strategy that we use has already been implemented in the computation of the Australian quasigeoid, which is an island nation and approximately 2.3 times larger than India. Other methods can also be tested, but it is left for the time when sufficient marine and airborne gravity data along with reliable terrestrial gravity data will be available over India. In the adopted method, the GETECH-derived free-air anomaly grid is superimposed over the EGM2008 (d/o 900) derived gravity anomalies. The gravity anomalies of the latter dataset at the overlapping grid nodes are replaced by the gravity anomalies from the former dataset. As a result, a 0.02°x0.02° grid of gravity anomalies on the land is obtained.
To concatenate the land and marine gravity anomaly data, 1'x1' gravity anomalies in the ocean are clipped (or separated) from the complete SCRIPPS dataset, i.e., on both ocean and land. It is then block averaged to the 0.02°x0.02° grid and is superimposed with the land gravity anomaly grid. The former values were replaced by the latter at overlapping nodes to obtain the 0.02°x0.02° grid of the merged gravity anomalies. Figure 2 shows the merged free-air gravity anomaly map. To check for any discontinuities at the edges of the merged datasets, we computed and plotted the arctangent ( Figure 3a) and logarithmic ( Figure 3b) values of the gradients of the merged data. We observe no clear visual indication of any discontinuities at the boundaries of the merged data, but also partially due to the ruggedness of the dataset in our study area that can be obscuring.

Digital Elevation Model
The Digital Elevation Model (DEM) is another important input in geoid computation. It is mainly used to compute the topographical effects (e.g., Forsberg 1984). Thus, a precise high-resolution DEM should be used. We would like to mention here that DEM is generally used synonymously with a < A c c e p t e d M a n u s c r i p t > Page 9 of 38 Digital Surface Model (DSM) (e.g., SRTM, ASTER), but this should be avoided. Quantification of the differences in the topographical effect with the use of DEM and DSM has been investigated by Yang et al. (2019). Since India does not have a national DEM, therefore, after a DEM/DSM analysis (Goyal et al. 2021a), it was decided to work with the best available DEM over India, i.e., the MERIT 3"x3" DEM (Yamazaki et al. 2017), for our computations. Though the accuracy of the MERIT DEM varies considerably (± 11.7 m to ± 47.3 m) over different landforms in India, an overall estimate for the whole of India is ± 17.3 m (Goyal et al. 2021a).

GNSS-Levelling
India has different horizontal and vertical control networks. Therefore, presently there are only a limited number of ground control points where we have the geodetic coordinates (latitude, longitude, ellipsoidal heights) and levelled heights. Moreover, due to several restrictions on the datasets, only a few of these available data points were available to us (Figure 4). The datasets in the Uttar Pradesh west (UP west) and Uttar Pradesh east (UP east) regions were procured from SoI, while the datasets over Hyderabad and Bangalore have been retrieved from Mishra (2018), who also used the SoI dataset.
According to Mishra (2018), horizontal and vertical precisions of GNSS data are within ±12 mm to ±26 mm and ±31 mm to ±53 mm, respectively. The vertical precision of the levelling heights is not known to us, but they are from the high precision first level net of India. These heights are from the Indian Vertical Datum 1909 (Burrard, 1910) and are based on the normal-orthometric height system, while those on Indian Vertical Datum 2009 (G&RB 2018) are based on Helmert's orthometric height system. We have not been provided with a clear indication on which heights have been provided to us, and therefore, due to this anonymity of the height system, we consider the levelling heights to be in the normal-orthometric height system (Jekeli 2000;Featherstone and Kuhn 2006).

< A c c e p t e d M a n u s c r i p t >
Page 10 of 38

METHODOLOGY AND RESULTS
An overview of the CUT methodology for computing the geoid undulations is shown by a flowchart in Figure 5. The CUT method primarily computes the quasigeoid using the analytical continuation solution (Moritz 1971;) of Molodensky's problem (Molodensky et al. 1962). Moritz (1971) showed that Molodenksy's G1 term can be approximated by the planar terrain correction (TC), which also needs an additional term that is equal to the first-order indirect effect (FOIE). We could not adopt the full CUT method-based reconstruction of Faye anomalies (Featherstone and Kirby 2000) because we already have gridded data, whereas CUT grids point Bouguer anomalies. Instead, we added the block averaged 0.02˚× 0.02˚ grid of TCs (Figure 6a) to the free-air gravity anomaly grid to calculate area-mean Faye anomalies (Figure 6b). The block-averaged 0.02˚× 0.02˚ TC grid was constructed from the 3"× 3" TC grid computed with the MERIT DEM using the Optimal Separating Radius (OSR) in the spatial-spectral combined method suggested by Goyal et al. (2020). This method of TC guarantees the full convergence of the TC solution, i.e., down to <1μGal.

PLACE FIGURE 5 NEAR HERE
A different approach is used in the CUT method to apply ellipsoidal correction. Unlike other geoid computation strategies considered (UNB or KTH; cf. Huang et al. 2003;Ellmann 2005), the CUT method computes ellipsoidal area-mean free-air gravity anomalies on the topography using a Global Geopotential Model (GGM) (Featherstone et al. 2018). These are subtracted from the mean Faye gravity anomalies to obtain residual gravity anomalies (Figure 6c), which are then Stokesintegrated with the Featherstone-Evans-Olliver (FEO) modified kernel ) to obtain the residual height anomalies. The FEO kernel, a deterministic modifier, is the combination of the Meissel (1971) and Vaníček and Kleusberg (1987) modifiers that simultaneously reduces the truncation error and improves the rate of convergence to zero of the series expansion of the truncation error (cf. , Featherstone 2003. Additionally, the spherical reference radius in

< A c c e p t e d M a n u s c r i p t >
Page 11 of 38 the Stokes integration is set equal to the geocentric ellipsoidal radius of the computation point, and this negates the need for further ellipsoidal corrections to Stokes's integral (Claessens 2006   in Heiskanen and Moritz (1967, p 328). The difference in the quasigeoid-geoid separation term from the two methods is shown in Figure 6f. Acknowledging that the number and distribution of the GNSS/levelling data points are not sufficient for reliable fitting (Kotsakis and Sideris 1999; where n is the total number of discrete GNSS/levelling data points. It is important to acknowledge that absolute accuracy is only an assumption. This is principally because the levelled heights that refer to the local vertical datum are not necessarily coincident with the geoid. This has been discussed in detail by Featherstone (2001). The descriptive statistics of abs i  are in Table 1.

< A c c e p t e d M a n u s c r i p t >
Page 13 of 38 The relative testing of geoid and quasigeoid (Eqn. 11) is an analysis tool to investigate their gradients. This type of analysis is of more interest to land surveyors who use relative GNSS baselines and a geoid/quasigeoid gradients as a replacement for the more time-consuming differential levelling.
  , 1, 2,3,....... ; The descriptive statistics of rel ij  ), and the ratio of mean differences to the mean baseline length in parts per million (average ppm in mm/km) for the geoid and quasigeoid are in Table 2. PLACE

DISCUSSION, RECOMMENDATIONS AND CONCLUSIONS
Though the number (119) and the distribution (Figure 4) of the GNSS/levelling data points are insufficient to draw concrete conclusions about the quality of the computed gravimetric geoid and quasigeoid models, the following are some major observations from our experimental results: 1. Since the study area comprises the most complex topography varying from the Himalayas to the Gangetic plains and a long peninsular coastline, Figure 6 possibly depicts the extreme (maximum and minimum) values of planar TC, Faye gravity anomaly, and quasigeoid-geoid separation on the planet.
2. From the viewpoint of the "cm-level accurate" geoid, Figure 6f suggests that a more rigorous method (e.g., Flury and Rummel 2009) should be preferred for calculating the quasigeoid-geoid separation over a simple approximate formula (e.g., Heiskanen and Moritz 1967). There exist other formulas for the quasigeoid-geoid separation term (e.g., Sjöberg 2010; Foroughi and Tenzer 2017), but they are not tested here. 4. Generally, standard deviations versus GNSS/levelling are large for lower modification degrees and larger integration radii (Featherstone et al. 2018;Claessens and Filmer 2020). However, Figure 7 shows an opposite trend in India, with smaller standard deviations for lower modification degrees and larger integration radii. This is primarily attributable to the north-south tilt in the India height datum (cf . Table 1). However, he smaller number of GNSS/levelling data and their poor distribution are also likely to contribute to this observation. Figure 7 shows that the Indian levelling heights are marginally better referred to the quasigeoid (std = ± 0.389 m) than the geoid (std = ± 0.396 m). However, Table 1 shows that the geoid has an equal or better precision estimate than the quasigeoid (in terms of standard deviation) in each of the four regions individually. The difference in the standard deviations of the quasigeoid and geoid comparison for the whole of India seems to be a consequence, mostly, of the smaller mean < A c c e p t e d M a n u s c r i p t > 6. Table 2 indicates that the largest misclosures in Figure 8 are probably due to the tilt in the Indian height datum and the relative closeness of data points in Hyderabad and Bangalore, which also Bangalore, an improvement in the standard deviation beyond ± 0.01 m is observed only for UP east. However, the standard deviation of gravimetric geoid in UP west is degraded by ± 0.03 m as compared to the EIGEN-6C4. A degradation in the standard deviation of the gravimetric geoid is also observed in Featherstone and Sideris (1998). This was, and similarly is, attributed to errors in either one or more of terrestrial gravity data, GGMs and the GNSS/levelling data.

5.
There is little to no improvement with the inclusion of the terrestrial gravity data with the CUT method because it makes use of the highest available degree-order GGMs. Also, the GETECH data is possibly already included in the high degree-order GGM (e.g., EGM2008, Pavlis et al., 2012;2013)

< A c c e p t e d M a n u s c r i p t >
Page 16 of 38 8. The Faye gravity anomaly (Figure 6b), geoid (Figure 9a), and contour map (Figure 9b)  other geoid/quasigeoid computation strategies should also be tested over India.

Acknowledgments
We are thankful to the GETECH Pty. Ltd. and the Survey of India for providing the gravity data and GNSS/levelling data, respectively. The Indian SPARC scheme is thanked for providing partial funding to procure the GETECH gravity data. The National Centre for Geodesy, established at Indian Institute of Technology Kanpur, is duly thanked for partial funding to procure the GETECH gravity data, and full funding for a) procurement of the GNSS/levelling data and b) the additional page

< A c c e p t e d M a n u s c r i p t >
Page 17 of 38 charges for this article. We are also thankful to the three anonymous reviewers and the editor (X. Li) for their prompt and constructive comments on an earlier version of this manuscript.

< A c c e p t e d M a n u s c r i p t >
Page 37 of 38

< A c c e p t e d M a n u s c r i p t >
Page 38 of 38