Regional Geoid for Nepal using Least-Squares Collocation

Regional Geoid for Nepal using Least-Squares Collocation Sushmita Timilsina , Martin Willberg, Roland Pail 1 Ministry of Land Management, Cooperatives and Poverty Alleviation, Land Management Training Centre, Dhulikhel, Kavrepalanchok, Nepal 2 Institute of Astronomical and Physical Geodesy, Technical University of Munich, Munich, Germany This manuscript is submitted to “Terrestrial, Atmospheric and Oceanic Sciences”


INTRODUCTION
Geoid is a complex surface and the most relevant physical figure of Earth. Its significance is comprehended from former studies and its determination has been performed for several decades in engineering and science (Bernhard& Moritz2005). The basic task of the geoid in physical geodesy is to serve as a reference for orthometric heights, which are usually determined by spirit levelling.
However, the error in spirit levelling increases with distance of the levelling network, so levelling requires a lot of effort, time and money.
The huge progress made in Global Navigation Satellite System (GNSS) to determine horizontal and vertical positions on Earth's surface has provided an opportunity to determine the height of any point above the reference ellipsoid everywhere on the Earth. However, the reference ellipsoid is a smooth mathematical approximation of physical Earth's surface and does not coincide with its true figure or to the geoid surface. If proper determination of the geoid height (height of geoid above ellipsoid measured along the ellipsoidal normal) is possible, then any ellipsoidal height h measured on Earth's surface above reference ellipsoid could be converted to an orthometric height H (height of a point P on Earth's measured above geoid along the plumb line) by applying . (1) Determination of the geoid height allows us to transform ellipsoidal heights to physical heights, which are associated with Earth's gravity field. Hence, the geoid is used in mapping and engineering works.
Being home to more than eight highest mountains in the world with elevation ranging from about 60 meters to 8848.86 meters (Nepal Gazette, 2021), geoid determination for Nepal is a crucial job.
This study aims to determine a regional geoid for Nepal employing airborne gravity data collected in the year 2010 with the application of a recent high-resolution GGM (Global Geopotential Model), i.e. XGM2019e (Zingerle, Pail, et al, 2019). XGM2019e is a combined global gravity field model represented through spheroidal harmonics up to d/o 5399, corresponding to a spatial resolution of 2' (~4 km). Determination of geoid height from the reduced airborne gravity anomalies was done using Least Squares Collocation with the proper determination of covariance functions (Moritz, 1980). Least Squares Collocation takes the statistical behaviour of the gravity quantities through the covariance function into account to interpolate/estimate the geoid height in least square sense. The comparison of thus determined geoid, named as NPG20, was compared with previously computed geoid NPG11 (Forsberg et al, 2014), where NPG11 was computed using same airborne gravity data with the application of spherical Fast Fourier Transformation (FFT). The comparative analysis of these models shows that NPG20 in the framework of XGM2019e has proven to be a better geoid model for Nepal.

< A c c e p t e d M a n u s c r i p t >
Page 4 in 18 pages

