The determination of Thailand Geoid Model 2017 (TGM2017) from airborne and terrestrial gravimetry

The Royal Thai Survey Department and Chiang Mai University developed the Thailand geoid model 2017 (TGM2017) with a 1’ × 1’ grid to support the transformation between Global Navigation satellite System (GNSS) ellipsoid heights and Kolak-1915 vertical datum orthometric heights. TGM2017 was based on Thailand gravimetric geoid model 2017 (THAI17G) and 299 GNSS ellipsoidal heights co-located with Kolak-1915 heights. All terrestrial gravity data used for geoid computation came from the new national gravity network, consisting of 87 absolute and 9929 relative gravity stations at 10 – 25 km intervals, mostly along with existing roads. From 2016 to 2017, airborne gravity surveys were conducted at a 4000m-flight altitude and 10 km along-track spacing to acquire the gravity data over mountainous and inaccessible areas, including coastal and marine areas, at an estimated accuracy of 3.0 mGal. Long-wavelength geoid structure was controlled by the GOCE-EGM2008 combined model (GECO) and the Technical University of Denmark’s global marine gravity model 2013 (DTU13). All gravity data were combined and downward, using least-squares collocation with the residual terrain model reductions from a digital terrain elevation data level 2 (DTED2). THAI17G was determined by multi-band spherical Fast Fourier Transform and converted to TGM2017 with the 38.2cm root-mean-square (rms) fit of 299 GNSS/leveling co-points and a mean offset of 37.0 cm. This value represents the separation between Kolak-1915 and a global mean sea level. The evaluation of TGM2017 at 100 GNSS/leveling checkpoints shows the rms of 4.9 cm, consequently leading to reliable orthometric heights at a 10-cm accuracy level or better.


INTRODUCTION
A height modernization system allows us to determine reliable and accurate orthometric heights using the Global Navigation Satellite System (GNSS) technology with an accurate geoid model.Such a modern height determination could be conducted with low cost and time, compared with spirit leveling.Accurate geoid models correspond to accurate orthometric heights for engineering works and natural disaster management, for instance, the risk assessment of landslides, the suitability of residential construction, and infrastructure rehabilitation and maintenance.Obtaining an orthometric height depends on an accurate geoid undulation that is determined from the required accuracy and resolution of gravimetric data on or near the Earth's surface.By these requirements, determining a geoid model for Thailand is difficult, mainly due to insufficient coverage and distribution of gravimetric data in the country and no possibility of accessing the data in neighboring countries.
In 2012, the Royal Thai Survey Department (RTSD) and Chiang Mai University released two local geoid models of Thailand, i.e., THAI12G and THAI12H, for public uses (Dumrongchai et al. 2012).THAI12G is the gravimetric geoid containing the long and medium wavelength information of the Earth Gravitational Model 2008 (EGM2008) (Pavlis et al. 2012).Most of the 3979 land gravity data stations used for the geoid computations were in the relative gravimetry Terr. Atmos. Ocean. Sci., Vol. 32, No. 5, Part II, 857-872, October 2021 campaigns from the 1960s to 1970s and tied to the International Gravity Standardization Net 1971 (IGSN71).The coordinates of the stations, including heights, were coarsely measured with unidentified positioning accuracies.As such, the gravity values may not represent the actual gravity field over areas of interest.Compared to EGM96 (Lemoine et al. 1998), EGM2008 shows an improvement of THAI12G accuracy in a few centimeters.Significant errors range from +8.7 to +119.6 cm, mainly in mountainous areas, bordering areas, and inaccessible areas with no terrestrial gravity measurements.THAI12H is the local hybrid geoid model, consistent with Kolak-1915vertical datum (Kolak-1915) through fitting 200 GPS/leveling co-points using a standard least-squares collocation (Moritz 1980).THAI12H yields an overall accuracy of 16 cm at 53 GPS/leveling checkpoints.
Although THAI12H provides potential accuracies as close as 5 cm in the Bangkok Metropolitan Region, significant errors up to 30 cm are found in other areas, especially, Chao Phraya basin, northern parts of Thailand (Dumrongchai and Duangdee 2019).These errors are mostly from the low point intensity of gravity measurements in those areas.The gravity values may be out of date and could not correctly represent the gravity field of the mass in the areas.Furthermore, the LaCoste and Romberg land relative gravimeters (Model G-1092) used by RTSD are outdated and may misread the gravity values measured.
In 2015, the Royal Thai Survey Department initiated Precise Geoid Model Project under Short-Term National Water Strategy Plan (fiscal years 2015 to 2017).The project aimed to produce a local geoid model covering the country's territory at centimeter-scale accuracy levels.RTSD planned to integrate the new geoid model into the national real-time kinematic GNSS network to support Thailand's height modernization system.Since the country had outdated gravity values, RTSD decided to establish a new fundamental absolute gravity network as a reference gravity frame for the first, second, and third-orders relative gravity networks.There were many people involved in the project and tremendous efforts in conducting terrestrial and airborne gravimetry across the country and numerical computation tasks.
Thailand Geoid Model 2017 (TGM2017) represents Thailand's first local geoid model using terrestrial and airborne gravity data.Our mission aims to produce TGM2017 using the most updated gravity data and the GNSS/leveling co-points in the national control networks.The methodology of geoid determination follows a classical Molodensky's approach (Heiskanen and Moritz 1967;Forsberg et al. 2000;Forsberg and Tscherning 2014;Featherstone et al. 2018;Vu et al. 2019).The contents of this paper consist of the overall results of terrestrial and airborne gravity gravimetry and relevant data in sections 2 and 3.All different types of gravity data were stochastically combined using least-squares collocation and downward continuation (Forsberg et al. 2007) in section 4. All areas outside Thailand were filled with GECO (Gilardoni et al. 2016) for lands and DTU13 (Andersen et al. 2015) for marine areas at an optimal gridding resolution.Next, the modified Stokes' algorithm was used for gravimetric geoid modeling, as described in section 5.The conversion surface was computed using the 299 GNSS/leveling co-points through least-squares collocation.In section 6, the accuracy of the hybrid geoid model, TGM2017, was assessed using the 100 GNSS/leveling checkpoints.Finally, the conclusions were summarized in section 7.

