A detailed quasigeoid model of the Hong Kong territories computed by applying a finite-element method of solving the oblique derivative boundary-value problem

1 Introduction

In the Hong Kong territories, the geodetic datum currently used for height control is the Hong Kong Principal Datum (HKPD). The heights of HKPD levelling benchmarks were determined by adjusting levelling measurements but without involving gravity observations along levelling lines. To eliminate systematic errors due to disregarding the gravity information, Nsiah Ababio and Tenzer (2022a) used terrestrial and marine gravity measurements to interpolate gravity values along levelling lines and then computed the normal and orthometric corrections to levelled height differences to determine the normal and orthometric heights of HKPD levelling benchmarks. The result was presented as the Vertical Control Network 2022 (VCN2022). Since the current coverage of levelling network throughout the Hong Kong territories is not suitable for many surveying applications, Nsiah Ababio and Tenzer (2022b) compiled the detailed gravimetric geoid model HKGEOID-2022, with the primary purpose of converting the geodetic (ellipsoidal) heights measured by the Global Navigation Satellite System (GNSS) to the orthometric heights defined in the updated HKPD. Nevertheless, both products (i.e. the VCN2022 and the HKGEOID-2022) do not meet the requirements for many proposed or approved engineering projects to be realized between Hong Kong and mainland China, because the normal heights and the quasigeoid model are used for a practical realization of vertical geodetic control in mainland China, whereas the orthometric heights and geoid models are used in Hong Kong. To achieve consistency between both vertical data, a new detailed quasigeoid model was prepared, and normal heights of levelling benchmarks were determined at the Hong Kong territories. As stated previously, the normal heights of HKPD levelling benchmarks were already determined by Nsiah Ababio and Tenzerr (2022a). In this study, we present a new detailed quasigeoid model compiled according to the method developed by Minarechová et al. (2021). Their method is based on applying the finite-element method (FEM) to solve the oblique derivative boundary condition (BC).

Traditionally, gravity anomalies have been used to determine geoid/quasigeoid models by means of solving the geodetic boundary-value problems (BVPs). The gravity anomaly is defined as the actual gravity (typically measured on or above the topographic surface) and the normal gravity computed as a function of the geodetic latitude and the physical (normal or orthometric) height (above the sea level). Nowadays, GNSS techniques provide information about the vertical position of gravity points with respect to the reference ellipsoid in terms of the geodetic (ellipsoidal) rather than physical heights. The gravity disturbance is then obtained as the difference between the observed and normal gravity (both referred to the same point at or above the topographic surface). The availability of this gravity data resulted in development of various methods that have been applied to determine geoid/quasigeoid models from gravity disturbances.

Gravity disturbances define the oblique derivative BCs of the fixed gravimetric boundary-value problem (FGBVP) that represents an exterior oblique derivative geodetic BVP for the Laplace equation (Koch and Pope 1972, Bjerhammar and Svensson 1983, Holota 1997). A detailed overview of various procedures for solving the oblique derivative BVP can be found, for instance, in Minarechová et al. (2021). Here, we only mention numerically efficient approaches, most notably the FEM, the finite-volume method (FVM), and the boundary element method (BEM). Principal reasons for using these methods in gravimetric geoid/quasigeoid modelling, instead of applying classical procedures based on solving Stokes and Poisson’s integrals, rely on the possibilities of applying a direct refinement of the discretization by using a very detailed gravity data grid and considering a real topographic relief.

Meissl (1981) and Shaofeng and Dingbo (1991) applied the FEM to solve gravimetric problems. Klees (1995) and Lehmann and Klees (1999) developed the indirect BEM approach that was later improved by Klees et al. (2001). Čunderlík et al. (2008) and Čunderlík and Mikula (2010) introduced the direct BEM approach. Later, Čunderlík et al. (2012) discussed the applicability of the BEM to solve the oblique derivative problem. Fašková et al. (2007, 2010), Šprlák et al. (2011), and Mráz et al. (2016) applied the FEM. Nevertheless, the oblique derivative BC was not considered in these studies. Fašková (2008) applied for the first time the FVM in gravity field modelling by solving the geodetic BVP with the Neumann BC. Macák et al. (2012) applied the FVM to the oblique derivative BVP, and Macák et al. (2014, 2015) applied it to solve FGBVPs on uniform grids. Medľa et al. (2018) presented the FVM for solving the oblique derivative BVP on 3-D unstructured meshes above the real topography. Droniou et al. (2019) further improved this method by treating the oblique derivative BC so that its tangential component is considered an advection along a topographic relief regularized by a carefully designed surface diffusion term. Minarechová et al. (2021) applied the FEM to solve the oblique derivative BC. In this method, the oblique derivative is incorporated directly into computational nodes by using two tangential vectors for each node in order to improve numerical accuracy.

