Determination of the geopotential value on the permanent GNSS stations in Vietnam based on the Geodetic Boundary Value Problem approach

SUMMARY

In the realisation of the International Height Reference System, the determination of the geopotential value and its variations in time plays an important role. In this study, the geodetic boundary value problem approach is applied for direct determination of the gravity potential value using a GOCE global gravity field model enhanced with terrestrial gravity data. This determination is carried out on the Global Navigation Satellite System-Continuously Operating Reference Stations (GNSS-CORS) stations towards the realisation of the International Height Reference System in Vietnam. First, the effects of the GOCE global gravity field model omission error, the indirect bias term on the disturbing potential and the systematic cumulative errors in levelling data are estimated. These errors affect the estimated geopotential value. The results calculated on the GNSS/levelling points show that the effect of the GOCE DIR-R5 (up to degree/order 260) omission error on the offset potential value is quite significant. This effect was eliminated using high-resolution terrestrial gravity data using the remove-compute-restore technique. The indirect bias term on the disturbing potential can be safely neglected by using a GOCE global gravity field model for degrees higher than 60 for this study region. The systematic cumulative errors in levelling data can be modelled and removed using a third-order polynomial model. Then, the mean zero-height gravity potential of the Vietnam local vertical datum is estimated equal to |${\rm{W}}_0^{{\rm{LVD}}}$| = 62 636 846.69 m² s^–2 with standard deviation of 0.70 m² s^–2 based on the proposed approach. Finally, the geodetic boundary value problem approach was used to determine the geopotential on the surface of three GNSS-CORS stations in Vietnam. Based on time-series of the vertical component derived from the GNSS observations as well as InSAR data, temporal variations in geopotential are also estimated on these permanent GNSS stations. The purpose is to monitor deformation of the vertical datum. The results indicate that the geopotential value needs to be monitored and determined with the time-dependent component on the three Vietnamese permanent GNSS stations for a vertical datum. These stations may contribute to increase the density of reference points in the International Terrestrial Reference Frame, which is being researched and implemented by the International Association of Geodesy.

Geopotential theory, Reference systems, Satellite geodesy, Satellite gravity, Time variable gravity

1 INTRODUCTION

Until relatively recently, reference datums in the world were invariable in time, and ‘control points’ in such reference datums were generally considered stable and reliable over decades. With the advent of space geodesy, we have the observing capability to understand the physical changes occurring on the Earth and to realise a dynamic reference datum. A dynamic International Terrestrial Reference System (ITRS) and its realisation, the International Terrestrial Reference Frame (ITRF), based on a global network of Satellite Laser Ranging (SLR), Lunar Laser Ranging (LLR), Very Long Baseline Interferometry (VLBI) and Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS) stations colocated with Global Navigation Satellite System (GNSS, Altamimi et al. 2002; Altamimi & Collilieux 2009), indicated change in such control points at the level of centimetres per year. At the request of the International Association of Geodesy (IAG) resolutions (Drewes et al. 2016), the International Height Reference Frame (IHRF), a realisation of the International Height Reference System (IHRS) is also defined and realised following the same structure as that of the ITRF (Ihde et al. 2017). Complementary to the IHRS/IHRF initiative, an International Gravity Reference System (IGRS) and its realisation, the International Gravity Reference Frame (IGRF), have also been recently defined by IAG for the gravity and geopotential datum (Wilmes et al. 2017; Wziontek et al. 2021). The IHRS/IHRF includes a geometric component given by a coordinate vector (X) in ITRS/ITRF and a physical component given by the determination of the potential value (W) or geopotential number (C) with respect to the conventional value (W₀) at (X). According to Resolution No. 1 of the IAG (Drewes et al. 2016), the international conventional reference gravity potential is equal to 62 636.853.4 m² s² (Sánchez et al. 2016) for the purpose of implementing the IHRS. Moreover, besides estimation of the physical component, its variation in time (⁠|$\dot{W}$|⁠) should also be determined like the variation in the geometric component (⁠|$\dot{X}$|⁠) in the ITRF, to continuously monitor deformation as well as regularly update the height reference frame. Temporal geometric height changes can be accurately determined from observations on permanent GNSS stations whilst deformation of the topography around these GNSS stations can be estimated by Interferometric Synthetic Aperture Radar (InSAR). Thus, the variations of the geopotential value in time can be determined (even if imperfectly because of mass change in the crust, e.g. aquifer depletion) on these GNSS Continuously Operating Reference Stations (CORS). Therefore, the realisation of the IHRS should be conducted on the GNSS-CORS stations as reference points, and the geopotential value and its variations at such stations should be determined on the Earth's surface and not on the geoid or quasigeoid (Sánchez et al. 2019). In 2017, the first proposal for the IHRF reference network included 163 stations colocated with GNSS worldwide, but with poor coverage in particular over Africa and Asia (Sánchez 2017 and Sánchez 2019). Presently, the implementation of the IHRF is still under discussion (Ihde et al. 2017; Sánchez 2019). At the national or regional scale, despite vertical movement of the Earth, most countries do not model the vertical datum change. The U.S. National Geodetic Survey (NGS), a pioneer in the matter, is in the process of building a dynamic vertical datum for North America by providing a dynamic component of the geoid through the Geoid Monitoring Service (report of NGS 2019). Together with the realisation of the IHRS on a global scale, the implementation of the IHRS on the GNSS-CORS stations on a national or regional scale is an important task in which the determination of potential plays a key role. It may contribute to increase the density of reference points in the IHRF, especially over poorly covered areas.

There are two categories that should be considered in potential changes with time: (1) changes in potential related to Mean Sea Level (MSL) rise and (2) changes in mass redistribution leading to potential variations. The geoid is the equipotential surface which best fits to Global Mean Sea Level (GMSL, Gauss 1828; Listing 1873). However, presently GMSL is rising at a rate of approximately 3.2 mm yr^–1 (IPCC 2014), so the value of the gravity potential (W₀) should also change to maintain the best fit to GMSL. The variation of gravity potential caused by sea level rise was studied by Burša et al. (1997), Ardalan et al. (2002), Burša et al. (2007) and Dayoub et al. (2012). Sánchez et al. (2016) recommended using a potential value obtained for a certain epoch as the reference value W₀ and monitoring the changes of the MSL, that is changes in potential value at the sea surface W_S. When large differences appear between W₀ and W_S (e.g. > 2 m² s^–2), the adopted W₀ should be updated. At national or regional level, Local Vertical Datum (LVD) refers to local MSL. Analogously, the geopotential of a LVD (⁠|${\rm{W}}_0^{{\rm{LVD}}}$|⁠) should also change in order to accommodate the change in local MSL. The second category is the focus of this study. Change in elevation, for example due to land subsidence or uplift, will lead to changing redistribution of masses within the Earth, and its gravity field will change as well. One has to remeasure gravity to determine these changes, but like levelling this is a costly and time consuming effort. Therefore, an alternative method is applied in this study, which is based on forward modelling of residual terrain effects to determine temporal gravity variations as well as disturbing potential variations due to observed changes in topography.

