The 2021 gravimetric and hybrid geoid models for Taiwan and offshore islands: application to a seamless vertical datum

Wenhsuan Huang, Ching-Chung Cheng, Cheinway Hwang, Daocheng Yu, Te Wei Chen and Yen-Ti Chen Dept of Civil Engineering, National Yang Ming Chiao Tung University, Hsinchu, Taiwan, ROC This manuscript is submitted to “Terrestrial, Atmospheric and Oceanic Sciences” Running title: The 2021 Taiwan geoid models for a seamless vertical datum Key Points 1. A new hybrid geoid for Taiwan and its islands for elevations over land and seas 2. Latest gravity and elevation data for geoid modelling by the condensation method 3. Heights of objects in Taiwan Strait from the geoid are in the TWVD2001 vertical system * Correspondence to: Cheinway Hwang (cheinway@nctu.edu.tw; cheinway@gmail.com)


Introduction
A geoid model is one of the most important infrastructures in modern geodetic surveys. It is needed for converting GNSS-derived ellipsoidal heights to orthometric heights. This modern method of orthometric heighting can reduce the cost of mapping projects and contribute to economic development. Many countries have invested significant resources on constructing high-precision and high-resolution geoid models. For instance, in the United States, the National Geodetic Survey (NGS) has released several gravimetric-only geoid and hybrid geoid models since 1990 (https://geodesy.noaa.gov/GEOID/). Recent US gravimetric-only geoid models have been enhanced by new gravity data from the Grav-D program over 2008-2022. The latest gravimetric-only geoid model (USGG2012) was released in 2012, while the latest hybrid geoid model (GEOID18) was released in 2018 (https://geodesy.noaa.gov/GEOID/models.shtml). In 2014, Taiwan released a geoid model (Hwang et al., 2020), which was derived from many gravimetric datasets, including those from terrestrial, airborne, shipborne and satellite altimetry measurements. Many scholars and institutions around the world have been developing geoid models. For example, the Chinese Academy of Surveying and Mapping has released various geoid models for defining normal heights in China. In one study, Li et al. (2015) determined a quasi-geoid model for the mainland China and its vicinity areas based on the Molodensky theory. Huang and Veronneau (2013)

developed the 2013 Canadian
Gravimetric Geoid (CGG2013), which is expected to be updated in 2024 (https://reurl.cc/3NQE19). Miyahara et al. (2014) presented the latest geoid model for Japan, GSIGE2011. Gatchalian et al. (2016) upgraded the Philippine Geoid Model 2014 (PGM2014) to PGM2016 using re-processed and densified land gravity data. Saadat et al. (2018) and Ramouz et al. (2019) separately determined a regional geoid model of Iran (IRG2016; IRG2018) based on radial basis functions and the least squares collocation method. Featherstone et al. (2018) have been continuously improving the geoid model for Australia, with the latest geoid model being released in 2018 and it contains the error estimates for the geoidal heights in the model. Yildiz et al. (2021) released a hybrid geoid model for Turkey, called Turkish Geoid Model-2020 (TG-20), computed by using a modified Stokes' kernel with an additive correction method. Işık et al. (2021) published the latest geoid model for the US state of Colorado using a least-squares modification of Stokes' and Hotine integral formulations. Furthermore, there have been a number of geoid projects in the EU (Denker et al., 2009), Africa, Southeast Asia, and South America (see (http://www.isgeoid.polimi.it)).
A land vertical datum is the zero surface from which an elevation (typically orthometric height) is defined. On the main island of Taiwan, the zero surface for land corresponds to a vertical datum called TWVD2001 (Yang et al., 2003) and is at the mean sea level of Keelung in northern Taiwan based on 50 years of tide gauge records. In 2014 and for height modernization, the surface of TWVD2001 was released and specifically defined by a hybrid geoid on a 30"× 30" grid containing gridded geoidal (ellipsoidal) heights of the hybrid geoid covering only for the island of Taiwan (Hwang et al., 2020). On the offshore islands of Taiwan and mainland China, the earlier zero surfaces (before 2010) were the mean sea levels defined by short (few years) tide gauge records at the key tide < A c c e p t e d M a n u s c r i p t > gauge stations of these islands, but were updated in the 2010s using longer and better tidal records.
On the other hand, the zero surface for ocean depths is normally a surface defined by the lowest astronomical tide (LAT), which is used and needed for marine navigation safety. The difference between the zero surfaces for land elevations and ocean depths poses a significant challenge for managing marine resources. As such, a seamless vertical datum has been proposed for managing spatial information over a region covering both land and seas (Zhao et al., 2006;Ziebart et al., 2007;Slobbe et al., 2014). A seamless vertical datum consists of a set of models for the geoid, the LAT surface, the mean sea surface (MSS) and many other needed surfaces around a given region (see also Section 5). In practice, all such models contain ellipsoidal heights on a regular grid with the same ellipsoidal parameters (semi-major axis and flattening). If the ellipsoidal height of an object is measured (typically by GNSS), the object's vertical distance to any of the surfaces in the seamless vertical datum can be determined. Obtaining the vertical distance of the object to a new surface requires only the difference between the ellipsoidal heights of the original surface and the new surface at the position of the object.
In view of the unique role of a geoid model in a seamless vertical datum for Taiwan, and the newly collected gravity data on some of the offshore islands of Taiwan and mainland China, the objective of this study is to (1) update the geoid model on the island of Taiwan, and (2) extend the geoid model to the Taiwan Strait and offshore islands. The offshore islands in this paper include Kinmen, Matzu, Penghu, Liuqiu, Lyudao and Lanyu (from west to east), where we have collected new GNSS/leveling data to validate our new gravimetric and hybrid geoid models. We will also use ocean depths and land elevations around Taiwan and the offshore islands to compute the terrain corrections needed for Faye gravity anomalies. All the GNSS/leveling data are used to construct a hybrid geoid that covers the listed islands and the Taiwan Strait. The role of the 2021 hybrid geoid model in a future seamless vertical datum of Taiwan will be discussed, with a special emphasis on orthometric height determination in offshore wind farms.

New gravity anomalies for offshore islands
The current official gravimetric and hybrid geoid models of Taiwan were released in 2014.
The two models have been used for height modernization, cross-island height datum connection and LiDAR mapping of orthometric heights (Hwang et al., 2020). The two models cover the island of Taiwan and were computed using land-based and airborne gravity observations on the island of Taiwan, which are used for the new geoid models presented in this paper. In order to extend the coverage of the geoid to Kinmen and Matzu (offshore mainland China), Penghu (in the Taiwan Strait), Liuqiu (offshore Taiwan), and Lyudao and Lanyu (in the Pacific Ocean east of Taiwan), we obtained gravity data over these offshore islands. Figure 1 shows the distributions of the gravity sites where relative gravity measurements were collected. These gravity sites are on the first-order and secondorder leveling benchmarks (Table 1). Except Liuqiu, an absolute gravity site was established for each of the islands. For each island, the relative gravity measurements were adjusted using the program "gravnet" (Hwang et al., 2003) by holding fixed the gravity value at the absolute gravity site. Table   1 shows the information about the point gravity measurements at the offshore islands.  Table 2 shows the changes in gravity anomalies (a constant for each island) due to the changes in height datums.

Gravity anomalies in offshore waters
The marine gravity anomalies for our geoid modelling are from the latest global gravity grid of Sandwell et al. (2019), constructed from the sea surface height measurements of the Geosat, ERS-1, Jason-1, Jason-2, Cryosat-2 and ALtiKa altimeters (Figure 2a). We compared the altimeter-derived gravity anomalies (Sandwell et al., 2019) with the gravity anomalies from shipborne gravity surveys made between 2006 and 2011 (Hwang et al., 2014). Figure 2b shows the differences between the gravity anomalies from Sandwell et al. (2019) and those from shipborne measurements around the coast of Taiwan. The root-mean-squared difference is 3.58 mgal. In addition to the land and altimeter-< A c c e p t e d M a n u s c r i p t > derived gravity anomalies, we also used gravity anomalies from three airborne surveys in Taiwan (Hwang et al., 2020). Figure 2c shows the gridded gravity anomalies from all possible gravity data sources, which were used for the FFT computation of the new geoid models in this paper (Section 3).
The gravity field in Figure 2c is largely correlated with the terrains of land and the ocean topography (see Section 2.3).

Digital elevation model and depth model
A digital elevation model (DEM) is needed for computing terrain corrections in this paper (see Section 3). This paper uses the same DEM as that used for the island of Taiwan (Hsiao and Hwang, 2010;Hwang et al., 2020) and adds DEMs for the offshore islands as follows.  (Table 2). These offshore island DEMs were constructed from data based on several years of photogrammetric surveys. The original DEM database was in the .tif image format, which was converted to the .txt format using QGIS. Their spatial resolution is 20 m.
(b) The DEM for Matzu is from the satellite mission TanDEM-X, which has a spatial resolution of 50 m (Figure 3f; Martone et al., 2016;. There is no photogrammetry-based DEM for Matzu. Our terrain correction computation also requires ocean depths. In this paper, the needed ocean depths are from the Ocean Data Bank, Ministry of Science and Technology, Taiwan and the depths are given on a 500-m grid. Figure 4 shows the land topography and ocean topography from these elevation and depth data around Taiwan. The terrains over Penghu and Matzu are smooth with the largest elevations below 100 m. The terrain over Kinmen is somewhat rugged in a small part of this island. The terrains of Liuqiu, Lanyu and Lyudao are typical of a volcanic island and can result in rapid geoidal variations over a short distance. The RMS differences between the elevations from leveling and from the DEMs at the offshore islands (except Matzu) range from 0.473 m to 0.862 m.

Observed geoidal heights at leveling benchmarks
An observed geoidal height at a leveling benchmark is the difference between a GNSS-derived ellipsoidal height and the orthometric height defined in a local vertical datum. Geoidal heights can be used to assess the accuracy of a geoid model and can be used to create a hybrid geoid model from a gravimetric geoid model. On the island of Taiwan, we used the same observed geoidal heights as those used in the 2014 geoid models (Hwang et al., 2020). On the offshore islands, we obtained new observed geoidal heights by differencing the ellipsoidal heights and orthometric heights at the first and second-order benchmarks (see Table 1). As stated in Section 2.21, the orthometric height at a benchmark in an offshore island is the vertical distance from the benchmark to the zero surface defined by the island's mean sea level. As shown in Table 2, in 2019 the islands of Matzu, Penghu and Lanyu underwent a vertical datum change. When determining the observed geoidal heights, we have considered such datum changes for these three islands.
In Taiwan, the first-order leveling requires a 2.5 mm√ misclosure for a double-run leveling session, where k is the distance (in km) between two neighbouring benchmarks. The leveling of the offshore islands in this paper was carried out using a Dini 12 precision electronic level. The differential heights from the precision leveling were corrected for the effect of orthometric correction and then adjusted by holding fixed the height at a benchmark near the tide gauge that defines the local mean sea level. As such, the typical standard error of an orthometric height is few mm. In addition, the session lengths for the GNSS data vary from one station to another. A typical formal error in the ellipsoidal heights is about one cm or less. Figure 5 shows the observed geoidal heights over Taiwan and the offshore islands. The observed geoidal heights over the high mountain region of Taiwan ( Figure 5) are distributed only along the major cross-island roads, but they provide a crucial contribution to the geoid accuracy assessment and the creation of the hybrid geoid.

Method and steps for constructing the 2021 gravimetric and hybrid geoid
The method for constructing the gravimetric and hybrid geoid is the method that was presented in Hwang et al. (2020). The method is based on Helmert's condensation theory and uses the procedure of remove-compute-restore (RCR). A summary of the computational steps and related formulae are given below.
Step 1: Compute the terrain corrections for the gridded free-air anomalies (Figure 2c) where ′ , are the elevations at the computational point and at the contributing point (both from the DEMs in Section 2.3, including depths), 0 is the horizontal distance between the points with ′ and , G is the gravitational constant, ρ is rock bulk density (2.67 g/cm3), dσ is differential spherical surface area and R is the mean earth radius (about 6371 km).
Step 2: Compute the residual gravity anomaly where ∆g is full gravity anomaly and given in Figure 2c, and is the gravity anomaly from the EGM2008 model to harmonic degree 2160 (Pavlis et al., 2012).
Step 3: Compute the residual Faye gravity anomaly = ∆ + (4) Step 4: Compute the residual height anomaly where is the spherical distance is latitude, is the modified Stokes' kernel (Wong and Gore, 1969), with M =180 and 1 and 1 −1 are the direct and inverse one-dimensional fast Fourier transform each operating over a latitudinal band. We have experimented with different M values for in Eq. (5). The use of M =180 leads to the best agreement between the modelled gravimetric geoidal heights and the observed geoidal heights (Section 4; Table 3 and Fig. 7).
Step 5: Compute the quasi-geoidal heights corrected for the indirect effect where  is height anomaly from the EGM2008 model to degree 2160 and the last term is the indirect effect, which is a function of the elevation (H) and normal gravity .
Step 7: Use the observed geoidal heights to construct a hybrid geoid In this step, the differences between the observed geoidal heights ( Figure 5) and the gravimetric geoidal heights at the leveling benchmarks are first computed. Then, the differences are interpolated onto a 30" × 30" grid using the "surface" command of GMT (Wessel et al., 2019) with a tension factor of 0.25. This difference grid is added to the gravimetric geoid grid to form the hybrid geoid grid. Figure 6a shows the gravimetric geoid created in Step 6 and Figure 6b shows the hybrid geoid created in Step 7. Both models cover 118˚E-125˚E and 21˚N-27˚N and are on a 30"×30" grid. On a visual inspection, the two models are similar in pattern and magnitude. However, the two models can differ by few dm in several parts of Taiwan due to the use of the reference geoidal surface (EGM2008), varying gravity data densities and a number of other reasons. Section 4 will show a quality assessment of these two geoid models, and a comparison between the geoidal heights and selected altimeterderived mean sea surface heights in the Taiwan Strait.

Assessment on land
For the assessments of the two geoid models on the island of Taiwan, we used the observed geoidal heights along 14 leveling routes (Figure 7). The expected accuracies of these observed geoidal heights are at the cm level because the sessions of the related GNSS observations are 12 to 24 hours (Hwang et al., 2020). On the six offshore islands, we inspected the documents of GNSS measurements and leveling data to choose the best observed geoidal heights for the assessment. The numbers of the chosen observed geoidal heights are as follows (Figure 7): 7 for Kinmen, 3 for Lanyu, 12 for Penghu, 3 for Liuqiu, 3 for Lyudao and 4 for Matzu. First, the difference between the observed ( ) and the modeled ( ) geoidal height at the ith leveling benchmark is: The mean, standard deviation root-mean-squared (RMS) value of the differences for n benchmarks are:

< A c c e p t e d M a n u s c r i p t >
The statistics of the differences between the observed geoidal heights and the model values from the gravimetric geoid and the hybrid geoid are shown in Table 3a (Taiwan) and 3b (offshore islands) and Table 4a (Taiwan) and 4b (offshore islands). Table 3a (Taiwan) shows that, the mean differences vary from 2.92 cm to 28.86 cm along the 14 leveling routes, despite just few cm of standard deviations in the differences between the geoidal heights from the observations and from the gravimetric model. On the six offshore islands (Table 3b), the magnitudes of the mean differences and standard deviations are similar to those on Taiwan. This suggests that the gravity data density, terrain steepness and other factors can introduce biases in the gravimetric geoid, which cannot be directly used for GNSS determination of orthometric heights. On the other hand, Table 4a and 4b show that the means and the standard deviations of the differences between the geoidal heights from the observations and from the hybrid geoid are only few cm, suggesting that the blending of the gravimetric geoid with the observed geoidal heights has resulted in a hybrid geoid suitable for direct orthometric heighting at the cm-level on the island of Taiwan and the six offshore islands.

Assessment the hybrid geoid in the Taiwan Strait
The hybrid geoid will be used for determining the orthometric heights of infrastructures in the Taiwan Strait under the TWVD2001 system, thus its accuracy is of a great concern for users such as builders of offshore wind farms in the Taiwan Strait. Unlike a geoid model on land, there are no observed geoidal heights at sea for assessing the accuracy of a given geoid model. However, a Baltic Sea geoid model (Nordman et al., 2018) has been evaluated by shipborne GPS measurements. In this case, the dynamic oceanic topography (DOT) was removed from a Baltic Sea DOT model from the ocean-tide-corrected sea surface heights from GPS to determine geoidal heights at the Baltic Sea, which were compared with the model geoidal heights. In the Taiwan Strait, we carried out two assessments for the 2021 hybrid geoid as follows.
In the first assessment, we determined the differences between the ellipsoidal heights from the geoid models and from the DTU 2018 mean sea surface model (DTU18) (Andersen et al., 2018), as shown in Figure 8. The differences are mostly between 0 to 20 cm in the Taiwan Strait. Off the coast of southeast mainland China, a region centered at 24.8°N and 119.1°E contains relatively large differences up to 50 cm. East of Taiwan at 24°N (coastal Hualien County), there are large negative differences (up to -20 cm). The large positive differences east of 122°E probably result from the oceanic dynamic topography generated by the north Pacific oceanic gyre and altimeter data errors. In a spot off southwest Taiwan and another spot off northwest Taiwan, there are also relatively large differences. The differences in Figure 8 suggest that, if the DTU18 model is to be the starting model for a best mean sea surface model around Taiwan, there is a need to investigate the causes of the large difference in Figure 8.
In the second assessment, we compared the along-track mean sea surface heights (MSSHs) from the TOPEX-Jason altimeters  with the geoidal heights from the 2021 gravimetric and hybrid geoid models and the geoid models released in 2014 (Hwang et al., 2020), as shown in Figure 9a (the various ellipsoidal heights) and 9b (the differences between ellipsoidal heights). Pass 164 (a descending track) crosses the Taiwan Strait from the east coast of Fujian Province of mainland China to offshore Taichung. Because the coastal observed geoidal heights on Taiwan and its offshore islands were used for constructing the hybrid geoid and because the fact that the vertical datums are the mean sea levels at the respective main tide gauge stations on Taiwan and offshore islands, the geoidal heights (specifically ellipsoidal heights) from the hybrid geoid should be consistent with the mean sea surface heights around Taiwan and the offshore islands. Figure 9a shows that, along Pass 164 the ellipsoidal height differences in all cases increase with latitude, suggesting increasing inconsistences between the sea surface heights from DTU18 and the model values from the four geoid models (2021 gravimetric geoid, 2014 gravimetric geoid, 2021 hybrid geoid and 2014 hybrid geoid).
The differences (blue line in Figure 9b) between DTU18 and the 2021 hybrid geoid are the least, particularly near the west coast of Taiwan. The comparisons in Figure 9a and 9b point out potential improvements for a global mean sea surface model like DTU18, which has benefited the constructions of seamless vertical datums world-wide (Andersen et al., 2018; see also Section 5).

The hybrid geoid for a seamless vertical datum and heighting in offshore wind farms
The hybrid geoid developed in this paper is part of a seamless vertical datum that provides various reference surfaces for elevations on land and depths at sea (Ziebart et al., 2007;Slobbe et al., 2014). Figure10 shows the various surfaces that constitute such a seamless datum over land and seas. Numerically, any surface in Figure 10 is represented by ellipsoidal heights at regular intervals along the latitude and longitudinal directions. For example, the grid intervals for the hybrid geoid in this paper are 30". As explained in Section 4.2, the hybrid geoid is also applicable to the waters in the Taiwan Strait because it incorporates the observed geoidal heights on the island of Taiwan and the six offshore islands.
It can be difficult to construct a numerical mean sea surface model that provides accurate ellipsoidal heights near coastal areas (see the mean sea level in Figure 10). This is because the largest source of measurements for a coastal mean surface is satellite altimeters, which normally have poor height measurement accuracies near coasts due to waveform contaminations and large footprints (at few hundreds of meters or larger) that allow only low-resolution coastal height observations (see < A c c e p t e d M a n u s c r i p t > Figure 9). Such poor radar altimeter measurements are worsened by low-accuracy corrections for the effects of ocean tides and wet tropospheric delays.
The surface of the lowest astronomical tide (LAT) in Figure 10 is commonly used as the chart datum, which is the zero surface of ocean depths. A numerical chart-datum model is formed by a numerical mean sea surface model and a numerical model of gridded lowest tides, whose accuracies can be affected by the low-accuracy data for the model constructions. The depth of seafloor is the difference between the ellipsoidal height of the chart-datum model and the ellipsoidal height of the seafloor; the latter can be obtained by differencing the GNSS-determined ellipsoidal height of the ship and the range from the ship to the seafloor. Figure 10 also shows a sample picture of a GNSS RTK base station, where a continuous GNSS receiver is installed. The broadcasting device at this base station sends GNSS signals to a receiver on a ship, helping to guide the machine on the ship to install a wind turbine. In Taiwan, the required vertical accuracy in the machine guidance is at the cm level in the TWVD2001 system (Yang et al, 2003).
In Taiwan, it is required that the orthometric height of a man-made structure such as the RTK base station in Figure 10 over the waters around Taiwan is based on the same vertical datum as that used on land, i.e., TWVD2001. If a target at sea is visible from a leveling benchmark, its TWVD2001based orthometric height can be measured by trigonometric leveling from this benchmark, but the error of this height can easily exceed 10 cm. If a target is not visible from land, its TWVD2001-based orthometric height can be obtained by subtracting the geoidal height (from the hybrid geoid) from a GNSS-measured ellipsoidal height. If both the accuracies of the geoidal height and the GNSSmeasured ellipsoidal heights are at the cm level, the resulting orthometric height will have a cm-level accuracy.
In Taiwan, wind is an emerging source of green energy. In 2020, the power generated by wind is 2289.3 Giga-watt-hour (GWh; https://en.wikipedia.org/wiki/Wind_power_in_Taiwan). Figure 11a shows the wind farms under construction and under plan, and four GNSS RTK stations in offshore Taiwan and in the Taiwan Strait. Figure 11b shows some of the installed wind turbines off the coast of Hsinchu County and Miaoli County. The horizontal and vertical coordinates of all wind turbines, RTK base stations and future man-made structures related to the wind farms in Figure 11 should be given in the horizontal coordinate frame of Taiwan (TWD1997) and in the vertical datum TWVD2001. The hybrid geoid is the unique model that can be used for determining TWVD2001based orthometric heights using GNSS devices (see also Figure 10).

Conclusion
The 2021 Taiwan geoid models on a 30"×30" grid have been developed in this paper to update the previous 2014 Taiwan geoid models. The 2021 models are based on the existing gravity data used for the 2014 geoid models on the island of Taiwan, and the new gravity data on the offshore islands Kinmen, Matzu, Penghu, Liuqiu, Lyudao and Lanyu. A notable example of re-computation of free-air gravity anomalies using a new vertical datum is given in this paper using the gravity data and vertical datum changes on the islands of Penghu, Matzu and Lanyu. The vertical datum updates result in 0.1-0.2 mgal of gravity changes. However, such gravity changes are inconsequential for a better match between the geoidal heights from the gravimetric geoid and the observations. Still, a hybrid geoid from a combination of the gravimetric geoid and sufficiently dense observed geoidal heights is the optimal model that can be directly used to determine an orthometric height from GNSS at a new point. Such an orthometric height is directly in the vertical datum that defines the orthometric heights at the points where the observed geoidal heights have contributed to the hybrid geoid. For example, the 2021 hybrid geoid in this paper can be used to determine orthometric heights in the TWVD2001 vertical datum of Taiwan and in the updated vertical datums of the six offshore islands.
In addition to the assessments on land, the 2021 hybrid geoid is also assessed in the Taiwan Strait where we compare the geoidal (ellipsoidal) heights from the hybrid geoid and those from the DTU18 mean sea surface model and TOPEX-Jason altimeters. The assessment result indicates increasing discrepancies between these three sets of ellipsoidal heights from the west coast of Taiwan toward the southeast coast of mainland China. The 2021 hybrid geoid has blended the newly observed mean sea levels (via the observed geoidal heights) in Taiwan and the six offshore islands, so it should result in model sea surface heights that conform to the true sea surface around these islands. Thus, the assessment result in the Taiwan Strait may provide a clue for advancing a global model like DTU18 to benefit the construction of a seamless vertical datum over regions covering both land and sea.