In this study, we applied the method developed by Minarechová et al. (2021) to compile a detailed gravimetric quasigeoid model in the Hong Kong territories. Since there are some inconsistencies in definitions and realizations of vertical and gravimetric data (such as different datum origins and application of different tidal systems and reference parameters), the gravimetric quasigeoid model is usually (optimally) fitted or combined with the geometric quasigeoid heights determined at the GNSS-levelling benchmarks from their geodetic and normal heights. In the method applied here, this fitting is implicitly carried out by compiling a quasigeoid model. The method applied for gravimetric quasigeoid modelling is briefly reviewed in Section 2. Input data acquisition is discussed in Section 3. Results are presented and analysed in Section 4. The study is summarized and concluded in Section 5.

2 Theory

A formulation of the oblique derivative of the BVP is explained in this section. We then briefly discuss its solution by means of applying the FEM according to the numerical procedure developed by Minarechová et al. (2021).

3 Data acquisition

We used terrestrial and marine gravity data, GGMs, levelling data, GNSS geodetic heights, and topographic and bathymetric models to compile a new gravimetric quasigeoid model and combine it with GNSS-levelling data. This section provides a brief summary of datasets and models and their use for a discretization of the Earth’s surface and a preparation of input BCs, namely, the Dirichlet BC generated from GGMs and the gravity disturbances as the oblique derivative BCs considered on the discretized Earth’s topography.

3.1 Terrestrial and shipborne gravity data

From the regional gravity survey conducted by Electronic and Geophysical Services Ltd in 1991, 640 gravity observations were done in the Hong Kong territories (Electronic and Geophysical Services Ltd 1991). These observations comprise 133 shipborne measurements, and 507 terrestrial measurements on land spaced at approximately 2 and 2–4 km apart, respectively. Terrestrial sites were accessed using either vehicles or helicopters, and measurements were performed by the Lacoste and Romberg “G” relative gravimeter. The seaborne gravity measurements were performed using the marine model H/U seabed gravimeter. The distribution of gravity data is shown in Figure 1. The gravity database has an accuracy of ±0.03 mGal deduced from its connection to the International Gravity Standardization Net 1971 (Evans 1990).

Figure 1

Distribution of the terrestrial and shipborne gravity data.

3.2 Discretization of the computational domain

The Hong Kong territories, covering a total area of about 1,100 km², are characterized by large topographic elevation changes with the maximum height reaching 957 m (Tai Mo Shan). Our 3D computational domain for the FEM numerical scheme was considered on such complicated topography. To get a reliable FEM numerical solution in this territory, we had to consider an extended area, as shown in Figure 2. Its size is approximately two times bigger than the size of the inner-zone area with the input terrestrial and shipborne gravity data. In this way, we minimized the impact of the Dirichlet BC prescribed on the side boundaries of the FEM solution inside the inner zone, as discussed by Fašková et al. (2010). Finally, the chosen computation area was bounded by meridians of 113.5°E and 114.8°E longitudes and parallels of 21.8°N and 22.9°N latitudes and discretized with the high-resolution 0.0005 deg × 0.0005 deg topographical grid.

Figure 2

The discretized Earth’s topography in the Hong Kong territories and surrounding areas as the bottom boundary of the computational domain (the resolution 0.0005 deg × 0.0005 deg).

Within the Hong Kong territories, we used the HK_DTM_5m digital terrain model compiled by the Lands Department (www.landsd.gov.hk/en/spatial-data/open-data/kf_dtm.html). We also used the HKGEOID-2022 geoid model Nsiah Ababio and Tenzer (2022b) to obtain the geodetic (ellipsoidal) heights of computational nodes. In surrounding areas, we used the multi-error-removed improved-terrain (MERIT) digital elevation model (DEM) (Yamazaki et al. 2017) combined with the EGM2008 geoid model (Pavlis et al. 2012). The ellipsoidal heights of computational nodes offshore were interpolated from the DTU21 mean sea surface (Andersen et al. 2021). In this way, we constructed the bottom boundary of the 3D computational domain (Figure 2).