There are more than 100 LVDs in the world nowadays. The main mission of the IHRS realisation is to integrate these existing local height systems into the IHRS by determining for each LVD the geopotential value (⁠|${\rm{W}}_0^{{\rm{LVD}}}$|⁠) or potential difference (δW^LVD = W₀–|${\rm{W}}_0^{{\rm{LVD}}}$|⁠, the IHRS vertical coordinates, Ihde et al. 2017) with respect to the equipotential surface of the Earth's gravity field realised by a conventional value (W₀) (Sánchez et al. 2016). This procedure is known as Height System Unification (HSU, Rapp & Balasubramania 1992; Ihde et al. 2000; Grigoriadis et al. 2014; Vergos et al. 2018) and its products are the vertical datum parameters, which are the potential and/or the potential differences between local and global reference level. In Vu et al. (2020), the zero-height geopotential value was estimated for the LVD of Vietnam using a local gravimetric-only quasigeoid and height anomalies derived from GNSS/levelling data. However, tilts, caused by accumulated systematic levelling errors and land subsidence, are known to have affected the GNSS/levelling data in Vietnam. The latter effects have been partly corrected through modelling based on the land subsidence derived from InSAR, but residual effects, due to the inconsistency between levelling and GNSS measurements [e.g. differences in accuracy level and tide systems (Sánchez 2012)], remain in the GNSS/levelling data. Based on the Geodetic Boundary Value Problem (GBVP) approach (Hofmann-Wellenhof & Moritz 2006), the geopotential value on the surface can be directly determined using the gravity data without height anomalies derived from GNSS/levelling data. As a result, the inconsistency that existed between the levelling and GNSS measurements is completely removed.

In this study, we first perform direct determination of the zero-height geopotential value of the LVD on the GNSS/levelling points using the GBVP approach to validate the proposed method as well as to show its advantages by comparing with the results derived from the GNSS/levelling data and gravimetric-only quasigeoid in Vu et al. (2020). Then the GBVP approach will be used for determination of the geopotential values including variations in time on the GNSS-CORS stations. This work aims to prepare the integration of the existing national height and gravity systems of Vietnam into the IHRS/IGRS. It also contributes to enhancing the density and distribution of reference stations in this part of the world, where the IHRF/IGRF has few reference stations.

2 DETERMINATION OF THE GEOPOTENTIAL VALUE

2.1 A review of approaches used for determination of the geopotential value

The realisation of the IHRS, that is the determination of the geopotential value of a reference point on the Earth's surface, is studied extensively in the geodetic literature. There are basically three different approaches for this realisation (Sánchez 2019). In this section, these approaches are briefly reviewed:

1. Using levelling and gravity reductions, the geopotential of the point P on the surface can be determined as follows (Fig. 1):

   $$\begin{eqnarray*} {W_P} = W_0^{LVD} - {C_P}, \end{eqnarray*}$$

   (1)

   where |${W_P}$| and |$W_0^{LVD}$| are the geopotential values on the Earth's surface and the LVD, respectively, |${C_P}$| is the geopotential number which can be estimated by combining spirit levelling with gravity measurements along the levelling line between the two benchmarks (Hofmann-Wellenhof & Moritz 2006).
   To determine W_P using eq. (1), it is necessary to know |${{W}}_0^{{\rm{LVD}}}$|⁠, which refers to the equipotential surface of the Earth's gravity field realised by a conventional value (W₀) and can be determined using GNSS/levelling data as follows (Tocho & Vergos 2016; Vergos et al. 2018):

   $$\begin{eqnarray*} W_0^{LVD} = {W_0} - \frac{{\mathop \sum \nolimits_{i = 1}^m \Delta {C_i}}}{m} = {W_0} - \delta {W^{LVD}} \end{eqnarray*}$$

   (2)

   where ∆C_i is the geopotential number of the GNSS/levelling points, relative to the W₀, and m is the number of the GNSS/levelling points. ∆C_i is given by the differences between height anomalies from GNSS/levelling measurements and those derived from quasigeoid, and the mean normal gravity value.

   The accuracy of this approach depends on the quality of the levelling data (Sánchez 2012; Sánchez 2019). Additionally, it cannot determine the variation in geopotential because continuous levelling measurements to determine the variation in levelling height are unavailable. It is well known that levelling contains distortions caused by vertical crustal movement, and systematic cumulative errors when levelling over long distances (Entin 1959; Vanicek et al. 1980). Moreover, to determine the geopotential of the point P on the Earth's surface (W_P) using eq. (1), the potential difference (C_P) determined using levelling data is once again affected by the errors of these levelling data. Therefore, this approach is unsuitable for the precise determination of geopotential on the surface, that is the realisation of the IHRS. It is suitable to determine the difference (⁠|$\delta {W^{LVD}}$|⁠) as well as the zero-height geopotential value of LVD (⁠|$W_0^{LVD}$|⁠) to integrate existing local height systems into the IHRS.

2. The availability of high-resolution Global Gravity field Models (GGMs), such as EGM2008 (Pavlis et al. 2012) or EIGEN-6C4 (Förste et al. 2014), which makes it possible to directly compute the geopotential of the ITRF coordinates X of any point via the spherical harmonic expansion equation (Hofmann-Wellenhof & Moritz 2006). However, biases in the LVDs lead to biases in the gravity data that were used in computing the GGMs For this reason, GGMs cannot be used to directly compute the geopotential value. More specifically, according to Rummel et al. (2014) these models have the expected mean accuracy of about ±0.4 to ±0.6 m² s² (equivalent to ±4 to ±6 cm) in well-surveyed regions, and about ±2 to ±4 m² s² (±20 to ±40 cm) with extreme cases of ±10 m² s² (±1 m) in sparsely surveyed regions. Vietnam is a representative of the sparsely surveyed regions, where fill-in data were used, and the estimated accuracy of EGM2008 and EIGEN-6C4 is 19 and 29 cm in standard deviation (STD) when compared to the GNSS/levelling data, respectively (Vu et al. 2019). These accuracies are insufficient to estimate the geopotential value in Vietnam.

3. Directly determine the potential value using the normal potential (U_P) and the disturbing potential (T_P) as follows (Hofmann-Wellenhof & Moritz 2006):

   $$\begin{eqnarray*} {W_P} = {U_P} + {T_P}.\end{eqnarray*}$$

   (3)