NEW GRAVITY DATA
Our goal is to achieve a local geoid model of Thailand at ten-centimeter accuracy or better anywhere if possible.According to Jekeli et al. (2009), the assessment of gravity data requirements, particularly data resolution in South Korea, shows that a three-arcminute (5.5 km) data resolution possibly causes 3.5-cm omission error in geoid undulations in rough areas.The error can be reduced by a few centimeters if topographic areas are relatively smooth, similar to the case of Thailand for smooth areas with lower 100m-elevation (Dumrongchai 2012).The computation of an accurate geoid model for Thailand requires dense gravity observations and precise and reliable gravimeters.All gravity data were collected across the country by terrestrial and airborne gravity surveys to generate an optimal gravity data set for geoid determination.Approximate 10000 terrestrial gravity points had been measured since the initiation of the geoid project by RTSD.Although the ground gravity surveys enable detailed characterization and high-resolution of the survey areas, their disadvantages are the limitations of inaccessible and remote areas, such as mountainous areas and marine areas.The airborne gravity surveys, on the other hand, provide better coverage and resolution in those areas.This section briefly discusses how to acquire the gravity data set, contributing to the geoid's short and medium wavelength component.

Terrestrial Gravity Data
The first gravity station, g 0 , located inside Gravity Measurement Building, Kalayanamitri Rd., Bangkok, was referred to Potsdam Gravity Datum in 1937 and converted to IGSN71 in 1971.More than 3979 relative gravity stations, referred to g 0 , were measured along main roads throughout the country using old-fashioned LaCoste and Romberg gravimeters.Each station was remotely scattered in a wide area of Thailand, having a total size of 513120 square kilometers.Some stations were destroyed, and might not be possible to be reobserved.Furthermore, the gravity data were too old, and the quality of them was obscure.RTSD decided to re-establish the new Thailand gravity networks for geoid determination by following the standards and specifications for geodetic control networks of 1984, Federal Geodetic Control Committee (FGCC1984) (Bossler 1984).The fundamental gravity network of Thailand has served as a reference gravity frame for all national gravity measurements.The gravity survey was carried out from April to July 2015 with 93 absolute gravity stations using Micro-g LaCoste A10 portable absolute gravimeter (Microg LaCoste Inc. 2008), as shown by the blue square dots in Fig. 1.The interval between stations was about 100 kilometers (km) and, mostly, co-located at the first-order leveling stations.Using the g9 software for data processing (Micro-g LaCoste Inc. 2012), we also performed a relevant series of corrections, such as transfer height correction, barometric height correction, earth tide and ocean load corrections, and polar motion correction.The absolute gravity values of the stations have the measurement accuracies ranging from 10.7 to 12.0 μGal.In 2020, we revisited these stations to ensure continuing stabilities and found that the long-term drift values changed less than 20 μGal (at least five years), following FGCC1984.
The 405 first-order relative gravity stations were extended from the 87 absolute gravity stations during the nationwide survey campaigns from June 2015 to May 2017, as shown in Fig. 1.We used Scintrex CG-5 portable relative gravimeters to acquire the gravity values of the stations.The survey lines mostly followed existing roads in the firstorder leveling network with loop misclosure errors of less than 30 μGal.The interval between the two stations ranged from 25 to 50 km, depending on topographic elevations.We assumed no significant effects of instrument drift, atmospheric mass attraction, and earth tides on gravity measurements to simplify adjustment computation.Thus, we treated the gravity values of the stations as the parameters to be estimated.The least-squares adjustment was run on the entire first-order gravity network, constrained by 87 absolute gravity control points.The results indicated that the propagation errors of the estimated gravity values in the network ranged from 3.4 to 13.3 μGal.
Due to the limited time of survey operations in the geoid project, we considered the third-order gravity survey network campaign extending from the first-order network, conducted from November 2015 to July 2017.The campaign aimed to establish gravity points in the gridded pattern at 10 to 25 square kilometers per station density.Each survey line started at one first-order gravity point and ended with another first-order point within one day or less.We defined the rejection criterion of 0.1 mGal misclosure for blunder detection and found that at least 20 of the total 1056 survey lines exceeding the criterion, requiring re-observations.During the field operations, the main problem was the longer 50-km survey lines, mostly in mountainous areas in the west part of the country, probably resulting from dynamic drifts of the gravimeters.We needed to reduce the drifts by implementing at least two second-order gravity points along the survey lines between two first-order points.Other problems were inaccessible areas such as rough terrains, swamps, rural areas, and restricted areas.A simple adjustment of each line was carried out to ensure that the observed gravities and the known (first-order) gravities of closing the controlling points agreed, similar to junction closures of a leveling line.The amount of correction was prorated along the line in proportion to the travel time from one endpoint.We finally had 9524 out of 9625 gravity measurement points in the network, as shown by pink square dots in Fig. 1.Most of the points were uniform grids in the Central Plain, where the Chao Phraya River watershed experienced many floods, especially the great flood of 2011.Other parts had the points less dense, mostly in existing road networks and public areas.All location points were positioned using the real-time kinematic GNSS survey with THAI12H for orthometric height determination.The horizontal and the vertical accuracies were about 5 and 10 cm, respectively, sufficient for geoid determination.