The upper boundary was chosen at the constant altitude of 200 km above the reference ellipsoid and was discretized with the same horizontal resolution of 0.0005 deg × 0.0005 deg. In the radial direction, the 3D computational domain was discretized non-uniformly depending on altitude. The radial size of finite elements on the Earth’s surface was set to 5 m while increasing linearly with altitude and exceeding 1 km for elements on the upper boundary. The whole 3D computational domain then consisted of 2,600 × 2,200 × 800 (longitude × latitude × height) = 4,576,000,000 finite elements (5,720,000 computational nodes on the discretized Earth’s topography).

3.3 Dirichlet BCs generated from GGMs

The Dirichlet BCs on the upper and side boundaries were generated from the GRACE/GOCE-based GGMs. In particular, we used the GO_CONS_GCF_2_DIR_R6 satellite-only GGM up to a degree/order of 300 (Bruinsma et al. 2014) to generate the disturbing potential on the upper boundary at an altitude of 200 km (Figure 3). The FEM solution was fixed to the gravity field information detected by the satellite missions CHAMP, GRACE, and GOCE. On the side boundaries, the disturbing potential was generated from the combined EIGEN-6C4 model up to a degree/order of 2160 (Förste et al. 2014).

Figure 3

Disturbing potential at an altitude of 200 km above the reference ellipsoid as the Dirichlet BC on the upper boundary generated from the GO_CONS_GCF_2_DIR_R6 model up to a degree/order of 300.

3.4 Gravity disturbances on the discretized Earth’s surface

On the bottom boundary, the input gravity disturbances were generated directly at computational nodes that discretize the Earth’s surface. Due to different sources of gravity data, we divided the bottom boundary into three zones, as shown in Figure 4. The Hong Kong territories covered by available terrestrial and shipborne gravity measurements represent the inner zone (“Zone_IN”). Here, we used the method of a reverse reconstruction of the gravity acceleration from the complete Bouguer anomaly (CBA) map. Outside this area, we generated input data from a combined GGM refined by the residual terrain model (RTM) methodology (“Zone_OUT_lands”). Over open sea (“Zone_OUT_sea”), we interpolated gravity disturbances from the DTU21_GRAV dataset, providing the altimetry-derived gravity data.

Figure 4

The bottom boundary divided into three zones with respect to gravity data sources.

3.4.1 Reconstruction of gravity from the CBA map

In the inner zone “Zone_IN” (Figure 4), we used the method of a reverse reconstruction of the gravity from the CBA map. The computation of the gravity value at an arbitrary point through a reverse reconstruction from the CBA map is considered to be the most accurate method (cf. Majkráková et al. 2016, or Zahorec et al. 2021). The reason is simple, an interpolation from the CBA map is more appropriate than, for example, from the free-air anomaly map. The first numerical step involved the computation of the CBA from the gravity measurements described in Section 3.1. We used a standard approach based on precise evaluation of the terrain/bathymetric corrections to a spherical distance of 166.7 km. For this purpose, we used TopoSK software (Zahorec et al. 2017) and the following DEMs: For the inner most zone (up to 250 m), we used HK_DTM_5m as the most detailed and available DEM of the Hong Kong territories. For the middle zone (250–5.2 km), we used HK_DTM_5m added by the Advanced Land Observing Satellite World 3D model version 2.1 (AW3D30) (Tadono et al. 2014; Takaku et al. 2018) and resampled to a resolution of 25 m. For the closer outer zone (5.2–28.8 km) and more distant outer zone (28.8–166.7 km), we used the MERIT DEM with resolutions of 3 and 30 arcsec, respectively. Bathymetric corrections were calculated using the SRTM15plus v.2.3 model (Tozer et al. 2019).

We used the CBA map for a correction density of 2,670 kg/m³ (Figure 5) to interpolate the CBA values on a regular grid of 50 m × 50 m. By the identical reverse calculation, including the evaluation of terrain/bathymetric corrections, we obtained values of the reconstructed gravity values at the FEM computational grid. From these values, we generated gravity disturbances at all grid points of the “Zone_IN” whose ellipsoidal heights had been reconstructed before (Section 3.2).