The normal potential is computed with the reference level ellipsoid (Hofmann-Wellenhof & Moritz 2006), and the disturbing potential is determined based on the GBVP approach using the gravity data. In this approach, the core is to determine the disturbing potential value by the GBVP approach applying the Remove-Compute-Restore (RCR) technique (Sansò & Sideris 2013; Barzaghi 2016), which will be described in Section 2.2. Thanks to the Gravity field and steady-state Ocean Circulation Explorer (GOCE) mission (Drinkwater et al. 2003), global geoids with an accuracy of 1–2 cm and gravity field models with an accuracy of 1 mGal at a spatial resolution of approximately 100 km are available to compute potential values. In this study, terrestrial gravity data are used to reduce the omission error of the GOCE model by increasing its resolution. Using the modern theory of Molodensky (Molodensky et al. 1962) in the GBVP approach, the determination of W_P is straightforward on the Earth's surface, and is unaffected by levelling errors or inconsistency between levelling and GNSS measurements. Additionally, from the potential values (W_P) on the reference points we can determine the zero-height geopotential value of LVD (⁠|${{W}}_0^{{\rm{LVD}}}$|⁠) based on the physical height of these reference points.

Figure 1.

Relation between global and local gravity potential.

From the geopotential values on the surface (⁠|${W_i}$|⁠) of each point i, the zero-height geopotential value for the LVD on this point |${( {W_0^{LVD}} )_i}$| can be determined based on physical height (⁠|$H_i^*$|⁠) as follows:

$$\begin{eqnarray*} {\left( {W_0^{LVD}} \right)_i} = {W_i} + H_i^* {\bar{\gamma }_i}, \end{eqnarray*}$$

(4)

where |$\bar{\gamma }$| is the mean normal gravity value along the plumb line from the reference ellipsoid to the telluroid (Hofmann-Wellenhof & Moritz 2006, eqs 4–60), and H* is the normal height. For each GNSS/levelling point we calculate a geopotential value, and the mean zero-height geopotential value for the LVD is obtained by:

$$\begin{eqnarray*} W_0^{LVD} = \frac{{\mathop \sum \nolimits_{i = 1}^m {W_i} + H_i^* {{\bar{\gamma }}_i}}}{m}, \end{eqnarray*}$$

(5)

where m is the number of the GNSS/levelling points. The zero-height gravity potential value is calculated from the physical heights, but not the height anomalies. This is advantageous since only the normal heights (⁠|${H^*}$|⁠) are taken into account, thereby avoiding any errors due to inconsistency in the difference between GNSS ellipsoidal height (⁠|$h),$| physical height derived from levelling and quasigeoid height (ζ). Moreover, this approach is impervious to the effects of land subsidence which occurred during the temporal separation between GNSS and levelling measurement.

2.2 Determination of the disturbing potential based on the GBVP approach

As described in Section 2.1, in the GBVP approach, the disturbing potential needs to be estimated. The spherical approximation of the fundamental boundary condition is as follows (Hofmann-Wellenhof & Moritz 2006):

$$\begin{eqnarray*} \frac{{\partial {\rm{T}}}}{{\partial {\rm{r}}}} + \frac{2}{{\rm{R}}}{\rm{T}} + \Delta {\rm{g}} = 0, \end{eqnarray*}$$

(6)

where r is the radius vector; R is the reference radius. Stokes’ formula is used to determine the disturbing potential from gravity data:

$$\begin{eqnarray*} T = \frac{{G\delta M}}{R} + \frac{R}{{4\pi }}\iint_\sigma {\Delta g\ S\left( \psi \right){\rm d}\sigma } \end{eqnarray*}$$

(7)

where S(ψ) is the Stokes function. |$\delta M$| is the difference between the mass of the Earth and that of the ellipsoid.

The Free-air terrestrial gravity anomalies are derived from measured gravity (g) and normal gravity (γ). The Free-air reduction in spherical approximation is calculated from the normal height (for the purpose of gravity reduction from the surface to the quasigeoid) according to Hofmann-Wellenhof & Moritz (2006) as follows:

$$\begin{eqnarray*} \Delta g = {\rm{ }}g - \frac{{\partial \gamma }}{{\partial {H^{\rm{*}}}}}{H^{\rm{*}}} - {\rm{ }}\gamma = g + \frac{{2\bar{\gamma }}}{r}\frac{C}{{\bar{\gamma }}} - \gamma = g + \frac{{2C}}{r} - \gamma . \end{eqnarray*}$$

(8)

The normal height (H*) is biased due to datum offset (⁠|$\delta {H^{\rm{*}}}$|⁠) between LVD and the global equipotential surface, so the gravity anomalies will also be biased by:

$$\begin{eqnarray*} \Delta {g^b} &=& \frac{{\partial \gamma }}{{\partial {H^{\rm{*}}}}} \delta {H^{\rm{*}}} = - \frac{{2\bar{\gamma }}}{r} \frac{{\delta {W^{LVD}}}}{{\bar{\gamma }}} = - \frac{{2\delta {W^{LVD}}}}{r} \nonumber \\ &=& 0.3086\, \delta {H^{\rm{*}}}. \end{eqnarray*}$$

(9)

So, it should also be pointed out that

$$\begin{eqnarray*} \Delta {{\rm{g}}^{{\rm{ub}}}} = \Delta {\rm{g}} - {\rm{ }}\Delta {{\rm{g}}^{\rm{b}}} = {\rm{ }}\Delta {\rm{g}} + {\rm{ }}\frac{{2\delta {W^{LVD}}}}{r} = \Delta {\rm{g}} - 0.3086\, \delta {H^{\rm{*}}} \end{eqnarray*}$$

(10)

is the unbiased (ub) gravity anomaly after correction for the datum offset.

Insertion of (10) into (7) gives:

$$\begin{eqnarray*} T = \frac{{\delta \left( {GM} \right)}}{R} + \frac{R}{{4\pi }}\iint_\sigma {\left( {\Delta g - 0.3086\delta {H^{\rm{*}}}} \right) S\left( \psi \right){\rm d}\sigma }. \end{eqnarray*}$$

(11)

Setting |${T_0} = \ \frac{{\delta ( {GM} )}}{R}$| is called the zero-degree term of the disturbing potential (Hofmann-Wellenhof & Moritz 2006), |$\ {T^{Stokes}} = \ \frac{R}{{4\pi }}\iint_\sigma {\Delta g\ S( \psi ){\rm d}\sigma }$|⁠, and |${T^{ind}} = \ \frac{R}{{4\pi }}\iint_\sigma {\frac{{2\delta {W^{LVD}}}}{{\rm{R}}}\ S( \psi ){\rm d}\sigma }$| is called the indirect bias term on the disturbing potential (hereinafter referred to as indirect potential bias term), which is affected by the offset from local (⁠|$W_0^{\rm LVD}$|⁠) to global datum (W₀). Eq. (11) can now be expressed as follows:

$$\begin{eqnarray*} {\rm{T}} = {{\rm{T}}_0} + {{\rm{T}}^{{\rm{Stokes}}}} + {{\rm{T}}^{{\rm{ind}}}}. \end{eqnarray*}$$

(12)