Airborne Gravity Data
Thailand has a shape like an ancient ax or a long trunk surrounding in a clockwise direction by the Andaman Sea, Myanmar, Laos, Cambodia, the Gulf of Thailand, and Malaysia.Over 30 % of Thailand's territory is mountainous and sea areas, where terrestrial gravity data are sparse or not available at all, as seen in Fig. 1.Consequently, determining a precise geoid model at centimeter accuracy levels may not be possible, particularly in those areas lacking gravity data.Airborne gravity technology plays a significant role in remotely collecting gravity data in inaccessible areas.RTSD conducted the airborne gravity survey campaign across the country from May 2016 to June 2017 (Dumrongchai et al. 2018).This survey was the first-ever airborne gravimetry in Thailand for geodetic purposes.However, no flights were allowed within 25 km from the border.The airborne surveys were carried out in seven block areas, based on the allocated time of the aircraft for supporting other missions and seasonal changes in different parts of the country, as seen in Fig. 2.
The Micro-g TAGS-6 gravimeter (serial number: S-120) was installed abroad Beechcraft Super King Air model B200 RTSD aircraft, equipped with a NovAtel DL-V3 GPS receiver as a timing and positioning unit (Micro-g LaCoste Inc. 2015).The TAGS-6 had a sampling data of 20 Hz, which allowed better GPS and gravity timing than the old-version, TAGS/Air III gravimeter, and faster aircraft speed to acquire gravity signals at 0.01 mGal resolution and 0.1 mGal accuracy.The TAGS-6 PiperPro software provided sophisticated data collection and processing during all flight operations.Three military airports were used and had gravity values at the airport parking spots tied to the first-order gravity stations by Scrintex-CG5 relative gravimeter.Each block area had at least two GPS/GNSS base stations used for the kinematic differential positioning of aircraft trajectories, determined with Waypoint GrafNav post-processing software (NovAtel Inc. 2014).The accuracy of the TAGS-6 GPS positioning results was about 10 cm.This value typically allowed vertical, filtered disturbing acceleration, and Eötvös corrections determined at better than 1 mGal (Forsberg et al. 2012).
The flight altitude of gravity survey was about 4000 m above mean sea level, theoretically resulting in the recoverable wavelengths of gravity signals at 3 km or shorter (Hwang et al. 2007).The aircraft, collecting data, cruised at a typical speed of about 200 knots, around 370 km per hour or 103 m per second.Since the sampling data was 20Hz, the TAGS-6 system recorded data at about 5m-interval along survey lines, enabling us to observe the effect of clear-air turbulence (CAT) on gravity data, which could deteriorate short-wavelength accuracy of a geoid model (however, the effect has still been under investigation and will not be addressed in this paper).The survey plane flew along tracks in a north-south direction with 10 km spacing and across tracks in an east-west direction at every 50 km for estimating crossover errors.The total number of along-and cross-tracks was 183 and 52 lines, respectively, including two repeated along tracks with bad data in Block I.The total flight distance was about 65000 km, and the total number of cross-over points was 999 points (Dumrongchai and Duangdee 2019).
The raw gravity data was denoised by a low pass fil-ter of 90-second using Aerograv post-processing software (Micro-g LaCoste Inc. 2011) for smooth flights with normal weather conditions.The 120-second filter was applied to the noisy data caused by unexpected conditions such as CAT, strong aircraft vibrations over very rough terrains, and monsoon weather.The processing scheme allowed gravity anomalies to be recovered in the range of 5.4 to 6.2 km spatial resolutions (half wavelengths) along tracks, with an airspeed of 200 knots.The root mean square (rms) misfit of cross-over differences and indicated noise level were 4.2 and 3.0 mGal, respectively.Only non-adjusted gravity data at cross-over points were used for geoid computation.More details in the cross-over adjustment of RTSD airborne gravity data can be found in Srimanee et al. (2020).In Fig. 2, the free-air anomalies varied from -88.7 to +55.2 mGal with a standard deviation of 16.2 mGal (see also Table 1).Most of the anomalies significantly changed in the northern region due to irregularly topographic features.It is also evident that the airborne anomalies correlate with the topography.

GNSS/Leveling Data
The RTSD horizontal geodetic networks are categorized into three levels as follows: (1) the zero-order network, (2) the primary network, and (3) the secondary network.The zero-order network consists of 7 GNSS stations as a local reference frame, realized in the International Terrestrial Reference Frame of 2008 (ITRF2008) at epoch 2013.81.The primary (first-order) network is extended from the zero-order network and comprised 19 GNSS stations with an interval of about 250 km for each station.
After the concurrence of the 9.2 Mw Sumatra-Andaman earthquake on the 26 th December of 2004, the network was readjusted and classified as class B by FGCC standard from 2005 (Satirapod et al. 2009).For the secondary (secondorder) network, there are 94 GNSS stations tied from the first-order network with a spacing interval of 50 to 100 km and an accuracy of about one part per million (ppm).RTSD has intensified the secondary network with 20 -50 km intervals with more than 600 stations since 2005.
The Kolak-1915 vertical datum has remained in the official vertical datum of Thailand.It was realized by averaging tidal observations between 1910 and 1915 from one tide-gauge station at latitude 11°47'42"N and longitude 99°48'58"E.The first-order vertical network contains 357 primary benchmarks of orthometric heights.There are 1428 secondary benchmarks extended from the primary control network.Because the vertical network was distributed in two spur shapes, RTSD independently performed the network adjustments in upper and lower areas by minimally constrained adjustments.The adjustments consequently caused datum inconsistency.The result of the network yielded a height accuracy of a few centimeters.
Only 412 GNSS stations are co-located at the firstorder leveling benchmarks.We exclude 13 GNSS stations due to large errors in the ellipsoidal heights.The possible causes may be, for instance, the height of instrument blunders, topographic or building obstructions, and 50 km or longer distances between reference stations and GNSS stations along borderlines.The remaining heights are accurate to better than 3 cm.Figure 3 shows the distribution of 299 GNSS/ leveling reference stations (dots) and 100 GNSS/leveling checkpoints (star dots) for geoid computations and evaluations, respectively.These stations are rather patchy, and their spacings ranged from 25 to 50 km.However, the stations are less dense in mountainous north and northwestern areas.