Figure 5

CBA map constructed from available gravity measurements (black dots).

3.4.2 GGM + RTM method

The input gravity disturbances at grid points of “Zone + OUT + lands” (Figure 4) were generated using the GGM + RTM methodology. We used the combined EIGEN-6C4 model up to a degree/order of 2,160. The RTM technique accounts for differences between the gravitational effect of the real terrain masses represented by high-resolution DEMs (the same used for CBA in Section 3.4.1) and smoothed mean elevation surface represented, for instance, by the DTM2006 model (Pavlis et al. 2007). The impact of RTM on the gravity generated from GGM is significant, as can be seen from the comparison with the gravity observed from terrestrial measurements in Figure 6.

Figure 6

Differences of the gravity from terrestrial/shipborne measurements and the gravity generated from GGM without RTM (black) and with RTM (red). Gravity differences are shown as a function of height differences between the detailed DEM and DTM2006.

3.4.3 Altimetry-derived gravity data

At grid points over the open sea in “Zone_OUT_sea” (Figure 4), the gravity disturbances were transformed from the altimetry-derived gravity anomalies interpolated from the DTU21_GRAV data set. For this transformation, we used the EGM2008 model (Pavlis et al. 2012). Due to the worse accuracy of the altimetry-derived gravity data in shallow waters and coastal areas, we did not use them in Lingding Bay (the Pearl River estuary in the northwest of Hong Kong) and Mirs Bay (in the northeast of Hong Kong).

Combining the gravity disturbances generated in all three zones, we finally obtained the input dataset on the bottom boundary as the oblique derivative BC for our FEM numerical solution of FGBVP. The spatial pattern of the gravity disturbances is depicted in Figure 7. We could recognize some minor discontinuities in data along borders between the selected three zones. Nevertheless, the magnitudes of these discontinuities are not that significant because we do not have at disposal any original gravimetric measurements outside the Hong Kong territories. Note that gravity data over mainland China are not publicly available.

Figure 7

Gravity disturbances at the Earth’s surface as the oblique derivative BC on the bottom boundary.

3.5 GNSS-levelling data

From the current vertical control network database in Hong Kong, there are only 16 benchmarks with high-quality GNSS measurements of geodetic (ellipsoidal) heights referred to the WGS84 reference ellipsoid at the ITRF96 epoch. From these 16 points, there are 10 first-order bedrock benchmarks and 6 second-order benchmarks. The normal heights of these benchmarks were obtained after applying the normal correction to levelled height differences and a subsequent readjustment of the whole levelling network that provided the VCN2022 solution. The new network adjustment and analysis revealed that the accuracy of VCN2022 normal heights is about ±1–2 cm, and a similar accuracy characterizes the GNSS vertical measurements Nsiah Ababio and Tenzer (2022b). The distribution of the GNSS-levelling benchmarks and their respective quasigeoidal heights are shown in Figure 8.

Figure 8

Distribution of GNSS-levelling benchmarks in the Hong Kong territories.

4 Results

A numerical solution of FGBVP based on using the FEM approach on a large 3D unstructured mesh with 4,576,000,000 finite elements (Section 3.2) required about 1.3 TB of internal memory. The large-scale parallel computations were performed on 176 cores of the cluster with the distributed memory and the NUMA (Non-Uniform Memory Access) architecture by using a hybrid parallelization (44 MPI processes, each with 4 Open MP threads). The BiCGSTAB linear solver converged after 5,970 iterations for a tolerance of 10⁻⁵. Hence, the parallel computations took about 35 h of CPU time.

A benefit of such large-scale computations is that they resulted in the disturbing potential obtained at every point of the entire 3D computational mesh, i.e. at all 4,576,000,000 finite elements. Hence, it is possible to derive different quantities of the local gravity field (e.g. the first, second, and higher derivatives in different directions) at every point. To get a local quasigeoid model, values of the disturbing potential T obtained on the bottom boundary (i.e. at the Earth’s surface; Figure 2) were converted to the quasigeoidal heights ξ by using the following formula:

where h is the geodetic (ellipsoidal) height, H ^N denotes the normal height, U is the normal gravity potential evaluated at the ith grid point at the Earth’s surface, γ is the mean normal gravity between the reference ellipsoid and the telluroid, and W ₀ is the global geopotential value adopted for a realization of the International Height Reference System (IHRS) (Sánchez et al. 2016). Parameters of the normal gravity filed were computed from the WGS-84 reference ellipsoid. In this way, the quasigeoidal heights have been expressed with respect to the WGS-84 reference ellipsoid and to the W ₀ value adopted for IHRS.

The local quasigeoid model over the Hong Kong territories (and surrounding areas) computed on a 0.0005 deg × 0.0005 deg grid is shown in Figure 9. To validate its precision, we compared the gravimetric quasigeoid heights with the corresponding geometric quasigeoid heights at GNSS-levelling benchmarks (Section 3.5). It is worth noting that the geometric quasigeoid heights were obtained from differences between the geodetic (ellipsoidal) heights (measured by the GNSS techniques) and the VCN2022 normal heights. In addition, we used 507 input terrestrial gravimetric measurements (Section 3.1) to validate the result of a gravimetric modelling. For this purpose, we converted the measured gravity anomaly data into the corresponding quasigeoid heights. The comparison of results is shown in Figure 10, with the statistics of differences summarized in Table 1.

Figure 9

Local quasigeoid model in the Hong Kong territories presented as the FEM numerical solution of FGBVP (with a resolution of 0.0005 deg × 0.0005 deg).

Figure 10

Validation of results: (a) Differences between the gravimetric and geometric quasigeoid heights at 16 GNSS-levelling benchmarks. (b) Differences between the quasigeoid heights obtained by solving FGBVP and computed from the measured gravity anomaly data at 507 gravity sites.

As seen from the comparison, there is a reasonable agreement between the quasigeoid models with both the geometric quasigeoid heights at GNSS-levelling points and with the quasigeoid heights computed from measured gravity data. The mean values of both differences (63.2 and 61.3 cm; Table 1) differ by less than 2 cm. This finding indicates good consistency between gravity, GNSS, and levelling data in Hong Kong. Moreover, this finding also assures that despite high-quality GNSS-levelling measurements are only at 16 benchmarks that are very irregularly distributed and clustered only at some locations, the validation of results should be quite realistic.

5 Summary and concluding remarks

We have compiled the first detailed quasigeoid model in the Hong Kong territories by applying the FEM to solve the geodetic BVP with oblique derivative BCs considered directly at computational nodes on the discretized Earth’s topography. Since the problem is defined for the oblique derivative BCs, we first converted the terrestrial and seaborne gravity measurements into the gravity disturbances at the Earth’s surface by utilizing information about the geodetic heights obtained from the digital elevation and geoid models (inland) and the mean sea surface topography (offshore). The absence of gravity data in mainland China was resolved by using the gravity information from a high-resolution GGM.

The result of gravimetric modelling revealed that the quasigeoid surface within the Hong Kong territories is mostly below the WGS84 reference ellipsoid, with a prevailing smooth pattern of decreasing quasigeoid heights in the northwest direction. Small irregularities in this prevailing pattern, mostly only a few millimeters, are mainly attributed to a rough topography in Hong Kong.

The validation of the gravimetric quasigeoid model revealed a relatively good agreement with the geometric quasigeoid heights at GNSS-levelling benchmarks, with the standard deviation of differences between the gravimetric and geometric quasigeoid heights of ±3.3 cm. Obviously, this accuracy specification was obtained from a minimal number of GNSS-levelling benchmarks distributed very irregularly over the territories, with most of them scattered only at a few locations. Moreover, these benchmarks are located at low elevations. Nevertheless, this value quite closely reflects the accuracy estimates of the geometric quasigeoid heights at GNSS-levelling benchmarks of about ±1–2 cm as well as expected uncertainties at the level of about ±2–3 cm that are attributed to the accuracy of gravity measurements (and a lack of gravity data in mainland China). This argument is supported by the fact that the standard deviation of differences between the quasigeoid heights obtained from gravimetric modelling and computed from the measured gravity data was found to be ±2.2 cm. In addition to these findings, a relatively good consistency was confirmed between the reference systems defining gravity and GNSS-levelling measurements.