The datum offset |$\delta {H^{\rm{*}}}$| is unknown prior to height unification and only the biased anomalies ∆g can be used to determine disturbing potential from terrestrial gravity data. So we cannot calculate T^ind. Amos & Featherstone (2009) have solved this issue by determining and correcting datum offsets for different local datums in New Zealand in an iterative manner. In this study, the effect of the indirect potential bias term will be assessed for Vietnam and its surrounding areas through simulation, presented in Section 4.

The RCR technique is used to calculate T^Stokes. It is realised by summation of three terms:

$$\begin{eqnarray*} {{\rm{T}}^{{\rm{Stokes}}}} = {{\rm{T}}_{{\rm{GGM}}}}{\rm{ }} + {{\rm{T}}_{{\rm{RTM}}}} + {{\rm{T}}_{{\rm{res}}}}, \end{eqnarray*}$$

(13)

where T_GGM is computed using a GGM, T_RTM expresses the Residual Terrain Model (RTM) effect (Forsberg 1984) and T_res is computed from residual gravity anomalies using the Stokes’ integral eq. (7). The residual gravity anomalies for determining T_res are computed as follows:

$$\begin{eqnarray*} \Delta {{\rm{g}}_{{\rm{res}}}} = \Delta {\rm{g}} - \Delta {{\rm{g}}_{{\rm{GGM}}}} - \Delta {{\rm{g}}_{{\rm{RTM}}}}, \end{eqnarray*}$$

(14)

where ∆g is the Free-air gravity anomaly, |$\Delta {{\rm{g}}_{{\rm{GGM}}}}$| is the gravity anomaly computed with a GGM and |$\Delta {{\rm{g}}_{{\rm{RTM}}}}$| is the RTM effect on the gravity anomaly. Stokes’ formula requires that gravity data are used over the whole Earth. However, in the RCR technique, Stokes’ integral is evaluated using regional gravity data over a spherical cap of limited radius around the computation point and a selected GGM is used for the outer zones. This causes truncation errors. Therefore, when computing the residual disturbing potential obtained from Stokes’ integral the Stokes kernel should be modified to reduce these truncation errors (Molodensky et al. 1962; Witte 1967; Wong & Gore 1969; Meissl 1971; Heck 1987; Vaníček & Sjöberg 1991; Featherstone et al. 1998). According to Wong & Gore (1969) the modified Stokes kernel should not contain degrees lower than N (N is usually chosen smaller than N_max of the GGM used in the removal procedure) by removing of the low-degree Legendre polynomials:

$$\begin{eqnarray*} {S^{mod}} \left( \psi \right) = S\left( \psi \right) - {S^N}\left( \psi \right). \end{eqnarray*}$$

(15)

A smoother Wong-Gore (WG) modification is commonly used by completely removing low harmonics up to degree N₁, that is the influence of the long wavelengths in the local data is eliminated, and then linearly tapered to N₂. We have developed a code for computation of the disturbing potential by means of Stokes’ integration using the WG modification of the Stokes’ kernel function using Free-air gravity anomaly data on the surface according to Molodensky's theory.

3 DATA FOR THE DETERMINATION OF POTENTIAL VALUES

3.1 Grid of residual gravity anomalies

In Vu et al. (2019), a 5 × 5′ grid of residual gravity anomalies was computed for Vietnam and its surrounding areas to estimate the gravimetric quasigeoid over this region based on new terrestrial gravity data. A set of 29 121 land gravity measurements provided by the Institute of Geophysics (IGP)-Vietnam Academy of Science and Technology (VAST), the Vietnam Institute of Geodesy and Cartography (VIGAC) and the Bureau Gravimétrique International (BGI; Bonvalot 2020), was used in combination with fill-in data where no gravity data existed, for example mountainous or coastal areas. A mixed model [GOCE DIR5 (Bruinsma et al. 2014) plus EGM2008 (Pavlis et al. 2012)], called the mixed DIR/EGM model, up to degree/order 2159, plus Residual Terrain Model (RTM) effects, and gravity field derived from altimetry satellites, were used to provide the fill-in information over land and marine areas, respectively. The mixed model up to degree/order 719 was also used for the removal of the long and medium wavelengths in the RCR procedure, whereas the short ones were removed using the RTM effects. The 90-m resolution SRTM3arc_v4.1 (Farr et al. 2007) was used as the detailed Digital Terrain Model (DTM) over land areas in computing the RTM effects, while the 15″ resolution Digital Bathymetry Model (DBM) SRTM15arc_plus (Becker et al. 2009) was used over sea, together called the mixed SRTM model (see Vu et al. 2021, for details). After the remove procedure, the residual height anomalies were interpolated to grid nodes using GRAVSOFT GEOGRID (Forsberg & Tscherning 2008) based on the Least-Squares Collocation (LSC) method. In order to avoid step-like discontinuities when merging the fill-in and measurement grids, a transition procedure was applied (see Vu et al. 2019; Vu 2020, for details). Fig. 2(a) shows the merged gravity grid which was used in this study for the computation of the residual disturbing potential on the Earth's surface. Similar to Vu et al. (2019), both mixed models, DIR/EGM and SRTM, are used in the RCR procedure in this study.

Figure 2.

(a) Data sources of the gravity grid and (b) GNSS/levelling data: yellow inverted triangles are first- and second-order levelling benchmarks, purple dots denote third-order ones whereas red triangles are the GNSS-CORS stations.

3.2 GNSS/levelling data

From 2001 to 2003, the Vietnam national levelling network was remeasured, and then it was readjusted in 2007 using the local MSL over 1950 to 2005 at the Hon Dau tide gauge station (Pham 2009). From 2009 to 2010, the Vietnam Department of Surveying and Mapping (VDSM) carried out GNSS observations on the levelling points. The GNSS baselines were observed using dual-frequency instruments in static mode with a minimum measurement time of 6 hr per session. The GNSS data were processed with the Bernese software to obtain WGS84 ellipsoidal heights. A total of 812 GNSS/levelling observations was produced in that study. The GNSS/levelling data include horizontal coordinates (latitude, longitude) and the computed height anomalies. Of the 812 GNSS/levelling points, 428 points are first- and second-order, and 384 points are third-order of the national levelling networks. First-, second- and third-order levelling in Vietnam allows misclosure of|$\ 5\sqrt k $|⁠, 12|$\sqrt k $| and |$25\sqrt k $| mm over a distance of k (km), respectively. Normal height is currently used in the national height system of Vietnam. In Vu et al. (2020), assuming a normal distribution of the residuals, 33 observations were rejected. The 779 good quality points, which are relatively well distributed over the entire country as shown in Fig. 2(b), are used in this study.

4 EFFECTS OF INDIRECT POTENTIAL BIAS AND OMISSION ERROR

