Geoid model validation and topographic bias

: Recently a number of geoid campaigns were per formed to verify di ﬀ erent types of geoid and quasigeoid mod eling techniques. Typically, GNSS - leveling was employed as an independent method, but in some cases zenith camera astronomic de ﬂ ection data were also used in astrogeodetic determinations of the geoid and/or quasigeoid. However, due to the uncertainty in the topographic density distribution data ( and thereby in orthometric heights ) , we conclude that neither GNSS - leveling nor astrogeodetic techniques can reliably verify di ﬀ erences between gravimetric geoid models at several centimeter levels in rough mountainous regions. This is because much the same topographic data are used both in the gravimetric geoid models and in their veri ﬁ cations by geometric and/or astrogeodetic geoid models. On the contrary, this is not a problem in verifying gravimetric quasigeoid models, as they are independent of the topographic density distribution, and so is the related normal height used in GNSS - leveling.


Introduction
During the last two decades several geoid campaigns have been performed, resulting in more accurate geoid and quasigeoid models than ever before. Frequently the effect of topographic mass density variations is significant but smaller than other geoid modeling uncertainties that one typically encounters in geoid validations. However, results from the geoid test networks in Auvergne, France (e.g., Yildiz et al. 2012), and more recently along a polygon in the Colorado Rocky Mountains with a rough topography reaching 4,000 m of elevation (e.g., Willberg et al. 2020, Wang et al. 2021), suggest that geoid or quasigeoid accuracies of the order of a few centimeters are possible. However, it is well known that the quality of a geoid model depends on the available knowledge of the topographic density distribution, while this is not needed for determining quasigeoid models. There are various methods used in studying the quality of gravimetric geoid models, such as determining standard errors in least squares modeling and comparison with more or less independent techniques. Here we study the effect of the erroneous data of topographic density distribution in the validation of geoid models and scrutinize some of the strategies.

Propagation of erroneous topographic data
Let us neglect the topographic error caused by erroneous topographic height and consider only the effect of uncertainty in the mean density ρ 0 of the Bouguer shell of thickness H. If ρ 0 is correct, the topographic potential of the Bouguer shell at a point P at radius r P = R + H P would be given by where μ 0 = Gρ 0 = gravitational constant × density, R and R s = R + H are the inner and exterior radii of the shell, respectively, σ is the unit sphere and = + − l r r r r ψ 2 cos P P P 2 2 with ψ representing the geocentric angle between computational and integration points. The total topographic correction of order H 2 , e.g., in applying the remove-compute-restore or KTH techniques, becomes (cf. Sjöberg 2007) where γ 0 is normal gravity at the reference ellipsoid. If ρ 0 = 2,670 kg/m 3 , this correction reaches the magnitudes of 11, 44 and 280 cm for H = 1, 2 and 5 km, respectively.
Considering an error dμ in the mean density, it propagates to the remaining geoid model error  By assuming an error of 10% (which may well be exceeded) the remaining geoid biases will be of the orders of 11, 44 and 280 mm, respectively. Hence, there is no doubt that such errors are most significant when verifying and comparing various precise geoid solutions in rough topography.