Data
The gravity dataset collected from the airborne gravity survey of Nepal in 2010 in cooperation between DTU-Space, Survey Department, Nepal and NGA, USA is used as the main input for this study. The main objective of that survey was to provide data for a new national geoid model, which will in turn support Global Navigation Satellite System (GNSS) surveying and national geodetic infrastructure, and provide gravity information for future global gravity field EGM2020 (Barnes, Beale,et al, 2015). The overall accuracy of the collected airborne gravity data was estimated to be 3.3 mGal (Forsberg, 2013). Figure 1 shows the observed gravity anomaly in mGal along the flight line. The gravity anomalies range between -200mGal to 350mGal. The high variation in gravity anomalies is the result of high variation in elevation range (~60-8848m) (Manandhar, 2010). Gravity observations in some flight lines are not complete due to high turbulences during the aircraft flight.
Three Global Geopotential Models EGM2008 (Pavlis,et al,2008), XGM2018 (Pail, 2018, and XGM2019e (Zingerle, Pail, et.al, 2019) were used for this study. The models EGM2008 and Field Model XGM2018 (Pail, et al, 2018) was developed with the optimal combination of the new and improved terrestrial data set of 15′×15′ area-mean gravity anomalies provided by the NGA, USA with the latest satellite gravity information (GOCO05s) (Pail, et.al, 2020). XGM2018 is complete up to spherical harmonic degree and order 760. XGM2019e is a combined global gravity field model that includes recent satellite model GOCO06s in the longer wavelength range with terrestrial measurements over land and ocean of gravity anomalies provided by NGA (identical to XGM2016, having a resolution of 15') for the shorter wavelengths. The terrestrial data itself is augmented with topographically derived gravity over land (Earth2014) (Zingerle et al, 2019).
XGM2019e is complete up to spherical harmonic degree and order 5539. The comparison for the best fitting global model was done by statistical analysis.

< A c c e p t e d M a n u s c r i p t >
Page 5 in 18 pages

Evaluation of input data and derivation of covariance functions
The determination of regional geoid for Nepal was done using Least-Squares Collocation approach in Remove-Compute-Restore framework. For removal of the long wavelength component, the GGM that results in the best statistical fit to gravity observations is considered the most suitable for modelling of the long-wavelength signal of the gravity field (Zhang 1997). Therefore, the removal or reduction of the long-wavelength component from observed gravity anomaly was done by using the best possible global model. The statistical analysis of the residuals of the airborne gravity data after the reduction from these global models is presented in Table 1.  Table 1 suggests that XGM models (2018 and 2019e) fit better to the input data than EGM2008.
Among the XGM models, XGM2019e provided slightly better results than XGM2018. The terrestrial data in XGM2019e is augmented with topographically derived gravity over land (Earth2014), and it is complete up to spherical harmonic degree and order 5539 (Zingerle et. al, 2019). So, the medium and long wavelength components of the gravity anomaly are reduced in a single step. Therefore, XGM2019e was chosen for further computations in this study.
The reduction step of gravity anomalies using XGM2019e undergoes following (2) where is the input observed gravity anomaly, is the gravity anomaly computed using XGM2019e at same position as the input data, and is the reduced gravity anomaly. The reduced gravity anomaly is presented in Figure 2. The reduced gravity anomalies are now in the range of ± 60mGal. After the reduction of gravity anomalies using Eq. 2, the determination of covariance function for least-squares collocation was done. A covariance function that fits the empirical gravity anomalies, i.e., Empirical Covariance Function (ECF) with correlation length of 4827.2 m, and a Model Covariance Function (MCF) derived from degree variances of the global model (XGM2019e) were computed and plotted against the distance (Fig. 3). The upper and lower limits of the spectral band of the XGM2019 variances were determined empirically, with the goal to obtain the best fit to the correlation length of the empirical covariance function. This is important, because the correlation length is the driving parameter for the spectral behaviour of LSC. The degree range for the model covariance was chosen so that it matches best with ECF. Here, degree range 455 to 2000 shows the best fit to for the empirical data. There was not much change in the covariance function for the degrees above 2000. In a second step, the variances (and thus the model covariance function) were scaled by a factor of 5.93 to obtain the same variances (covariance values at distance zero) for model and empirical covariance functions. This can be phyiscally justified, because the signal variance in the Himalaya region is significantly higher than the global average, which is represented by the XGM2019 variances. It should be emphasized that the scaling does not have an impact on the collocation result, but just on the (scaling of the) formal error estimates. Figure 3 shows the comparison of the Empirical Covariance Function and the Model Covariance Function.

Figure 3: Comparison ECF (Computed from input gravity anomaly) and MCF (Computed from XGM2019e) for degree
range 455-2000. The agreement between ECF and MCF is an important part of collocation. However, it is not possible to find a satisfactory fit for the behaviour of ECF. So, we focus on the realization of a realistic half width at half maximum and accept the differences for higher distances. Figure 3 shows quite good agreement of MCF and ECF for distance less than 20 km, while the fit around 40 km can be certainly improved. This MCF was then used as the covariance function to proceed to Least Squares Collocation.

Least Squares Collocation
Least Squares Collocation (LSC) has been widely applied in geodesy for estimating the gravity field of the Earth both locally and globally (Gaetani, 2016). The main aspect of this method is the statistical interpretation of proper covariance functions of the gravity data, which describe the spatial correlation of the observations, as a kernel function.
All the gravity quantities can be derived from the disturbing potential as the basis gravity functional, because there exists a relationship between them. For example, if T is the disturbing potential γ is the normal gravity, then the geoid height can be calculated as N=T/γ. For calculating the geoid height with input disturbing potential T, 1/ γ is the linear operator that links these two gravity quantities. Similarly, the covariance of each derived quantity can also be calculated using the same linear operator following covariance propagation rule.

< A c c e p t e d M a n u s c r i p t >
Page 8 in 18 pages Any measurement of gravity quantities like any other physical quantities consists of a mathematical model that defines the physical quantity, the signal component and noise component. Least Squares Collocation considers all these components of physical quantities (Ruffhead, 1987).
If is the signal and is the noise in the observation of , and is the linear functional of disturbing potential or geoid height, then the basic observation equation of LSC is given as follows, , where is the cross-covariance of estimated signal and input signal and sum of covariance matrices of noise and signal. These covariance matrices are based on the model covariance function.
Therefore, the choice of the model covariance function is of great significance.
The error of estimation is given by (4) where is a priori covariance matrix and is the a posteriori covariance matrix. The error estimate relates predicted signal quantities with input signal through cross-covariances of input and output signals ( ) and autocovariances of input signal .
The basic task of LSC is the prediction of signal quantities where the measurement has not been done with or without the transformation of input signal along with the optimal removal of noise.
LSC has been exercised in this study for gravity signal transformation, i.e. to determine geoid heights from gravity anomalies along with its error estimation.

< A c c e p t e d M a n u s c r i p t >
Using the reduced gravity signal (reduced by XGM2019e) at ellipsoidal coordinates (latitude, longitude, elevation) and MCF, LSC was then employed to compute geoid height from the input gravity anomalies dataset of 11105 points on a regular interval of 3'x3'. The interval was chosen as such in order to make it comparable with the previously computed geoid model. The latitude ranges between 26 degrees to 30.50 degrees and longitude ranges between 80 degrees to 88.50 degrees.
The boundary of Nepal fits properly inside these intervals. The geoid residuals obtained from LSC as presented in Figure 4 range in between ± 1 m inside the country's borders. Higher values are seen in the major mountainous regions of the country. Due to the lack of terrestrial data at the border and outside the country's border, LSC cannot derive realistic results, resulting in higher residuals ± 2.5 m.
The error estimates of the output of LSC are presented in Figure 5. They relate the input gravity anomaly points with output geoid height values through cross-covariance between input ∆g and output and auto-covariance of input ∆g.The error estimate of LSC in Figure 5 shows that the standard deviation in the estimation of geoid values ranges between 2 cm to 20 cm. The deviation is higher in the region where there is no input data, i.e. outside the country's border. Inside the country's border, distinction between the area with flight lines and without flight lines is visible.
The area without the flight lines has slightly higher standard deviation than the area with the flight lines, indicating missing data.
After obtaining the geoid residuals from LSC and its s error estimates, restoring of geoid height was done by ,

< A c c e p t e d M a n u s c r i p t >
where is the final geoid value obtained after adding the computed from the XGM2019e to the geoid residuals . The final geoid is presented in Figure 6. Figure 6 Regional geoid for Nepal using LSC (NPG20) employing airborne gravity data as input data and XGM2019e as the reduction model Nepal Geoid (NPG20).

Comparison with previously computed geoid model (NPG11)
For sake of simplicity, the output geoid of this study is named as NPG20, and it is compared with Nepal Geoid 2011 named as NPG11. The statistics of NPG11 are obtained from Forsberg et al, 2014. The geoid NPG11 was computed using the same airborne gravity data that has been used for this study along with 1114 surface gravity points (in a previous analysis, we found very large inconsistencies of these surface gravity with the global reference models, therefore, in this study we omitted the ground-based data and the global potential model with airborne gravity data was only used. The global potential model used was EGM2008. Downward continuation of gravity data was done using Least Squares Collocation and final gravimetric geoid was computed using spherical Fast Fourier Transform method. The comparison between outputs obtained from both geoid models is presented in   NPG11 (Nepal Geoid 2011) andNPG20 (Nepal Geoid 2020) At first glance, the values presented in Table 2 show that the geoid residuals of NPG11 and NPG20 seem to be very similar regarding their standard deviation. The use of 1114 surface gravity points in the determination of NPG11 do not seem to have greater impact in the geoid solution, because the airborne gravity data seem to dominate the surface gravity data in this solution. However, the range of restored geoid values (min&max) suggests that, there exists large differences of around -11 meters between these geoid models, which are depicted in Figure 7. In most parts of country, there is a bias of -3 m. In the north western region of the country (black ellipse), the differences range up to -20 m to 5 m. One possibility for such differences would be due to the presence of data in input gravity anomaly. There is a data gap in the same area where the geoid differences between NPG11 and NPG20 are higher (inside black box). However, the small deviation (6-8 cm) in estimation/interpolation from collocation procedure ( Figure 5) in missing flight lines do not suggest to be the reason for such high frequency values (-20 m) seen in Figure 7. The input airborne gravity dataset itself is not sufficient to account the full gravity signal. This is also indicated by the spectral behaviour of geoid residuals in Figure 4.

< A c c e p t e d M a n u s c r i p t >
Since the high-frequency differences were not likely to be the result of input data inconsistencies or data interpolation, i.e. from LSC method, the GGMs used for computation of geoid models (NPG11 & NPG20) was compared. The difference in the two GGMs is presented in Figure 8. Figure 8 clearly shows the significant difference between the reference models used for geoid computation.
There are differences between EGM2008 and XGM2019e in the data sources they are built on, and in the combination method. The discrepancies between geoid undulations computed from EGM2008 and those computed from independent GNSS/leveling data over areas covered with high quality gravity data (e.g., USA, Europe, Australia) are on the order of ±5 to ±10 cm (Pavlis et.al,

< A c c e p t e d M a n u s c r i p t >
Page 12 in 18 pages 2012). In contrast to this, XGM2019e provides slightly improved behaviour in the magnitude of a few mm RMS over land up to degree 720 (Zingerle et al., 2020).
NPG11 was computed by restoring EGM2008 up to degree 720 (Forsberg, et al., 2011). The difference between geoid values of NPG11 and the geoid values computed from EGM2008 up to degree 720 was computed and presented in Figure 9 (a). The difference in geoid values of NPG11 and the geoid values computed from XGM2019e up to degree 720 is presented in Figure 9 (b).

Figure 9 Geoid residuals computed from the differences of NPG11 from EGM2008 (a) and XGM2019e (b) up to d/o 720 respectively
The geoid residuals obtained using EGM2008 up to 720 presented in figure 9 (a) are higher than the geoid residuals obtained using XGM2019e up to 720 presented in figure 9 (b). The high frequency signal from topography in Figure 9 (a) shows up clearly in the same area where there exist large differences between the models as seen in figures 7 and 8. The long-wavelength structure dominates the difference between NPG11 and EGM2008. In contrast, the difference between NPG11 and XGM2019e as presented in Fig. 9 (b) shows shorter wavelength differences. Accordingly, XGM2019e fits much better to the airborne gravity observations than EGM2008, and this improvement in the background model can be considered as main advantage of the new NPG20 model. This shows that model EGM2008 is not sufficient to model the high frequency geoid signals in the highly mountainous areas.
On the other hand, the geoid residuals obtained from XGM2019e are smaller and smoother. From this comparison in Figure 9, one can say that high frequency signals coming from topography are modelled more suitably by XGM2019e than EGM2008. Also, NPG20 seems to be more consistent with the latest GGMs in the long-wavelength part due to the use of XGM2019e instead of EGM2008 as a background model.

< A c c e p t e d M a n u s c r i p t >
Page 13 in 18 pages

CONCLUSION
This paper focuses on the determination of a new geoid model for Nepal using only airborne gravity data by means of LSC, and the comparison with a previously computed geoid model.
The overall methodology of this study includes: (1) selecting an optimal GGM as the reference gravity field; (2) assessing the quality and accuracy of the airborne gravity data; (3) reducing the gravity data through the use of optimal GGM (XGM2019e); (4) determining of empirical covariance function from the reduced data; (5) fitting a model covariance function using a global covariance model; (6) applying Least Squares Collocation approach to interpolate, estimate and downward continue gravity anomaly data to determine geoid residuals; (7) restoring the geoid height computed from XGM2019e to the geoid residuals to obtain final geoid NPG20.
This study has provided the opportunity to deal with various GGMs and compare them for the best fit to the input gravity data. XGM2019e not only provided a better fit to the input data, but has also reduced one step in the reduction and restoring of gravity data, as it includes the topographyinduced gravity reduction as well. The application of the LSC method has provided an easy, efficient mathematically compliant way for determining geoid for Nepal.
The comparison of geoid models, NPG11 and NPG20 in Figure 7, shows that there exist significant differences in these geoid models. These differences are neither related to collocation procedure nor to the distribution of input data. The airborne gravity data themselves are not dense enough to capture such high frequency signals as they are revealed by Figure 7. Figure 8 suggests that the differences between these geoid models are caused by the difference in the GGM that has been used to determine these models. Figure 10 shows that the major differences of these models are in regions with elevations above 7000 m. These regions are Annapurna, Manaslu, Everest, Dhorpatan, and Kanchenjunga massifs. (these massifs are listed as the protected sites by Nepal Government).
Figure 10: Geoid effects due to Topography. The white triangular points are the peaks that are above 7000m plotted over geoid differences (NPG11 and NPG20) Figure 7 suggests that the input airborne gravity dataset itself is not sufficient to account the full gravity signal and the high-frequency differences between NPG20 and NPG11 are almost solely related to the different handling of the high-frequency gravity field component of the background model. The geoid residuals obtained using EGM2008 presented in figure 9 (a) are higher than the geoid residuals obtained using XGM2019e presented in figure 9 (b). In case of EGM2008, the highfrequency component suggests the presence of omission error, as it is neither represented by the background model nor the available input data. On the other hand, Figure 9 (b) shows that these high frequency components are well included in NPG11 via forward-modelled gravity in the frame of the XGM2019e model. Assuming that the forward-modelled high-frequency gravity signal is closer to reality than completely neglecting the high-frequency geoid component, NPG20 seems to be the better geoid model.
From this study, it can be concluded that the use of currently available airborne gravity data only is not sufficient to precisely determine geoid for Nepal, because their distribution is not dense enough to fully represent the high frequency gravity signal in the mountainous regions. The most straightforward way to capture these high-frequency signals could be to have denser sampled airborne data, ideally even in lower flight altitudes. High-quality surface gravity observations for filling the gaps of airborne data would be helpful in determining spectrally complete geoid signal.
Due to unavailability of reliable GNSS-levelling points, the external validation of NPG20 was not done. The establishment of reliable GNSS-levelling points would be helpful in determining accurate corrective surface for fitting geoid model and properly encounter the accuracy of geoid model, as well as to resolve the reference level issue that is indicated by Figure 7.