Several recent studies have been done on the effects of indirect potential bias on the HSU for different regions in the world, for example in Europe (Gerlach & Rummel 2013) and North America (Amjadiparvar et al. 2016). Its effect on the estimated disturbing potential for Vietnam is also determined. In Vu et al. (2020), the datum offset (⁠|$\delta {H^{\rm{*}}}$|⁠) was estimated at about 69 cm with the global equipotential surface realised by the conventional value W₀ = 62 636 853.4 m² s^–2. We used this offset value in a simulation, in which terrestrial gravity anomalies are biased according to: |$\Delta {g^b} = \ 0.3086\ \delta {H^{\rm{*}}} \approx 0.2$| (mGal). The simulation data are made as follows: the gridpoints derived from the land gravity points (red dots on Fig. 2a) are set to 0.2 (mGal), while the remaining gridpoints (fill-in and tapered grids) are set to 0 (mGal). Here, we assume that the fill-in gridpoints derived from a global model do not have a datum offset. It should be noted that the red dots lying on the surrounding areas, that is Thailand, China, Laos and Cambodia (Fig. 2a) are also set to 0. First, the indirect potential bias term is computed using the Stokes integral (eq. 7) with the original Stokes kernel. The results are shown in Fig. 3(a) and listed in Table 1. The effect computed using the original Stokes kernel is significant and ranges from 0 to 0.386 m² s^–2. Secondly, the indirect potential bias term is computed with different modified-degree kernels. The results are listed in Table 1 and shown for two cases in Figs 3(b) and (c). The effect decreases for higher degrees of truncation. The indirect potential bias term is on average less than 0.1 m² s^–2 (equivalent to 1 cm) for all truncation degrees higher than 60. Therefore, with the aim of determining the height reference system with cm level accuracy, the indirect potential bias term value can be safely neglected if the GGM is used in the computation of the disturbing potential to degree higher than 60. According to the study done by Gerlach & Rummel (2013) for Europe, this value is smaller than 1 cm when a GGM is used in the RCR technique and the residual geoid undulation is computed using a modified Stokes kernel with degree of GGM n > 200. Similarly for North America, Amjadiparvar et al. (2016) recommended using n > 180. The degree of truncation in this study can be smaller because the offset value of our study region is smaller, for example for Alaska it is 148 cm, for Belgium 232 cm, while for Vietnam it is 69 cm. Furthermore, here we use the smoother WG modification for the Stokes’ kernel (Amjadiparvar et al. 2016).

Figure 3.

Indirect bias effect on disturbing potential: (a) N₁ = 0; N₂ = 0, (b) N₁ = 50; N₂ = 60 and (c) N₁ = 60; N₂ = 70.

Figure 4.

Estimated disturbing potential: (a) T_RTM, (b) T_res (N₁ = 220; N₂ = 230) and (c) T^{Stokes_WG_GBVP} = |${{{T}}_{{\rm{GGM}}}} + {{{T}}_{{\rm{RTM}}}} + {{{T}}_{{\rm{res}}}}$|⁠.

Figure 5.

Gravity potential on the GNSS/levelling points: (a) W_i = U_i + T_i and (b) |${( {W_0^{\rm LVD}} )_i}$|⁠.

Table 1.

Indirect bias effect on disturbing potential (m² s^–2).

Truncation degree	Min	Max	Mean	STD
N1 = 0; N2 = 0	0	0.386	0.089	0.097
N1 = 50; N2 = 60	−0.020	0.110	0.015	0.033
N1 = 60; N2 = 70	−0.031	0.088	0.007	0.027
N1 = 100; N2 = 110	−0.055	0.041	−0.007	0.015
N1 = 150; N2 = 160	−0.030	0.030	−0.002	0.009
N1 = 220; N2 = 230	−0.026	0.026	0.002	0.006

A GOCE GGM is used for the computation of the geopotential values, and its omission error is estimated here. Degree/Order (D/O) 260 of GOCE DIR-R5 (Bruinsma et al. 2014) was found to be optimum for this study region (Vu et al. 2019), so we estimate the omission error of the model for that D/O. It is evaluated on 779 GNSS/levelling points, for GOCE DIR-R5 at D/O 260, and after extending with EGM2008 from D/O 261–2190, as well as for the gravimetric-only quasigeoid GEOID_LSC (Vu et al. 2019). Table 2 shows the gravity potential offset value for the LVD with respect to W₀ = 62 636 853.4 m² s^–2 calculated on 779 GNSS/levelling points (see Vu et al. (2020) for details). The potential offset value computed with mixed DIR/EGM up to its full resolution is 6.64 m² s^–2 and that computed with GEOID_LSC is 6.68 m² s^–2, while the GOCE DIR-R5 only and mixed DIR/EGM up to D/O 719 yield 6.12 and 6.45 m² s^–2, respectively. Thus, the effect of the GOCE GGM omission error on the offset value is estimated at 0.5 m² s^–2 (equivalent to 5 cm) and that of the mixed DIR/EGM to D/O 719 to about 0.2 m² s^–2 (equivalent to 2 cm). These results indicate that the contributions of the smaller scales of the gravity field are important when evaluating the potential as well as the datum offset. The offset values computed with DIR/EGM at full resolution and the GEOID_LSC model are very close. This does not necessarily mean that omission errors of both models are equal over the entire country. It may be due to the locations of the GNSS/levelling points, most of which are in lowland areas where the effect of omission error is smaller when compared to mountainous areas (Hirt et al. 2010). To clarify this issue for the study region, we have done tests using different subsets of the GNSS/levelling points according to an elevation threshold (1000 m). With 755 points having elevation <1000 m, the potential offset value is very similar (6.69 and 6.65 m² s^–2 for GEOID_LSC and DIR/EGM at full resolution, respectively). This value is 6.52 and 6.14 m² s^–2 for GEOID_LSC and DIR/EGM, respectively, on 24 points having elevation >1000 m. This means that the effect of the mixed DIR/EGM omission error at full resolution is 0.4 m² s^–2 when the comparison is done in the mountainous areas. The higher frequency information from the terrestrial gravity data in the GEOID_LSC also significantly reduces the STD of the estimated potential offset value (0.86 and 1.53 m² s^–2 for GEOID_LSC and DIR/EGM at full resolution, respectively). It is also clear that DIR/EGM at D/O 719 significantly reduces the STD of the estimated potential offset value when compared to the GOCE DIR-R5 model at D/O 260 (1.69 and 2.29 m² s^–2 for DIR/EGM and GOCE DIR-R5, respectively). Similar to Vu et al. (2019) we use the mixed DIR/EGM model at D/O 719 for the RCR procedure.

Table 2.

Gravity potential offset for LVD with respect to W₀ = 62 636 853.4 m² s^–2 on 779 GNSS/levelling points (m² s^–2).

Model	|$\delta {W^{\rm LVD}}$|	STD
GOCE DIR-R5 (D/O 260)	6.12	2.29
DIR/EGM (D/O 719)	6.45	1.69
DIR/EGM (D/O 2190)	6.64	1.53
GEOID_LSC	6.68	0.86
DIR/EGM (D/O 2190) (<1000 m, 755 points)	6.65	1.50
GEOID_LSC (<1000 m, 755 points)	6.69	0.82
DIR/EGM (D/O 2190) (>1000 m, 24 points)	6.14	2.26
GEOID_LSC (>1000 m, 24 points)	6.52	1.08

