On the Consistency between a Classical Definition of the Geoid-to-Quasigeoid Separation and Helmert Orthometric Heights

It is acknowledged that a classical definition of the geoid-to-quasigeoid separation as a function of the simple planar Bouguer gravity anomaly is compatible with Helmert’s definition of orthometric heights. According to Helmert, the mean actual gravity along the plumbline between the geoid and the topographic surface in the definition of orthometric height is computed approximately from the measured surface gravity by applying the Poincaré-Prey gravity reduction. This study provides theoretical proof and numerical evidence that this assumption is valid. We demonstrate that differences between the normal and (Helmert) orthometric corrections are equivalent to the geoid-to-quasigeoid separation differences computed for individual levelling segments. According to our theoretical estimates, maximum differences between these 2 quantities should be less than ±1 mm. By analogy, differences between the Molodensky normal and Helmert orthometric heights at levelling benchmarks should be equivalent to the geoid-to-quasigeoid separation computed from the Bouguer gravity data. Both theoretical findings are inspected numerically by using levelling and gravity data along selected closed levelling loops of the vertical control network in Hong Kong. Results show that values of the geoid-to-quasigeoid separation at levelling benchmarks differ less than ±0.1 mm from differences between the normal and orthometric corrections. Relatively large differences (slightly exceeding 2 mm) between values of the geoid-to-quasigeoid separation and differences between the normal and (Helmert) orthometric heights at levelling benchmarks are explained by errors in levelling measurements rather than by inconsistencies in computed values of the geoid-to-quasigeoid separation and (Helmert) orthometric correction.


Introduction
Orthometric and normal heights are the two most commonly used types of heights for a practical realization of geodetic vertical controls by physically establishing levelling benchmarks. Both types of heights are determined from precise levelling and gravity measurements (along levelling lines). Since the realization of levelling networks and their maintenance is extremely costly, geoid and quasigeoid models have been used in many countries around the world for a realization of geodetic vertical datums (either solely or together with a levelling network). Detailed geoid/quasigeoid models are also indispensable for the conversion of geodetic (geometric) heights measured by the Global Navigation Satellite Systems (GNSS) techniques (such as GPS) for officially adopted orthometric/normal heights. In essence, concepts of defining orthometric and normal heights (or equivalently geoid or quasigeoid heights) depend on a treatment of the topographic density. Whereas the topographic density is taken into consideration in the definition of orthometric heights, Molodensky [1,2] completely disregarded the topographic density in the definition of normal heights (see also [3]). Normal heights can then directly be determined from levelling and gravity measurements. His main argument was that knowledge of the actual expressions for the geoid-to-quasigeoid separation difference and its relationship with the normal and orthometric correction difference for a levelling segment are derived.

Normal and Orthometric Heights
The orthometric height H O is defined by (e.g., [32] and Equations (4)- (21)) where C is the geopotential number of a point at the topographic surface, and g is the mean actual gravity along the plumbline within the topography (see Figure 1).
OR PEER REVIEW 4 of 24