Use of least squares modeling
In the KTH method (Sjöberg and Bagherbandi 2017, Sect. 6.2.2), frequently named least squares modification of Stokes' formula, different types of data (such as a global gravitational model and terrestrial gravity data) are combined in a least squares procedure that also provides standard errors of the estimated geoid heights. For comparison, standard errors can also be computed from theoretical error estimates obtained from a priori error and signal degree variances that are used in the adjustment (without actual gravity data). Comparing a priori and posteriori solution standard errors yields an indication of the fitness of the used theory, degree variances and data, but the error sources of the erroneous topographic height and density models will not be controlled in this way.
One recent version of the UNB technique for geoid determination (Foroughi et al. 2019) employs a discrete least squares method for solving an overdetermined version of Poisson's equation in downward continuation of gravity anomalies from Earth's surface to the Bjerhammar sphere (cf. Sjöberg 1975, Bjerhammar 1975, Sjöberg and Bagherbandi 2017. Here the resulting error estimates of gravity anomalies on the internal sphere must be regarded as internal errors, which are not affected by an error in the topographic density, and the same conclusion may be drawn for the propagated errors in geoid heights after Stokes integration.

Verification by GNSS-leveling
A most common technique to verify a gravimetric geoid model N grav is to compare it with the so-called geometrically determined geoid height: where h is the geodetic height of the topographic surface and H is the orthometric height. However, as H cannot be determined without topographic gravimetric data, the errors in the two geoid heights (N grav and N geo ) are highly correlated. More precisely, if we write equation (4) in the form and note that h is a pure geometric quantity, it is obvious that any error in the density model causes an error dH in orthometric height, which implies an error dN = −dH in the geoid height. Sjöberg (2018) emphasized that the error in H caused by an erroneous density distribution of the topography yields the same error with opposite sign in the gravimetric geoid height. Hence, this error cannot be verified with GNSS-leveling.

Verification by astro-gravimetric leveling
If two points A and B on the geoid are sufficiently close, their geoid height difference can be approximated by where ε A and ε B are the deflections of the vertical along the tangent from A to B, and s is the distance between the two points. This technique for geoid determination was already reported by Helmert (1880), but traditional astronomic determination of the deflections of the vertical is not accurate enough to cope with today's precise geoid models. However, since a few decades, a new tool, the zenith camera, allows determining deflections of the vertical to an accuracy of the order of 0."1 Flury 2008, Hirt et al. 2010), so that the geoid height difference between two points separated by 1 km could be determined to the order of 1 mm from equation (6). By accumulating several such differences, the error grows to about 4-8 mm for a 38 km polygon. This technique was applied by Westrum et al. (2021) in an attempt to verify the accuracies of gravimetric geoid and quasigeoid heights determined in the Colorado Rocky Mountain test area. Despite the encouraging agreement between gravimetric and astrogeodetic geoid heights, the result must be put at doubt with respect to the uncertainties in topographic density and orthometric heights, which affect all such geoid estimates. This is because the deflections of the vertical are observed at the topographic surface, while equation (6) needs deflections at sea level. This problem was solved by Bodemueller (1957), who formulated the curvature reductions dξ and dη to the observed surface deflections in the north-south and east-west (x,y) directions: where H is the orthometric height, g and ḡ are surface and mean gravity along the vertical and β 1 and β 2 are the northsouth and east-west slopes of the terrain, respectively.
Alternatively, the orthometric reduction between points A and B can be provided from spirit levelling by the formula (Heiskanen and Moritz 1967, p. 168): which should be subtracted from the primary result of equation (6). Here δn represents the levelled height differences and H x , where x = A, B are the orthometric heights at end points. It is obvious that the corrections provided by equations (7a), (7b) and (8) utilize orthometric heights and mean gravity along the vertical from the geoid to topographic surface, i.e., quantities that cannot be precisely estimated without detailed information on the topographic density distribution.

Quasigeoid determination
Verifying quasigeoid models is much easier than geoid models. First, the quasigeoid can be determined without any information on the topographic density distribution. Second, the reduction in zenith camera data from the topography to the geoid is now replaced by the reduction in the telluroid. Hence, in this case, equation (8) is substituted by the normal correction (Heiskanen and Moritz 1975, p. 171): where the mean gravity along the vertical and orthometric heights at the end points A and B are substituted by mean normal gravity (γ) and normal heights (H N ), which can be determined without topographic density data. Also, GNSSlevelling, yielding the quasigeoid height by the equation is an excellent independent tool in verifying a quasigeoid model.

Discussion
Let us return to the equation and consider the following facts: As the problem of determining the geoid (as well as the orthometric height) is an inverse problem (due to the partly unknown density distribution of the topographic density), equation (11) is met not only for the true modeling of the topography but for any such model. However, if the topographic model differs for geoid and orthometric height models, equation (11) fails.
There are national and regional geoid and orthometric height models based on specifically defined topographic models, which may or may not agree. If they disagree, any GNSS-or astro-gravimetric leveling validation will, in principle, fail. As an example, Kingdon et al. (2005) rigorously determined the orthometric heights for Canada, which should (only) match a national geoid model based on the same advanced topographic model.

Concluding remarks
Today moderate or ultra-high degree global gravitational models in combination with regional gravity data are frequently compared for geoid determination and validated by GNSS-leveling and/or astrogeodetic leveling at the few centimeter-level. However, as shown in this article the two techniques using GNSS-or astrogeodetic-leveling cannot detect errors in the geoid models caused by erroneous topographic density distribution data. Already a topographic density error of order 10% may propagate to geoid errors exceeding several centimeters in mountainous regions with elevations reaching 2 km or higher.
This validation problem does not refer to quasigeoid determination, which does not depend on the topographic density distribution. More precisely, the height anomaly (= quasigeoid height) can be determined by the formula T P /γ Q , where point P is located on the topographic surface and normal gravity γ Q is located at normal height. Assuming that the Earth's surface is accurately known from geometric geodetic methods, γ Q replaced by normal gravity at Earth's surface yields the height anomaly good to within 3 mm at any place on Earth, and the accuracy can be further improved by iteration.
The geoid problem is an inverse problem as well as a free boundary value problem (BVP), while determining the quasigeoid is a forward problem and fixed BVP. It is not only the lack of knowledge of the topographic density distribution that complicates the geoid task, but (to some extent) this problem is also due to the missing knowledge of the zero-level of the topography, i.e., the unknown geoid.