5 ZERO-HEIGHT GEOPOTENTIAL FOR LVD ON THE GNSS/LEVELLING POINTS USING THE GBVP APPROACH

The residual disturbing potential (T_res) is computed for 779 GNSS/levelling points using the regular grid of residual gravity anomalies. The GBVP approach based on Molodensky's theory using the Stokes’ integral with the smoother WG modification of the kernel function, called Stokes_WG_GBVP, at degrees N₁ = 220 and N₂ = 230 is used to compute residual disturbing potential. These optimum degrees for the study region were determined in Vu et al. (2019), and as a result the effect of the indirect potential bias terms can be omitted in the estimation disturbing potential. The residual disturbing potential varies from −7.15 to 3.81m² s^–2. The disturbing potential (T^{Stokes_WG_GBVP}) is obtained by restoration of T_GGM (using the mixed DIR/EGM model up to D/O 719) and T_RTM (using the mixed SRTM from degree 720 to 216 000; degree 216 000 equivalent to 3″ resolution of DTM). The results are listed in Table 3. Fig. 4 shows the estimated disturbing potential, which varies from −329.21 to 48.79 m² s^–2.

Table 3.

Statistics of the disturbing potential (m² s^–2).

	Min	Max	Mean	STD
TDIR/EGM (D/O 719)	−323.70	47.90	−147.28	120.11
TRTM	−2.50	2.84	−0.24	0.63
Tres	−7.15	3.81	−0.13	1.34
TStokes_WG_GBVP	−329.21	48.79	−147.65	120.40
Wi	62 620 927.52	62 636 844.69	62 634 717.26	2931.87

From the disturbing potential on 779 GNSS/levelling points (T^{Stokes_WG_GBVP}), we can directly calculate geopotential values on the surface (W_i) by adding the normal potential U_i on the reference ellipsoid (WGS-84 here), and then we can determine the zero-height geopotential level |${( {W_0^{LVD}} )_i}$| based on the normal height (⁠|${H^{\rm{*}}}$|⁠) of the GNSS/levelling points using eq. (4). By simply averaging of |${( {W_0^{\rm LVD}} )_i}$| according to eq. (5) we get the mean zero-height gravity potential (⁠|${{W}}_0^{{\rm{LVD}}}$|⁠) for the LVD. The results are listed in Table 4 and shown in Fig. 5. The zero-height geopotential value varies from 62 636 843.89 to 62 636 849.28 m² s^–2, and the Vietnam Local Vertical Datum (VLVD) is equal to |${{W}}_0^{{\rm{LVD}}}$| = 62 636 846.69 m² s^–2 with STD of 0.82 m² s^–2. There is only a slight improvement in the STD when compared to |${{W}}_0^{{\rm{LVD}}}$| estimated from GNSS/levelling data and the GEOID_LSC model (0.86 m² s^–2 in Table 2). This will be clarified next.

Table 4.

Zero-height gravity potential |${( {W_0^{LVD}} )_i}$| on the GNSS/levelling points based on Stokes_WG_GBVP (m² s^–2 ).

	Min	Max	Mean (⁠|${\rm{W}}_0^{{\rm{LVD}}}$|⁠)	STD
Before detrending	62 636 843.89	62 636 849.28	62 636 846.69	0.82
After detrending	62 636 844.25	62 636 848.74	62 636 846.69	0.70

As mentioned in the Section 2, |${( {W_0^{\rm LVD}} )_i}$| is calculated from the normal heights, so without errors due to inconsistency between GNSS, levelling and gravimetric quasigeoid model for deriving height anomalies. Such errors are always present in the approach using GNSS/levelling and a gravimetric quasigeoid model for the determination of the potential value, for example effect of land subsidence on geometric height anomalies derived from GNSS and levelling data referring to different epoch times (Vu et al. 2020). GNSS/levelling is affected by land subsidence in the Mekong Delta, but this is not visible in Fig. 5(b). However, the systematic error in levelling data always affects the zero-height geopotential |${( {W_0^{LVD}} )_i}$|⁠, which were converted from the geopotential on the surface using physical height derived from levelling. This systematic error was analysed in Vu et al. (2020) based on the differences between geometric height anomalies derived from GNSS/levelling data and those derived from the gravimetric-only quasigeoid model. However, when using GNSS and levelling data to derive the geometric anomalies, besides systematic errors in levelling data, errors in the GNSS measurements, as well as inconsistencies between GNSS and levelling measurements, are always present. This is the motivation to re-analyse the effect of systematic errors in the levelling data based on the estimated gravity potential where only levelling data is used for estimation. This effect can be seen in Fig. 6, which shows the height residuals (e_i). The height residuals are computed based on the mean zero-height gravity potential for the LVD (⁠|${{W}}_0^{{\rm{LVD}}})$| in eq. (5) as follows (Grigoriadis et al. 2014):

$$\begin{eqnarray*} {e_i} = H_i^* - \frac{{{\rm{W}}_0^{{\rm{LVD}}} - {W_i}}}{{{{\bar{\gamma }}_i}}} \end{eqnarray*}$$

(16)

Figure 6.

Height residuals: (a) The north–south direction of linear regression of all points and (b) The east–west direction of linear regression of all points north of 21° latitude.

A trend of 0.006 m deg^–1 in the north–south direction may cause an error of about 8 cm over a distance of 14° (GNSS/levelling data located from 9° to 23° in latitude). The tilt is 0.030 m deg^–1 in the east–west direction in the north of Vietnam. These tilts are slightly smaller than those derived from GNSS/levelling and quasigeoid, but still significant over long distances. The trends in |${( {W_0^{\rm LVD}} )_i}$| due to the systematic errors in the levelling data are shown in Fig. 7. The orange triangles indicate potential residuals that were calculated corresponding to the height residuals.

Figure 7.

Potential residuals, before and after detrending, corresponding to height residuals in Fig. 6: (a) In north–south direction and linear regression of all points and (b) In east–west direction and linear regression of all points in the north (>21° in latitude).

The trends in height residuals are modelled and removed using the third-order polynomial model with two distinct distortion parameters for the two parts, one for the northern part (>17° in latitude) and another for the southern (see Vu et al. 2020). In the models developed by Kotsakis et al. (2012) to estimate the vertical datum offset, systematic errors can be corrected with the geographically correlated terms (i.e. tilts). But the bias in the height data may still remain because the constant component, including the bias term within the height data and the vertical datum offset, is not included in these models. In the model adopted in this study, both the constant and geographically correlated components are incorporated in the equation below:

$$\begin{eqnarray*} {e_i} &=& H_i^* - \frac{{{\rm{W}}_0^{{\rm{LVD}}} - {W_i}}}{{{{\bar{\gamma }}_i}}} - ({a_0} + {a_1}\varphi + {a_2}\lambda + {a_3}{\varphi ^2} + {a_4}\varphi \lambda \nonumber \\ &&+\, {a_5}{\lambda ^2} + {a_6}{\varphi ^3} + {a_7}{\varphi ^2}\lambda + {a_8}\varphi {\lambda ^2} + {a_9}{\lambda ^3}). \end{eqnarray*}$$