Difference between the Normal and Orthometric Corrections
Heights H obtained from levelling measurements (Equation (4)) are typically converted to either orthometric heights O H (Equation (1)) or normal heights N H (Equation (2)) by applying the orthometric or normal corrections to levelled height differences i H Δ , respectively. It is worth noting that heights H could eventually be first converted to normal-orthometric heights (Equation (3)) by using normal gravity values computed along levelling lines, and consequently to the normal heights by applying the cumulative normal to normal-orthometric height correction [51][52][53]. This two-step numerical scheme is obviously not beneficial if we are not particularly interested in determining this type of height. The orthometric correction , 1 i i OC + of a levelling segment between two benchmarks i and i + 1 is defined by (e.g., [32] (Equations (4)-(33)) The normal height H N was defined by [1,2] in the following form where the mean normal gravity γ along the ellipsoidal normal between the reference ellipsoid and the telluroid (i.e., the surface on which the normal potential is equal to the actual potential at the topographic surface) is evaluated according to the Somigliana-Pizzetti's theory of the normal gravity field [49,50].
The geopotential number C in definitions of orthometric and normal heights (Equations (1) and (2)) is practically computed from measured levelling height differences ∆H i and observed gravity values g i along levelling lines, i.e., C = ∑ i g i ∆H i . When normal gravity values are used instead of observed gravity values along levelling lines, the vertical geodetic datum is realized in the system of normal-orthometric heights H N−O defined by (e.g., [17]) where the normal geopotential number C N is computed from measured levelling height differences ∆H i and normal gravity values γ i along levelling lines, i.e., Heights of levelling benchmarks in some countries and territories (such as Hong Kong) have been determined only from levelling measurements, involving neither the actual nor normal gravity information. In this case, heights H of levelling benchmarks are directly computed from measured levelling height differences ∆H i so that

Difference between the Normal and Orthometric Corrections
Heights H obtained from levelling measurements (Equation (4)) are typically converted to either orthometric heights H O (Equation (1)) or normal heights H N (Equation (2)) by applying the orthometric or normal corrections to levelled height differences ∆H i , respectively. It is worth noting that heights H could eventually be first converted to normalorthometric heights (Equation (3)) by using normal gravity values computed along levelling lines, and consequently to the normal heights by applying the cumulative normal to normalorthometric height correction [51][52][53]. This two-step numerical scheme is obviously not beneficial if we are not particularly interested in determining this type of height.
The orthometric correction OC i,i+1 of a levelling segment between two benchmarks i and i + 1 is defined by (e.g., [32] (Equations (4)- (33)) where H i and H i+1 are heights of levelling benchmarks i and i + 1, respectively, g i and g i+1 are the corresponding mean gravity values (as defined in Equation (1)), and δH k are levelled height differences (at levelling setups k between two benchmarks i and i + 1 of a levelling segment, i.e., ∆H i,i+1 = i+1 ∑ k=i δH k ). The normal gravity γ 0 at the reference ellipsoid in Equation (5) is a constant value, meaning that it is computed for the same geodetic latitude, for instance ϕ = 45 • . The normal correction NC i,i+1 is given by (e.g., [32] (Equations (4)- (45)) where γ i and γ i+1 are the mean normal gravity values (as defined in Equation (2)). As seen from the comparison of Equations (5) and (6), the mean normal gravity values γ i and γ i+1 are used in the definition of the normal correction NC i,i+1 instead of the mean gravity values g i and g i+1 in the orthometric correction OC i,i+1 . From Equations (5) and (6), the normal and orthometric correction difference (i.e., the difference between the normal and orthometric corrections) of an individual levelling segment between two benchmarks i and i + 1 is found to be Equation (7) defines the normal and orthometric correction difference rigorously as a function of heights H i and H i+1 of benchmarks determined from levelling measurements. In addition, the functional relation involves differences between the mean actual and normal gravity values g i − γ i and g i+1 − γ i+1 (at locations of benchmarks i and i + 1).

Geoid-to-Quasigeoid Separation Difference
The geoid-to-quasigeoid separation χ can be defined as a difference between the normal and orthometric heights so that By analogy with Equation (8), the geoid-to-quasigeoid separation difference ∆χ i,i+1 for a levelling segment (between two benchmarks i and i + 1) can be defined in the following form: The orthometric height difference ∆H O i,i+1 (between two levelling benchmarks i and i + 1) is obtained by applying the orthometric correction OC i,i+1 to the levelled height difference ∆H i,i+1 so that The application of the normal correction NC i,i+1 to the levelled height difference ∆H i,i+1 yields the normal height difference ∆H N i,i+1 . Hence, By combining Equations (9)-(11), the following relation is found: As seen in Equation (12), the difference between the normal and orthometric height differences ∆H N i,i+1 − ∆H O i,i+1 directly equals the difference between the normal and orthometric corrections NC i,i+1 − OC i,i+1 .
Substitution from Equation (7) to Equation (12) yields By combining Equations (9) and (13), the following expression is obtained: Equation (14) defines the geoid-to-quasigeoid separation difference ∆χ i,i+1 for a levelling segment (between two benchmarks i and i + 1) in terms of the normal and orthometric correction difference NC i,i+1 − OC i,i+1 that was introduced in Equation (7).

Approximate Definition of the Geoid-to-Quasigeoid Separation
As stated in Section 1, orthometric heights in all countries are defined according to Helmert's theory [4,5]. By analogy, the geoid-to-quasigeoid separation χ is defined as a difference between the Molodensky normal height H N and the Helmert orthometric height H O . Despite this definition being well known and readily found in geodetic literature, Sensors 2023, 23, 5185 6 of 21 e.g., [32], it is worth recapitulating derivations of its approximate form χ = H N − H O from its rigorous definition χ = H N − H O (given in Equation (8)) in order to better understand the adopted assumptions concerning the accuracy.
Inserting from Equations (1) and (2) to Equation (8), the geoid-to-quasigeoid separation χ becomes As seen in Equation (15), the geoid-to-quasigeoid separation is defined rigorously as a function of the difference between the mean actual and normal gravity values. The mean normal gravity γ in the denominator of Equation (15) can be replaced by the normal gravity value γ 0 at the reference ellipsoid. If the mean normal gravity γ is defined as a function of the normal gravity value γ 0 and the normal (linear) gravity gradient ∂γ/∂h, i.e., γ ∼ = γ 0 + ∂γ/∂h (H/2), the difference between using γ 0 instead of γ introduces the following approximation error: where γ ∼ = GM/R 2 and, consequently, ∂γ/∂h ∼ = −2GM/R 3 are defined in terms of the Earth's mean radius R = 6371 × 10 3 m and the geocentric gravitational constant GM = 3.986 × 10 14 m 3 s −2 (i.e., the product of the Newton gravitational constant G and the total mass of the Earth M). For maximum elevations of H ≈ 9 km (in the Himalayas), the approximation error could reach maximum of ε γ−γ 0 ≈ 0.014 m s −2 (Equation (16)). In Hong Kong, these errors are obviously much smaller. The error ε γ−γ 0 propagates into the error in values of the geoid-to-quasigeoid separation ε χ γ as follows: Note that a standard error propagation was not applied to derive the expression in Equation (17), because only the mean normal gravity γ in the denominator of Equation (15) is approximated. Moreover, nonlinear gravity changes are disregarded, having no impact on the error analysis.
Substitution from Equation (16) to Equation (17) yields For maximum values of the geoid-to-quasigeoid separation χ within ±5 m, the error ε χ γ is less than ±7 mm in mountainous regions with extreme elevations (particularly in the Himalayas, Tibet, and the Andes). Elsewhere, this error is typically less than ±1 mm. The geoid-to-quasigeoid separation χ in Equation (15) can then be defined as follows: Approximations adopted in Helmert's definition of orthometric heights are applied to the difference between the mean actual and normal gravity values g − γ in Equation (19). As already mentioned in Equation (16), the mean normal gravity γ is described as follows: Note that in Equation (20), the normal height H N is used in the definition of the mean normal gravity according to Molodensky [1][2][3].
The mean actual gravity gradient in Equation (19) is further approximated the Poincaré-Prey gravity gradient so that The actual mean gravity in Equation (21) is first described as a function of the surface gravity g and the actual gravity gradient ∂g/∂H. The actual gravity gradient is then approximated by the Poincaré-Prey gravity gradient, ∂g/∂H ≈ ∂γ/∂h + 4 π G ρ T , which is defined as the normal linear gravity gradient ∂γ/∂h and the term 4 π G ρ T (i.e., the Poisson equation). An average upper continental crustal density of 2670 kg m −3 [54] is typically adopted as a topographic density ρ T in geodetic and geophysical applications. It is worth noting that the actual topographic density could vary substantially with respect to the average topographic density of 2670 kg m −3 . In the Hong Kong territories, for instance, the igneous and sedimentary rocks of lower densities represent most of the geological setting. Nevertheless, the density value of 2670 kg m −3 is, until now, exclusively used to define the (Helmert) orthometric heights in countries around the world where this type of height is adopted officially for a realization of the geodetic vertical control. Consequently, this density value is used to compute the geoid-to-quasigeoid separation.
Combining Equations (20) and (21), the following expression is found: Assuming that ∂γ ∂h Equation (22) further simplifies to Note that approximations applied in Equations (22)-(24) do not affect the accuracy. Inserting from Equation (24) to Equation (19), the geoid-to-quasigeoid separation is then obtained in the following form: From a definition of the free-air gravity anomaly ∆g FA , i.e., the expression in Equation (25) becomes The simple planar Bouguer gravity anomaly ∆g SPB is computed from the free-air gravity anomaly ∆g FA by applying the Bouguer gravity reduction. Hence, Substitution from Equation (28) to Equation (27) yields (e.g., [32] (Equations (8)-(103)) Equation (29) defines the geoid-to-quasigeoid separation χ approximately as a function of the simple planar Bouguer gravity anomaly ∆g SPB and the normal height H N . Note that heights H (obtained from levelling measurements) can be considered instead of H N and H O in Equations (29) and (28) without affecting the accuracy of χ .

Approximate Definition of the Geoid-to-Quasigeoid Separation Difference
In Equation (14), the relation between the geoid-to-quasigeoid separation difference and the normal and orthometric correction difference (for a levelling segment between two benchmarks i and i + 1) was rigorously formulated. To find the corresponding approximate definition of the geoid-to-quasigeoid separation difference, the rigorous definition of the geoid-to-quasigeoid separation was first used to derive the rigorous definition of the geoid-to-quasigeoid separation difference. From Equation (15), we can write Combining expressions defining the geoid-to-quasigeoid separation difference in Equations (14) and (30), the following relation is obtained Substitution from Equations (30) and (14) yields To estimate approximation errors caused by using heights H (obtained from levelling measurements) instead of orthometric heights H O in Equation (32), the error analysis is, firstly, written as follows: where Substitution from Equations (34) and (35) back to Equation (33) then yields The error ε ∆χ−NC+OC depends mainly on differences between the mean actual and normal gravity values at levelling benchmarks. Considering even very large differences g − γ in Equation (36) of the order of several hundreds of milligals, the error ε ∆χ−NC+OC is completely negligible (less than 0.1 mm). Consequently, it could be concluded that rigorous definitions of the geoid-to-quasigeoid separation difference according to formulas given in Equations (14) and (30) are equivalent.
Finally, the approximate definition of the geoid-to-quasigeoid separation in Equation (29) is used to introduce the approximate definition of the geoid-to-quasigeoid separation difference ∆χ i,i+1 in the following form Again, normal heights in Equation (37) can be disregarded. Instead, heights H obtained from levelled height differences can be used. Such an assumption introduces errors in values of χ typically less than ±1 mm that correspond to even smaller errors in values of ∆χ i,i+1 so that Consequently, Rigorous and approximate relations between the normal and orthometric correction difference and the geoid-to-quasigeoid separation difference were derived and presented above. The rigorous definition in Equation (30) was described by means of the difference between the actual mean and normal gravity values. The corresponding approximate relation in Equation (39) is described as a function of the simple planar Bouguer gravity anomaly values. If the Poincaré-Prey gravity gradient closely approximates the actual vertical gravity gradient inside the topography, both definitions should provide results that differ less than ±1 mm. In other words, the relation NC i,i+1 − OC i,i+1 ∼ = ∆χ i,i+1 in Equation (39) should be accurate enough to be applicable for a practical realization of Helmert orthometric heights H O and a conversion between the Helmert orthometric and Molodensky normal heights (i.e., χ ) by using the formula in Equation (29). This theoretical aspect is numerically investigated next.

Numerical Procedures
The accuracy of the geoid-to-quasigeoid separation differences was assessed at a vertical geodetic control in Hong Kong, practically realized by the Vertical Control Network 2022 (VCN2022). Since gravity values along levelling lines were not measured directly, detailed terrestrial and marine gravity measurements were used to interpolate gravity values along levelling lines [55]. Interpolated gravity values (at levelling benchmarks) were then used to compute the orthometric and normal corrections to measured levelling height differences, and the entire levelling network was readjusted. The newly determined normal and orthometric heights of levelling benchmarks were presented as the VCN2022 solution. The adjustment of the orthometric and normal levelling networks of the VCN2022 attained a −0.2 and 2.0 mm misclosure, respectively.
To keep the presentation simple but still instructive, we conducted the numerical analysis only along four closed levelling loops of the VCN2022 levelling network which were characterized by the largest topographic elevation changes in Hong Kong. The location of the selected VCN2022 levelling sections is illustrated in Figure 2. mined normal and orthometric heights of levelling benchmarks were presented as the VCN2022 solution. The adjustment of the orthometric and normal levelling networks of the VCN2022 attained a −0.2 and 2.0 mm misclosure, respectively.
To keep the presentation simple but still instructive, we conducted the numerical analysis only along four closed levelling loops of the VCN2022 levelling network which were characterized by the largest topographic elevation changes in Hong Kong. The location of the selected VCN2022 levelling sections is illustrated in Figure 2.

Gravity Data Interpolation
Since gravity data from mainland China are not publicly available, [55] applied a simple gravity data interpolation instead of more refined methods based on computing the complete spherical Bouguer gravity anomalies [56] and their downward continuation (by solving the inverse of the Poisson integral equation). First, they used the measured free-air gravity anomalies ∆g FA to compute the simple planar Bouguer anomalies ∆g SPB by applying the Bouguer gravity reduction. They then used the simple planar Bouguer gravity anomalies (at gravity sites) to interpolate the corresponding values at levelling benchmarks by applying the inverse distance weighted mean. This method was selected based on testing a performance of various interpolation techniques, particularly by applying the kriging, natural neighbor, least-squares collocation, nearest neighbor, and radial basis functions for the gravity data interpolation. According to their results, maximum differences in interpolated gravity values (i.e., the simple planar Bouguer gravity anomalies) from these methods at levelling benchmarks were within ±5 mGal. Such gravity differences correspond to differences in computed values of the normal and orthometric corrections less than ±1 mm (cf. [55]). Finally, they converted the interpolated simple planar Bouguer gravity anomalies to the free-air gravity anomalies at levelling benchmarks. A gravity data interpolation by applying only the Bouguer gravity reduction (while disregarding the terrain gravity correction) is obviously less accurate. A discussion of this aspect is postponed until Section 5. Nevertheless, it is worth already clarifying here that this factor does not affect the numerical findings in this study.

Orthometric and Normal Heights
Nsiah Ababio and Tenzer [55] used the free-air gravity anomalies ∆g FA at levelling benchmarks to compute the orthometric and normal corrections to measured levelling height differences. For 2 consecutive levelling benchmarks i and i + 1 of a levelling segment, the orthometric correction OC i,i+1 was computed according to [25] using the following equation: By analogy with Equation (40), the normal correction NC i,i+1 was computed from Adopting Helmert's definition of orthometric heights, the mean gravity in Equation (40) was computed from where the surface gravity g i was computed from the free-air gravity anomaly ∆g FA as follows: Substitution from Equation (43) to Equation (42) yields Equations (21) and (44) are equivalent. The last term on the right-hand side of Equation (44) is computed for the topographic density value of 2670 kg m −3 . To check the correctness of both results, the orthometric and normal corrections were also computed according to Equations (5) and (6), finding that both results are equal.
The orthometric height differences ∆H O i,i+1 and the normal height differences ∆H N i,i+1 between levelling benchmarks i and i + 1 were computed according to Equations (10) and (11) by applying the orthometric OC i,i+1 and normal NC i,i+1 corrections to measured levelling height differences ∆H i,i+1 , respectively. The newly determined orthometric ∆H O i,i+1 and normal ∆H N i,i+1 height differences were then used to readjust the entire levelling network and to compute the normal and orthometric heights of levelling benchmarks.

The Orthometric and Normal Correction Differences and the Geoid-to-Quasigeoid Separation Differences
Values of the orthometric corrections OC i,i+1 and the normal corrections NC i,i+1 prepared by [55] were used to compute the orthometric and normal correction differences OC i,i+1 − NC i,i+1 along four selected closed levelling loops. The geoid-to-quasigeoid separation differences ∆χ i,i+1 = χ i+1 − χ i were computed according to Equation (39) from the Bouguer gravity anomaly values ∆g SPB i+1 and ∆g SPB i . For completeness, the geoid-to-quasigeoid separation values χ i were computed according to Equation (29) and compared with the differences between the Molodensky normal heights and the Helmert orthometric heights at levelling benchmarks (obtained after the adjustment of the whole levelling network).

Results
The orthometric and normal corrections and their cumulative values along four selected closed levelling loops are plotted in Figures 3 and 4 Tables 1 and 2), respectively, where the topographic relief and height differences (between individual levelling segments) are also shown. Values of the orthometric correction are mostly within ±3 mm, and values of the normal correction vary largely within ±2 mm (Table 1).  Both corrections reach maximum (absolute) values at levelling sections characterized by the largest elevation changes. Cumulative values of the orthometric correction closely mimic a topographic relief, while this trend in cumulative values of the normal correction is less pronounced (see Figure 4, upper panels) so that their differences (see Figure 4, middle panels) are mostly attributed to the cumulative orthometric correction. Both corrections reach maximum (absolute) values at levelling sections characterized by the largest elevation changes. Cumulative values of the orthometric correction closely mimic a topographic relief, while this trend in cumulative values of the normal correction is less pronounced (see Figure 4, upper panels) so that their differences (see Figure 4, middle panels) are mostly attributed to the cumulative orthometric correction.

(with statistical summaries in
As discussed in Section 2, the orthometric and normal correction differences should be equal (or very similar) with the geoid-to-quasigeoid separation differences computed for individual levelling segments. This aspect was inspected in Figure 5 (with the statistical summary in Table 3), where values of the orthometric and normal correction differences were plotted and compared with values of the geoid-to-quasigeoid separation differences ∆χ i,i+1 . As can be seen, the geoid-to-quasigeoid separation differences fully agree with the orthometric and normal correction differences.  As discussed in Section 2, the orthometric and normal correction differences should be equal (or very similar) with the geoid-to-quasigeoid separation differences computed for individual levelling segments. This aspect was inspected in Figure 5 (with the statistical summary in Table 3), where values of the orthometric and normal correction differences were plotted and compared with values of the geoid-to-quasigeoid separation differences As can be seen, the geoid-to-quasigeoid separation differences fully agree with the orthometric and normal correction differences.  For completeness, we compared values of the geoid-to-quasigeoid separation χ i at levelling benchmarks with corresponding values computed cumulatively from the geoidto-quasigeoid separation differences ∆χ i,i+1 . The results are plotted in Figure 6 (with the statistical summary in Table 4).
As seen, cumulative and pointwise values (at levelling benchmarks) agree. Finally, we compared values of the geoid-to-quasigeoid separation χ i with the differences between the Molodensky normal heights and the Helmert orthometric heights at levelling benchmarks obtained after the readjustment of the whole network. The results are plotted in Figure 7, with the statistical summary of the results in Table 5. As seen, values of the geoid-toquasigeoid separation differ from differences between the Molodensky normal heights and the Helmert orthometric heights. The largest differences, exceeding even 2 mm, are seen along the L12/HK levelling loop in the Lantau Island. The existence of these relatively large differences is discussed in the next section.    As seen, cumulative and pointwise values (at levelling benchmarks) agree. Finally, we compared values of the geoid-to-quasigeoid separation i χ′ with the differences between the Molodensky normal heights and the Helmert orthometric heights at levelling benchmarks obtained after the readjustment of the whole network. The results are plotted in Figure 7, with the statistical summary of the results in Table 5. As seen, values of the geoid-to-quasigeoid separation differ from differences between the Molodensky normal heights and the Helmert orthometric heights. The largest differences, exceeding even 2 mm, are seen along the L12/HK levelling loop in the Lantau Island. The existence of these

Discussion
In Hemert's definition of orthometric heights, the Poincaré-Prey gravity gradient was adopted to approximate the actual gravity gradient. The mean gravity within the topography was then approximated by the surface gravity continuing downward to a midpoint by using the normal gravity gradient and by applying the Poisson equation to take the topography into consideration. To convert the Molodensky normal heights to the Helmert orthometric heights and vice versa, the geoid-to-quasigeoid separation was defined as a function of the simple planar Bouguer gravity anomaly . It is thus expected that the expressions for computing the normal and orthometric correction differences and the geoid-to-quasigeoid separation differences provide the same (or very similar) results.

Discussion
In Hemert's definition of orthometric heights, the Poincaré-Prey gravity gradient was adopted to approximate the actual gravity gradient. The mean gravity within the topography was then approximated by the surface gravity continuing downward to a midpoint by using the normal gravity gradient and by applying the Poisson equation to take the topography into consideration. To convert the Molodensky normal heights to the Helmert orthometric heights and vice versa, the geoid-to-quasigeoid separation was defined as a function of the simple planar Bouguer gravity anomaly ∆g SPB based on adopting approximations equivalent to those used in Helmert's definition of orthometric heights and further rearranging the expression in terms of g − γ to its final form described as a function of ∆g SPB . It is thus expected that the expressions for computing the normal and orthometric correction differences and the geoid-to-quasigeoid separation differences provide the same (or very similar) results.
As demonstrated in the numerical examples, this assumption is correct. Both computed differences are almost the same with differences reaching less than ±0.1 mm. It is worth noting that [55] inspected the reliability of intermediate numerical steps involved to compute the normal and orthometric correction differences and the geoid-to-quasigeoid separation differences, assuring that findings presented in this study are valid. Particularly, they demonstrated that the application of different gravity interpolation techniques does not affect the accuracy. Obviously, the selection of a gravity interpolation technique in mountainous regions with much higher elevation changes requires a careful analysis.
These findings ascertain that the computation of the geoid-to-quasigeoid separation from the Bouguer gravity data is fully compatible with Helmert's definition of orthometric heights. The expression for computing the geoid-to-quasigeoid separation from the simple planar Bouguer gravity anomaly provides the result that is equal to the geoid-to-quasigeoid separation differences.
Finally, we demonstrated that the geoid-to-quasigeoid separation differs from differences between the Molodensky normal and Helmert orthometric heights obtained after the levelling network readjustment (carried out individually for the normal and orthometric height differences), particularly along the L12/HK closed levelling loop. According to our results, the differences between the geoid-to-quasigeoid separation and differences between the Molodensky normal and Helmert orthometric heights there exceed even 2 mm. This inconsistency is explained mainly by the propagation of errors in measured levelling height differences. This was confirmed by the analysis of the adjusted levelled height differences. The results of this analysis revealed that levelled height differences along the L12/HK closed levelling loop located in the Lantau Island are systematically affected by errors of levelling measurements conducted along the bridge, which connect the island with the rest of the territories. Elsewhere, these differences are much smaller.

Summary and Concluding Remarks
We have demonstrated that the geoid-to-quasigeoid separation defined as a function of the simple planar Bouguer gravity anomaly is fully compatible with Helmert's definition of orthometric heights. Since both definitions involve the same assumptions regarding the approximation of the actual gravity gradient by the Poincaré-Prey gravity reduction, the computation of the geoid-to-quasigeoid separation in terms of the simple planar Bouguer gravity anomaly in Equation (29) introduces errors that should not exceed more than ±1 mm (except for extremely large topographic elevations particularly in the Himalayas, Tibet, and the Andes). Obviously, the mean gravity values computed by means of applying the Poincaré-Prey gravity reduction might still be quite inaccurate. Nevertheless, the approximately-computed geoid-to-quasigeoid separation (from the Bouguer gravity data) is consistent with the Helmert's definition of orthometric heights that has been, until now, exclusively used for a practical realization of geodetic vertical controls in countries where orthometric heights are officially adopted.

Data Availability Statement:
The data that support the findings of this study are available on request from the corresponding author.