The Global Gravity Model GECO
GECO is a global gravity model computed by incorporating the GOCE-TIM-R5 global model into EGM2008 (Gilardoni et al. 2016).It optimally combines two sets of GOCE and EGM2008 spherical harmonic coefficients by setting up the error covariance matrices of two models for least-squares adjustment.GECO is more informative than EGM2008 in the areas where no ground gravity data are available at the period of EGM2008 computation.It improves the accuracy of EGM2008 in the low to medium wavelengths, corresponding to the reduction of commission and omission errors up to degree and order 359.From de-gree 360 to 2190, the GECO coefficients are the same as EGM2008.Comparing GECO with EGM2008 indicates the geoid improvement in the very high mountainous areas, for example, the Himalayas, the Andes, and the Congo basin.In the Southeast Asian region, the difference between GECO and EGM2008 geoid error up to degree and order 359 are about 5 cm or larger.For this reason, GECO was chosen as the reference geoid model in all computations in Thailand and neighboring countries, providing its structures of long (and feasibly some medium) wavelengths.The range of spherical harmonic coefficients for local geoid computations is degree and order 2 to 2159 with additional coefficients up to degree 2190 (available for public uses at http:// icgem.gfz-potsdam.de/home).

The Global Marine Gravity Model DTU13
The DTU13 global marine gravity model is a satellite altimetry-derived gravity model (Andersen et al. 2015).The data used for developing the model consist of altimetry data from old multi-mission satellite altimeters from September 1992 to December 2009, Cyosat-2 launched April 2010 (3 repeats of 369 days), and Jason-1 from April 2012 to May 2013.DTU13 was computed by means of the remove-restore procedure using EGM2008 up to degree and order 1960, and DTU07/EGM2008 (up to degree and order 100) mean dynamic topography.Performing cross-over adjustment and decreasing filtering results in improved medium-to-short wavelengths of the gravity field in marine areas.The spatial resolution of DTU13 is one arcminute grid.The DTU13 free-air gravity data were used for data padding in open sea areas, where no airborne and shipborne gravity data were available, particularly in the Gulf of Thailand and the Andaman sea.The DTU13 global marine gravity model is available for free download via https://ftp.space.dtu.dk/pub/DTU13.

The RTSD Digital Terrain Model Level 2 RTSD-DTED2
The RTSD digital terrain elevation data level 2 (DTED2) was used to generate a grid of the residual terrain model or RTM (Forsberg 1984) for terrain reduction.National Geospatial-Intelligence Agency or NGA (formerly National Imaginary and Mapping Agency or NIMA) contributed DTED2 to RTSD under the Thailand Major Mapping project in December 1999.The one-arc-second DTED2 covers the country's land areas from latitude 5.5N° to 21.0°N and longitude 96E° to 106°E with some areas along the border missing.The areas having no data were filled-in using SRTM 1 Arc-Second Global elevation data (downloaded from http://ww.usgs.gov in 2016).DTED2 has been frequently updated using orthometric heights and photogrammetric digital terrain model (DTM), particularly in void areas and changing terrains.The modified version of the original DTED2, called RTSD-DTED2, relatively more represents the topography in Thailand than SRTM 1 and usable for geoid determination, as seen in Fig. 1.However, it has still contained some erroneously registered height values.The RTSD-DTED2 was generalized to coarser resolution grids, e.g., 0.5, 1, and 2 arcminute grids, using the GRAVSOFT SELECT program (Forsberg and Tscherning 2014) to be prepared for RTM terrain corrections.

Data Augmentation
The high-resolution gravimetry in Thailand was required before transforming to regular grids as the input data for geoid computations through the Fast Fourier Transform (FFT) technique.Extending geoid computation area mitigated edge effects in FFT operations.Otherwise, either interpolation or extrapolation of existing gravity points onto a desired regular grid in vast data gaps significantly introduced unwanted errors to the computations.We used GECO and DTU13 for data padding in neighboring countries and marine areas to reduce geoid errors caused by the lack of gravity data in those areas, especially border areas and coastal areas [in future, we shall use a more recent marine (coastal) gravity field based on the modern altimetry technology, e.g., DTU17 (Andersen and Knudsen 2019), that is better suited for coastal geoid modeling].The augmentation of the measured gravity data sets with GECO and DTU13 free-air gravity data covered 3N° to 23°N and longitude 95E° to 108°E using the GRAVSOFT SELECT program and the GMT grdlandmark function (Wessel and Smith 2015).Modeling either the geoid or quasi-geoid was possibly problematic from the inconsistency of different data sources and data discontinuity.We padded the data away from Thailand's borderline and airborne gravity survey areas at least 10 km onto a regular grid suited for accurate geoid modeling, for instance, 0.5' × 0.5', 1' × 1', and 2' × 2' grids.Figure 4 shows the example of the augmentation of land and airborne free-air gravity anomalies with data padding on a regular ground grid.

DATA COMBINATION
Airborne gravimetry differs from terrestrial gravimetry in view of accuracy gain.The latter provides better accuracy due to the higher sensitivity of the instrument.However, terrestrial gravimetry is limited in inaccessible areas such as mountainous, coastal, and water areas, where airborne gravimetry has more advantages.The integration of the different gravity data sources requires an optimal method to mitigate errors due to, for instance, data inconsistency and instrument errors accumulated in the gravity data sets.This section describes the combination of the gravity data sets with optimal gridding of gravity data using least-squares collocation.The GECO and the DTU13 gravity data were padded in neighboring countries and marine areas, respectively.
Airborne gravity data sets were measured at the flight altitude of 4 km above mean sea level.The airborne data were needed to combine with terrestrial gravity data.Some parts of Thailand, e.g., north and west areas, were mountainous and inaccessible.The terrestrial gravity resolutions in those areas were not uniform, with data mostly following existing roads.Large data gap areas were in the order of over 50 km (27.8 arcminutes), possibly increasing omission errors in the areas where the geoid was computed.The data distributions were more uniform in other areas, and their resolutions varied from 5 to 10 km (2.5 to 5.5 arcminutes).As seen in Figs. 1 and 2, integrating the different types of gravity data requires an approach that optimally combines these data.Figure 5 shows the steps of heterogeneous data combination using least-squares collocation.We applied the residual terrain model (RTM) for terrain corrections using 30" DTM, averaging from 1" RTSD-DTED2.For the areas outside Thailand, including marine areas, we filled in with GECO and DTU13.The least-squares collocation (Moritz 1980;Forsberg and Tscherning 1981;Forsberg 1987) was used for the downward continuation of all available gravity data onto regular ground grids, e.g., 30" × 30", 1' × 1', and 2' × 2' grids.
All available terrestrial and airborne gravity data sets, Δg, can be combined and computed to acquire the predicted   gravity data values, g D , onto a ground grid using the leastsquares collocation as where C g g D D is the covariance matrices between predicted values and observed values, and C ΔgΔg is the covariance matrices among observed values.The symbol D ΔgΔg is the (diagonal) covariance matrix of observation noises under white noise assumption.The elements of the covariance matrices are obtained using a planar logarithmic covariance model (Forsberg 1987;Forsberg and Tscherning 2014), where the gravity data set are at two different altitudes, e.g., h 1 and h 2 , as which is expressed by three free parameters: C 0 is a variance of gravity anomalies, D is the Bjerhammar sphere depth, T is the long-wavelength attenuation factor.The d symbol is the distance between the predicted point and the observation point and i a is the weight factor at i.The GPFIT program from GRAVSOFT can estimate these parameters.
Combining all available gravity data sets at different altitudes follows a remove-predict-restore step separating the observed free-air gravity anomalies into three parts where Δg GGM is the free-air anomalies from the global geoid model, and Δg RTM is the RTM terrain effects.In this study, we generated Δg GGM from GECO up to degree and order 359, which corresponded to 50 km resolution (a half wavelength of geoid) at the topographic heights where the gravities were measured.The last term of Eq. (3), i.e., Δg res , is the residual free-air anomalies as an input of Eq. ( 1); we replace Δg, C ΔgΔg , and D ΔgΔg in Eq. ( 1) with Δg res and its corresponding covariance matrices.We removed the gravity terrain effects using the RTM terrain reduction method to reduce the downward continuation error.The short-wavelength part of the terrain effects was computed using prism integration with the mean crust density of 2.67 g cm -3 .The reduction made a stable solution and diminished topographic aliasing into the airborne data and, particularly, the terrestrial data lacking random distribution relative to the topography (Forsberg et al. 2007;Zhao et al. 2018).We used 30" DTM averaging from 1" RTSD-DTED2 for producing RTM terrain effects on ground and airborne surfaces.To obtain a suitable smooth height sur-face, we chose 15' × 15' moving window (27 km × 27 km) using SELECT from GRAVSOFT.The assigned window avoided repeatedly removing the lower frequency contents of terrain effects, which had already been taken into account by GECO.The RTM data was computed using terrain correction (TC) program from GRAVSOFT and applied to the gravity data sets.Before performing downward continuation, the RTM terrain effects at the flight altitude of 3 -4 km were filtered by the same filter used in the airborne data processing, i.e., low pass filter half-widths of 90 -120 seconds.
Table 1 shows the statistics of free-air gravity anomaly data sets with RTM terrain reductions.The airborne gravity anomalies vary from -88.7 to +55.2 mGal with a mean bias of -11.8 mGal.The standard deviations of the residual anomalies decrease to 7.8 and 11.6 mGal for airborne and terrestrial anomalies, respectively, whereas 16.2 and 19.5 mGal before reductions.These results indicate reductions of long and short wavelength gravity components by GECO and RTM.We filled GECO(720) up to degree and order 720 and DTU13 free-air anomalies in outermost land and marine areas, respectively (see Fig. 4), to reduce gravity discontinuity effects in border and coastal areas during the downward continuation processing.In addition, we included RTM terrain corrections in DTU13 due to minor coastal effects in satellite altimetry.The residual terrestrial anomalies, g res Ter D , generally are less than 30 mGal in magnitude except that large values occur in the mountainous areas, mainly in the north and the northwest parts of Thailand.These values imply that GECO may not sufficiently represent the gravity field in those areas where the elevations are higher than 1000 m (see Fig. 1).For the residual airborne, g res Air D , the mean value decreases to 0.1 mGal with the standard deviation (s.d.) of ±7.8 mGal.The decreasing s.d.values of g res Air

D
and g res Ter D indicate the success of terrain reductions on short-wavelength information in the data sets.
By definition of random process in the collocation (Moritz 1980;Forsberg 1987), the inputs in Eq. ( 1) are the residual anomalies with bias-free.Thus, removing the mean biases of all available gravity data was done before the downward continuation process.Assuming no correlation among observation errors (white noise process) yielded the D ΔgΔg diagonal covariance matrix with the assigned standard errors of terrestrial and airborne gravities of ±0.1 and ±3.0 mGal, respectively (see section 2).We computed the empirical covariance function from terrestrial and airborne residual anomalies.The empirical function then was fitted by the planar logarithmic covariance model of Eq. ( 2) using EMPCOV and GPFIT from GRAVSOFT, respectively, with three free parameters: C 0 = (9.0) 2 mGal 2 , D = 1 km, and T = 32 km. Figure 6 shows the relationship between the covariance functions and distances.The combined gravity data correlation is 7 to 10 km and diminishes after 20 -25 km.The empirical covariance plot shows a hump at around 40 -60 km, implying irregular distribution of gravity measurements when airborne gravities and terrestrial gravities were combined, for instance, in the north and northwest parts of Thailand.The elements of C ΔgΔg and C g g D D were computed among the observed and the predicted points.
In the final step, we downward continued the airborne and terrestrial residuals to ground level by leastsquares collocation in 1° × 1° blocks expanded with a 0.6° × 0.6° overlap between blocks using GPCOL1 from GRAVSOFT.There were three candidate sets of the reduced free-air anomaly grids, i.e., 30" × 30", 1' × 1', and 2' × 2' grid.Figure 7 shows the example of the reduced free-air anomalies after downward continuation to ground level.In Figs.7a and b, the reduced airborne anomalies are smoother than the terrestrial anomalies.This spatial pattern of the airborne anomalies reflects the longer wavelength contents of the gravity field than the terrestrial anomalies.The airborne anomalies agree with the terrestrial anomalies within 10 -20 mGal in magnitude.In fact, the residual terrestrial anomalies are only interpolated and extrapolated on a ground grid.They could be used as reference values for evaluating the quality of the airborne downward continuation.It is noted that there are good data agreements, mainly in the central and the northeast parts of the country, where high densities and regular distributions of land gravity measurements exist, as shown in Fig. 7c (see also Fig. 4).However, significant differences are found in the rough topographic areas of mountainous north and northwest parts, where the terrestrial gravity points are sparsed and measured only along existing roads.The combined all available gravity anomalies after downward continuation is illustrated in Fig. 7d.This combination clearly shows the advantage of airborne gravimetry along with coastal areas and a good transition to DTU13 satellite-altimetry anomalies in sea areas.

COMPUTATION OF THAI17G
The THAI17G gravimetric geoid modeling is theoretically based on Molodensky's geodetic boundary value problem that requires gravities on the topographic surface of Thailand.In practice, not the entire geoid model area is covered with gravity data.Therefore, all available data conducted by different types of gravimetric surveys are stochastically combined.This section describes the gravimetric quasi-geoid determination using the gravity data sets from airborne and terrestrial gravimetry with GECO(720) up to degree and 720 for neighboring country areas.The DTU13 gravity anomalies were filled in marine areas.The GECO(359) up to degree and order 359 was used to create residual free-air anomalies with RTM terrain reductions and downward continued to a ground grid so that Stokes' algorithm was applicable.Finally, converting the quasi-geoid to the geoid provided THAI17G.
In this project, the height anomaly, g, was computed through the modified Stokes' integral (Heiskanen and Moritz 1967).All computations were in the non-tidal system.With the usual remove-and-restore procedure, the residual height anomaly is defined as follows where R is the Earth's mean radius and c is the normal gravity at the telluroid.The Wong and Gore modification of the Stokes' kernel function, S MWG , (Wong and Gore 1969;Forsberg and Tscherning 2014) is given by with } is Legendre function with the spherical distance, }, between the point of computation and the surface element, dv.The linear taper coefficient, ( ) n a , is applied to avoid the influence of the residual quasi-geoid at long wavelengths.Low harmonics can be removed from ( ) S } up to degree N 1 and then linearly tapered to N 2 , which should be less than or equal to the maximum degree, N max , of a reference geoid model (e.g., EGM2008).
For applying the remove-predict-restore technique, the height anomaly is split into three components as GGM R TM res The residual gravity anomalies is converted to the residual quasi-geoid, res g , based on performing the F Fast Fourier Transformation (FFT) operation We applied the multi-band FFT method (Forsberg and Sideris 1993) implemented in GRAVSOFT for Eq. ( 8).Restoring the quasi-geoid effects of RTM terrain and GECO(359) yielded the height anomaly, g in Eq. ( 7).Finally, we had g converted to the N geoid height or undulation by the approximation relation (Heiskanen and Moritz 1967, pp. 327-328) where Δg B is the Bouguer gravity anomaly, c is the mean normal gravity along the normal plumb line between ellipsoid and telluroid, and H is the orthometric height.In this work, we approximated Bouguer gravity anomalies from a regular grid of free-air anomalies (as described in section 4) with the Bouguer reduction using a digital terrain model, e.g., RTSD-DTED2.The term c was computed using H from the digital terrain model, instead of normal height.All geoid computation steps were done by GRAV-SOFT using Geodetic Reference System 1980 (GRS80) as a reference ellipsoid (Moritz 2000).We computed the derived geoid heights from 399 GNSS/leveling co-points to determine the suitable values of N 1 and N 2 in Eq. ( 5) that best fitted the gravimetric geoid in Eq. ( 9).Based on multiple tests, the values of N 1 = 150 and N 2 = 360 indicated some low wavelength contents of GECO(359) remaining in the residual anomalies and provided the best geoid precision.Therefore, we used these values for the final gravimetric geoid computation.We also experimented whether a few long-wavelength contents of GECO(359) might be in the residuals by increasing N 2 .The geoid differences were relatively large after N 2 > 360.The residual quasi-geoid model was constructed using Eq. ( 8) with 100 % zero-padding to avoid the cyclic convolution effect in the model due to the FFT algorithm by SPFOUR.The RTM terrain effects were restored using TC.The GECO(359) normal heights generated at the topographic surface were added back to obtain g according to Eq. ( 7).The conversion of Eq. ( 9) yielded the geoid model.
As stated in section 4, we looked for an optimized grid resolution that provided the smallest errors of the comparisons between the geoid models and the derived geoids from 399 GNSS/leveling co-points, plotted in Fig. 3.We performed tests on several spatial resolutions and found three possible grid candidates.Table 2 lists the statistics of the undulation differences at 30" × 30", 1' × 1', and 2' × 2' grid points.Griding on one-arcminute resolution yielded the smallest standard deviation of 0.095 m.Therefore, we obtained the final gravimetric geoid, THAI17G.
We compared THAI17G, EGM2008, and GECO at several maximum degrees and orders with the derived geoid undulations from 399 GNSS/leveling co-points.Theoretically, these GNSS/leveling co-points are orthometric heights.Thus, EGM2008 and GECO included the height anomaly-to-geoid undulation correction of Eq. ( 9) to make the comparisons consistent.In Table 2, the statistics of the geoid differences indicate that the 1' resolution THAI17G provides the best geoid precision with the standard deviation of 9.5 cm, 20% accuracy improvement with respect to the GECO(720) geoid (12.0 cm in s.d.).The values of rootmean-square (rms.)indicate the fitness of two datums (36.2 -38.2 cm).THAI17G is in a global mean sea level similar to GECO and EGM2008, whereas Kolak-1915 vertical datum lies on a local mean sea level.The datum shift between two different mean sea levels is about +0.37 m.Such a datum inconsistency probably results from the global sea-level rise of about 21 -24 cm over a century (Lindsey 2021) and land subsidence, relative to Kolak-1915 based on averaging the water records from 1910 to 1915.Other factors may come from several error sources of the data used, e.g., leveling network adjustment, gravity data sets, digital terrain model, and long-wavelength errors of the global geoid models, propagating to the geoid.
Figure 8 shows the plots of the differences THAI17G minus GECO(720), and EGM2008(720) with values ranging from -0.064 m (in dark blue) to +0.969 m (in white red) (see Table 2).The large values appear in mountainous areas, in particular the north and northwest parts of the country.
The statistics of the THAI17G comparisons show an average offset of +37 cm with a standard deviation of 9.5 cm as listed in Table 2.It is evident that THAI17G contains more high-frequency contents than GECO and EGM2008.The high standard deviations in EGM2008 comparisons (about 14.0 -14.5 cm) reflect the lack of updated high-frequency data, which results in omission errors in the model.Figure 8b indicated larger values of the differences occurring along the areas located around longitude 100°E and latitude 14°N than the THAI17G and GECO(720) comparison in Fig. 8a.These disagreements reveal the improvement of some long and medium wavelength structures of THAI17G, including GECO(359), relative to EGM2008.There are high-frequency differences in marine areas, as seen in Fig. 8a (compared with Fig. 8b).This result is due to the inclusion of DTU13 altimetry-derived gravity anomalies in computing THAI17G.Although THAI17G gains the most accuracy improvement among other models evaluated, theoretically, it refers to a global reference system, i.e., an equipotential surface.The model differs from Kolak-1915 by a mean offset of 37 cm.Therefore, to be consistent with existing leveling, a stochastic process is needed to convert the THAI17G gravimetric geoid to the local sea level.

COMPUTATION OF TGM2017
In section 5, the validation of the THAI17G gravimetric geoid model with respect to 399 GNSS/leveling copoints indicates the existence of tilts and distortions of the model.There exists a mean offset of about 37 cm.Thus, to  be consistent with the Kolak-1915 vertical datum, a conversion surface is needed that related the model to the datum.We chose 299 GNSS/leveling co-points for constructing the surface using least-squares collocation.Finally, we had a hybrid geoid model that referred to Kolak-1915.The accuracy of the model was assessed by comparing its geoid undulations with the derived geoid undulations from 100 GNSS/ leveling checkpoints in Fig. 3.This hybrid model was officially named as Thailand Geoid Model 2017 (TGM2017) in August 2017.
Based on the previous results in section 5, the comparison of THAI17G and GNSS/leveling data set provides a means of estimating and removing possible systematic errors in the geoid model, leveling, and GNSS measurements.The model of combining THAI17G geoid undulation, N, GNSS derived ellipsoid height, h, and Kolak-1915 orthometric height, H, is defined by Since the least-squares collocation generally requires centered quantites, all systematice errors such as a bias and a tilted plane have to be removed from e.The zero-mean residuals are expressed as an input vector of zero-mean observed signals, i.e., x, with noise, n.The detrended signals, s J , (instead of g D ) is predicted on a 1' × 1' grid according to Eq. (1).We model a tilted plane using a simple form of polynomial surface (Dumrongchai et al. 2012).The relationship between e and tilts (considered as systematic errors) is given by where the parameters, a, b, and c, correspond to the datum shift and tilts of the geoid with respect to Kolak-1915 datum.These parameters are estimated using a simple least-squares adjustment as listed in Table 3.The mean offset between Kolak-1915 and THAI17G is +34.1 m.A tilt (-3.215 ppm) occurrs in the east-west direction, while a north-south tilt (+1.444 ppm) is relatively small.The covariance matrix elements were computed using a Gaussian (exponential) covariance function with the d distance between points, i.e., C 0 exp(-d/L), that optimally fit the empirical covariance function of x with the correlation length (L) = 35 km and the function variance at zero distance (C 0 ) = (12.6) 2 cm 2 .The trend surface was computed on a 1' × 1' grid using the parameters listed in Table 3.This surface was restored to the 1' × 1' grid of s J , which was predicted using least-squares collocation, similar to Eq. ( 1), to provide the conversion surface.Removing the (conversion) surface from THAI17G produced the final hybrid model, TGM2017. Figure 9a shows the conversion surface, which is relatively smooth and similar to the residual geoid undulations in Fig. 9b, but it is not reliable, along with the border stripe areas of Thailand.
Finally, TGM2017 was assessed using 100 GNSS/leveling checkpoints in Fig. 3. Table 4 lists the statistics of the differences between TGM2017 and Kolak-1915 geoid undulations.It shows the 4.9-cm rms of fit, which agrees with the 5.1-cm rms from the 299 co-points used for constructing the conversion surface.This comparison confirms that the conversion process is successful.The histogram of the geoid differences at 100 checkpoints is plotted in Fig. 10.It is apparent that the differences vary in the bound of ±10 cm and mostly cluster within -5 to +7 cm.However, a few differences are out of the bound values, around the border and the mountainous areas.The largest value of 14.9 cm is found in mountainous northwest Thailand (the purple dot at around latitude 17°N and longitude 98°E).This is likely due to a lack of gravity measurements, accumulation errors in leveling networks, and ellipsoid height errors.These make TGM2017 more uncertain in the areas near the border.The conversion surface is generally constructed to fit THAI17G to GNSS/leveling data and transfer Kolak-1915 across the country.The surface seemingly accommodates systematic errors, for instance, datum distortions and long-wavelength ellipsoid height errors, in TGM2017.As described in section 3, the 299 GNSS ellipsoid heights used for the surface computation are in the national horizontal network.Thus, if we measure the GNSS ellipsoid heights by differential positioning ways, these systematic height errors will be significantly reduced by referring to the same erroneous network.Consequently, TGM2017 leads to reliable orthometric heights in Kolak-1915 vertical datum at a 10-cm accuracy level or better.Checkpoints 100 -10.3 14.9 1.1 4.5 4.9 Computing points 299 -17.6 17.0 0.0 5.0 5.1 Table 4.The statistics of geoid differences between TGM2017 and GNSS/leveling co-points (units in cm).405 first-order and 9524 third-order relative gravity stations with the estimated accuracies of about 20μGal, 30μGal, and 0.1 mGal, respectively.Thailand was covered by airborne gravity data from the airborne survey until June 2017 from May 2016 at 10 km along-track spacing and a nominal flight altitude of 4000 m above mean sea level.The overall accuracy of the final airborne data set is 3.0 mGal with no crossover adjustment performed.TGM2017 was created in August 2017 based on the THAI17G gravimetric geoid and 299 GNSS heights on leveled benchmarks, referred to Kolak-1915 vertical datum.The 1" RTSD-DTED2 was used to compute RTM terrain effects on land and airborne gravity data.All available gravity data were downward continued onto a one-arcminute ground grid using least-squares collocation with GECO and DTU13, filled in neighboring countries and marine areas.A quasi-geoid model was constructed by multi-band spherical FFT and then converted to THAI17G referenced in a global mean sea level.THAI17G was fitted to 299 GNSS/leveling co-points.All computation works were done by GRAV-SOFT.TGM2017 was achievable and assessed by 100 GNSS/leveling checkpoints.The differences GNSS checkpoints minus TGM2017 vary in the bound of ±10 cm and mostly cluster within -5 to +7 cm with 4.9-cm rms.Large errors mostly appear in mountainous areas, particularly the northwest areas, where gravity data are not available or insufficient.With an overall accuracy of 10 cm or better, TGM2017 has been planned to be applicable as a reference for a height modernization system of geodetic infrastructure and, subsequently, supports the GNSS continuously operating reference station network of Thailand.

Fig. 1 .
Fig. 1.Thailand gravity network stations.Fig. 2. Seven airborne gravimetry block areas and free-air anomalies at flight lines in Thailand (units in mGal).
model for Thailand has been successfully computed at RTSD and Chiang Mai University.The model is a hybrid geoid model that has officially been called Thailand Geoid Model 2017 (TGM2017), covering the area 3° to 23°N latitude and 95° to 108°E longitude.The model is based on two new gravity data sets from terrestrial and airborne gravimetry campaigns from April 2015 through June 2017.All land gravity data refers to Thailand gravity networks consisting of 87 absolute gravity stations and Parameter Estimated value Datum shift or bias +34.1 cm tilt in East-West direction -3.315 ppm tilt in North-South direction +1.444 ppm Fig. 9. (a) The conversion surface relating THAI17G gravimetric geoid model to the TGM2017 hybrid geoid model; (b) the differences 299 GNSS/ leveling geoid minus THAI17G (contour interval of 0.02 m).

Table 1 .
Statistics of gravity data sets (units in mGal).

Table 2 .
The evaluations of THAI17G and global geoid models at 399 GNSS/ leveling co-points (units in m).

Table 3 .
The estimation of trend parameters.