(17)

On each GNSS/levelling point we get an equation like eq. (17), and least square (LS) adjustment has been performed with equal weights. After detrending, the corrected normal heights are used to compute the gravity potentials and the results are listed in Table 4. The mean zero-height gravity potential value (⁠|${{W}}_0^{{\rm{LVD}}}$|⁠) does not change, but STD improves significantly (0.70 m² s^–2 after de-trending, decreasing by 0.12 m² s^–2). Fig. 7 shows the improvements of the potential residuals after detrending of the levelling data (green dots). In Vu et al. (2020), the gravity potential of the VLVD was estimated at |${{W}}_0^{{\rm{LVD}}}$| = 62 636 846.81 ± 0.70 m² s^–2 after applying corrections for GNSS/levelling data due to effect of land subsidence in the Mekong Delta, which improved the consistency between GNSS and levelling data significantly. Now we have two estimates for the gravity potential values for VLVD, with a small difference of 0.12 m² s^–2 (equivalent to 1 cm), which is acceptable. This proves that the applied method is reliable for directly calculating the gravity potential value on the Earth's surface. In the following section, the gravity potential value and its variations will be determined based on this approach using GNSS-CORS stations as realisation of the IHRS, to monitor deformation of the height reference system for Vietnam.

6 GEOPOTENTIAL AND ITS VARIATIONS ON THE GNSS-CORS STATIONS

6.1 GNSS and InSAR data

To avoid site motion, geodetic control points, horizontal as well as vertical, were usually set into structures with deep foundations. Such control points were considered as fixed in time. However, even when deeply embedded into the bedrock, these points are susceptible to vertical motion. Thanks to space geodetic techniques, the displacement can be measured. The coordinates of control points can be given at the reference epoch together with their variations in time, allowing continuous monitoring of the deformation of the reference datum. When differences exceed a given threshold, the reference coordinates should be updated. For this purpose, a vertical reference datum, at global as well as national or regional scale, should be realised based on the GNSS-CORS stations. For Vietnam, continuous time-series of heights from 11 permanent GNSS stations over long periods of time were used in Vu et al. (2020) to determine the vertical land motion. From these 11 stations, we selected 3 stations about equally spaced from north to south (named PHUT, QNRS and BACL) to calculate the gravity potential used as the reference points for the vertical datum. Fig. 2(b) shows the locations of these 3 stations (red triangles). Variations in geometric height derived from these GNSS stations can be used to infer change in the normal geopotential. Vertical land motion rates from InSAR reflect topography changes around the GNSS station and can be used to calculate time variations of the disturbing potential. From time variations in normal and disturbing potential, the gravity potential change on the GNSS stations is obtained as follows:

$$\begin{eqnarray*} {\dot{W}_P} = {\dot{U}_P} + {\dot{T}_P}. \end{eqnarray*}$$

(18)

For topography change, we benefited from estimations of the annual average land subsidence rates over 2015–2018 derived from Sentinel-1 imagery time-series covering most of the south of Vietnam, which are available in the frame of the project ‘EMSN062: Assessing changes in ground subsidence rates, Mekong Delta, Vietnam’ (Report of COPERNICUS 2019).

For the purpose of identifying the permanent GNSS stations as the reference points for the modern height system in Vietnam, we first determine the gravity potential value on the surface of these points, and then their variations are calculated.

6.2 Estimated geopotential and its time variations

The Stokes_WG_GBVP approach is used for computation of the residual disturbing potential in 4 × 4⁰ blocks surrounding each GNSS-CORS station. The disturbing potentials are obtained by restoring T_GGM and T_RTM. The geopotential value on the surface is determined by adding the normal potential on the reference ellipsoid. The computation procedure was performed in the same way as for the GNSS/levelling points. The results, listed in Table 5, were calculated using the ellipsoidal height at epoch 2018.0.

Table 5.

Estimated potential on the GNSS stations at epoch 2018.0 (m² s^–2).

Name	Tres	TStokes_WG_GBVP	Ui	Wi
PHUT	−0.20	−264.84	62 636 982.30	62 636 717.46
QNRS	0.44	−89.37	62 636 797.04	62 636 758.62
BACL	−0.87	−38.43	62 636 734.76	62 636 645.40

The ellipsoidal heights time-series derived from continuous measurements on the GNSS-CORS stations are shown in Fig. 8. A notable subsidence rate of −28 mm yr^–1 was observed for the BACL station, located in the Mekong Delta, whereas that of the PHUT stations was only 2 mm yr^–1.

Figure 8.

Ellipsoidal heights time-series on the GNSS-CORS stations: (a) PHUT station, (b) QNRS station and (c) BACL station.

These time-series of GNSS heights reveal signals induced by vertical deformation, but also outliers, especially for the QNRS station. The outliers were determined assuming a normal distribution of the residuals, and three sigma (3σ) rejection led to elimination of 90, 89 and 19 epoch times for the PHUT, QNRS and BACL stations, respectively. The height residuals are computed using the linear model as follows:

$$\begin{eqnarray*} {e^t} = {h^t} - ({a_0} + {a_1}t), \end{eqnarray*}$$

(19)

whereh^t is the ellipsoidal height at epoch time t; a₀ and a₁ are coefficients which are determined using LS adjustment.

The filtered ellipsoidal heights time-series are used to calculate the normal potential time-series, assuming that it does not change on the reference ellipsoid (U₀). From time-series of the normal potential, the movement rates of the normal potential can be determined. The results are shown in Fig. 9. Rates in normal potential of 0.02 ± 0.00(12), 0.04 ± 0.00(15) and 0.28 ± 0.00(52) m²s⁻² yr^–1 were observed for the PHUT, QNRS and BACL stations, respectively. There clearly is a large change in the normal potential in the BACL station, located in the Mekong Delta (Fig. 2b), while it is insignificant for the other two. As the Mekong Delta is known to be deforming for decades, this result is not surprising.

Figure 9.

Normal potential time-series on the GNSS-CORS stations: (a) PHUT station, (b) QNRS station and (c) BACL station.

To determine the time variations in the disturbing potential for the BACL station, we first determine the gravity change due to the changed topography derived from Sentinel-1 InSAR data (Fig. 10a) via forward modelling of residual terrain effects with constant density. Then the Stokes’ integral is applied to this gravity change grid to determine the change in disturbing potential. This approach was used for the determination of the change in height anomalies caused by the 2016 Kaikoura earthquake (McCubbine et al. 2020). The gravity anomaly changes are computed as follows:

$$\begin{eqnarray*} \Delta {g^{\rm change}} = \Delta H\frac{{\partial \gamma }}{{\partial h}} - 2\pi G\rho \Delta H + TC\left( {\Delta H} \right), \end{eqnarray*}$$

(20)

where ΔH is the change in elevation given as the annual average land subsidence derived from Sentinel-1, |$\frac{{\partial \gamma }}{{\partial h}}$| is the linear free-air gravity gradient, |$\rho $| is density of crust, G is the universal gravitational constant and |$TC( {\Delta H} )$| is terrain correction due to the change in topography. The latter two terms on the right-hand side in eq. (20) are the RTM effects due to the change in topography determined by the InSAR data. Eq. (20) can be rewritten as follows:

$$\begin{eqnarray*} \Delta {g^{\rm change}} = \Delta H\frac{{\partial \gamma }}{{\partial h}} + \Delta {g_{\rm RTM}}\left( {\Delta H} \right) \end{eqnarray*}$$

(21)

|$\Delta {g_{\rm RTM}}( {\Delta H} )$| has been computed with the same mixed SRTM model that was used as the detailed DTM. The reference DTM is created by adding the annual average land subsidence rates over 2015–2018 derived from Sentinel-1 (Report of COPERNICUS 2019) into the mixed SRTM model. The grids of subsidence rates are not specified over sea (set to 0). Fig. 10(a) shows the annual average land subsidence rates derived from Sentinel-1. Average subsidence rates range from 1 to 4 cm yr^–1 in built up areas whereas over agricultural areas it is less than 1 cm yr^–1. Fig. 10(b) indicates the forward modelled gravity anomaly changes using eq. (21). The gravity anomaly changes vary from −0.07 to 0.20 mGal.

Figure 10.

(a) Annual average vertical deformation over the 2015–2018 period in the Mekong Delta determined from Sentinel-1 InSAR, project EMSN062 (Report of COPERNICUS 2019), orange triangle is the BACL station and (b) gravity anomaly changes computed from the annual average vertical deformation [min: −0.07 (mGal), max: 0.20 (mGal)].

From the grid of gravity anomaly changes, the Stokes’ integral is used for determining the change in disturbing potential at the BACL station. The change due to significant subsidence in the Mekong Delta is only 0.004 m²s⁻² yr^–1, and is very small when compared to the change in normal potential. Hence, the time variations in geopotential are considered equal to the change in normal potential, that is time variations in the disturbing potential can be omitted in the estimation of geopotential change. According to the recommendation in Sánchez et al. (2016), an update of W₀ is required if the cumulative change reaches 2 m²s⁻². This value will be reached within 8 yr at the BACL station at the current rate.

7 CONCLUSION

The effects of the GGM omission error, the indirect potential bias term and the systematic levelling error on the estimation of the gravity potential were investigated. The results indicated that the indirect potential bias term is less than 0.1 m² s^–2 (equivalent to 1 cm) if the degree of truncation is higher than 60. The objective of this study is to determine the height reference system with cm level accuracy, and so the indirect potential bias term can be safely neglected when using GOCE DIR-R5 up to D/O 260 in combination with the smoothing WG modification, which fully removed low harmonics up to degree N₁ = 220, then linearly tapered to degree N₂ = 230. The effect of the GOCE DIR-R5 (at D/O 260) omission error on the offset value, based on the GNSS/levelling data, was estimated by extending with EGM2008 from D/O 261–2190. It is at 0.5 m² s^–2 (equivalent to 5 cm). This proves that the GOCE omission errors should be taken into account in determination of the geopotential value for this region. The remaining quasigeoid signal of the mixed DIR/EGM above D/O 2190 was also estimated by the GEOID_LSC model. This signal is still significant (0.4 m² s^–2), especially in the mountainous areas.

After removing the trends in levelling data due to the systematic cumulative errors, the mean zero-height gravity potential of the VLVD was estimated equal to |${{W}}_0^{{\rm{LVD}}}$| = 62 636 846.69 m² s^–2 with STD of 0.70 m² s^–2 based on the GBVP approach. This value is very similar to the result in Vu et al. (2020) and the difference of 0.12 m² s^–2 (equivalent to 1 cm) is within the limits of the study objective. This proves that the applied method is reliable for directly calculating the gravity potential value on the Earth's surface. Hence, the GBVP approach was used for determination of the geopotential on the surface of three GNSS-CORS stations in Vietnam. Based on time-series of the vertical component derived from the GNSS observations as well as InSAR data, time variations of geopotential were also estimated on these permanent GNSS stations. The purpose is to monitor deformation of the vertical datum. The geopotentials on the surface as well as the velocities of the three GNSS stations were estimated at epoch 2018.0 equal to W^PHUT = 62 636 717.46 + 0.02 m² s^–2, W^QNRS = 62 636 758.62 + 0.04 m² s^–2 and W^BACL = 62 636 645.40 + 0.28 m² s^–2. The cumulative change will reach 2 m² s^–2 within 8 yr at the BACL station at the current rate. Sánchez et al. (2016) recommended an update of W₀ if the cumulative change reaches this threshold value.

This study confirmed again that the geopotential values need to be monitored and determined with the time-dependent component on the permanent GNSS stations in order to plan updates. It should be noted that for the GNSS-CORS stations the zero-height geopotential value should also be known. We recommend that the levelling height on these stations also needs to be determined for conversion of the gravity potential on the Earth's surface to the reference of a LVD.

AVAILABILITY OF DATA AND MATERIALS

The data of this study: gravity anomaly and height anomaly grids are available via BGI website (http://bgi.obs-mip.fr/).

ACKNOWLEDGEMENTS

We are very grateful to the VDSM, the Institute of Geophysics (IGP)—Vietnam Academy of Science and Technology (VAST), the Vietnam Institute of Geodesy and Cartography (VIGAC) and the Bureau Gravimétrique International (BGI) that supplied data for this study. The SRTM3arc_v4.1 is available via http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp. The SRTM15arc_plus model is available via https://topex.ucsd.edu/pub/srtm15_plus/. All used GGMs in this paper are available via ICGEM website (Ince et al. 2019) http://icgem.gfz-potsdam.de/tom_longtime. The DTU15 gravity is available via https://ftp.space.dtu.dk/pub/. The processing result of Sentinel-1 is available from project EMSN062 via https://emergency.copernicus.eu/mapping/list-of-components/EMSN062. In this paper, we used the GMT for producing some of figures. We are grateful to the two anonymous reviewers who helped us at improving our manuscript.

Author contributions: Dinh Toan Vu, Sean Bruinsma and Sylvain Bonvalot designed the research. Dinh Toan Vu developed the software and performed all the data processing. Dinh Toan Vu drafted the manuscript. All authors analysed and discussed the preliminary results. Sean Bruinsma and Sylvain Bonvalot supported the observations. All authors read and approved the final manuscript.

FUNDING

Dinh Toan VU receives funding of the Vietnamese Government's 911 project and Université Toulouse III-Paul Sabatier (UPS)-GET. This study is also supported by CNES, CNRS and IRD.