Geoid-to-Quasigeoid Separation Computed Using the GRACE/GOCE Global Geopotential Model GOCO02S - A Case Study of Himalayas and Tibet

The geoid-to-quasigeoid correction has been traditionally computed approximately as a function of the planar Bouguer gravity anomaly and the topographic height. Recent numerical studies based on newly developed theoretical models, however, indicate that the computation of this correction using the approximate formula yields large errors especially in mountainous regions with computation points at high elevations. In this study we investigate these approximation errors at the study area which comprises Himalayas and Tibet where this correction reaches global maxima. Since the GPS-leveling and terrestrial gravity datasets in this part of the world are not (freely) available, global gravitational models (GGMs) are used to compute this correction utilizing the expressions for a spherical harmonic analysis of the gravity field. The computation of this correction can be done using the GGM coefficients taken from the Earth Gravitational Model 2008 (EGM08) complete to degree 2160 of spherical harmonics. The recent studies based on a regional accuracy assessment of GGMs have shown that the combined GRACE/GOCE solutions provide a substantial improvement of the Earth’s gravity field at medium wavelengths of spherical harmonics compared to EGM08. We address this aspect in numerical analysis by comparing the gravity field quantities computed using the satellite-only combined GRACE/GOCE model GOCO02S against the EGM08 results. The numerical results reveal that errors in the geoid-to-quasigeoid correction computed using the approximate formula can reach as much as ~1.5 m. We also demonstrate that the expected improvement of the GOCO02S gravity field quantities at medium wavelengths (within the frequency band approximately between 100 and 250) compared to EGM08 is as much as ±60 mGal and ±0.2 m in terms of gravity anomalies and geoid/quasigeoid heights respectively.


InTROdUCTIOn
The orthometric and normal heights are the most commonly adopted concepts for realizations of the vertical reference frames worldwide.These height systems are related to the definition of the geoid and the telluroid (or more commonly used quasigeoid).An appropriate method for the evaluation of the orthometric heights has been discussed for more than a century.The first theoretical attempt is attributed to Helmert (1890).In Helmert's definition of the orthometric height, the Poincaré-Prey gravity gradient is used to evaluate the approximate value of mean grav-ity from gravity observed at the surface point.Later, Niethammer (1932) and Mader (1954) took into account the mean value of the gravimetric terrain correction within the topography.More recently, Hwang and Hsiao (2003) introduced additional corrections due to vertical and lateral variations in the topographic mass density distribution.A rigorous method of computing the orthometric heights using gravity data and digital terrain and density models was introduced by Tenzer (2004) and Tenzer et al. (2005).The comparison of various types of the orthometric heights was done by Santos et al. (2006).Asserting that the topographic mass density and the actual vertical gravity gradient inside the Earth could not be determined precisely, Molodensky (1945Molodensky ( , 1948) ) formulated a theory for normal heights.According to his concept, the mean actual gravity within the topography is replaced by the mean normal gravity between the reference ellipsoid and the telluroid.An approximate formula relating the normal and orthometric heights (i.e., the geoid-to-quasigeoid correction) can be found in Heiskanen and Moritz [1967, Eq. (8-103)].It was demonstrated by several authors that this approximate definition of the geoid-to-quasigeoid correction is not accurate especially in mountainous regions.To improve the computational accuracy, Tenzer et al. (2006) formulated the explicit relation between normal and orthometric heights.An alternative expression was given by Sjöberg (2006).His formulation includes a topographic roughness term (topographic bias) and a term related to the lateral mass density variation within the topography for the practical application of the height conversion.The explicit formula for the geoid-to-quasigeoid correction in Tenzer et al. (2005) requires the discretised numerical integration which is very time consuming.For a fast computation, the spectral expression developed by Sjöberg (2006) can then be applied instead and achieve sufficient numerical accuracy.
The practical computation of the geoid-to-quasigeoid correction according to Sjöberg (2006) requires three types of global models: the global geopotential model (GGM), global terrain, and topo-density models.Whereas a compilation of the global topo-density models is still restricted by the absence of reliable global density data, the global terrain models are available to a very high accuracy and resolution.The global topographic model DTM2006.0 (Pavlis et al. 2007) is provided with a spectral resolution up to degree 2160 of spherical harmonics.The same spectral resolution is used to describe the gravitational field by the coefficients of the Earth Gravitational Model 2008 (EGM08; Pavlis et al. 2012).The DTM2006.0 and EGM08 were made available by the US National Geospatial-Intelligence Agency EGM development team.The EGM08 and DTM2006.0 coefficients complete to degree/order 2160 could be used to compute the geoid-to-quasigeoid correction (Sjöberg and Bagherbandi 2012).
In association with the precise modeling of gravity field two satellite missions provide information which considerably increased the accuracy of existing GGMs and significantly improved its applications in various scientific fields: the Gravity Recovery and Climate Experiment (GRACE) launched in 2002 and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) launched in 2009.These missions are somewhat complementary.The GOCE gravity gradiometry provide the gravity field to a spatial resolution of ~100 km or better but is relatively inaccurate at long wavelengths (above 700 -1000 km).The GRACE Kband ranging inter-satellite observations provide the precise information on a low-frequency gravity spectrum.The combined GOCE/GRACE solutions thus provide a more accu-rate estimation of the Earth's static gravity field (as well as its temporal variations).Since the EGM08 coefficients were generated without using GOCE data, the accuracy of the computed gravity field quantities at medium wavelengths should improve when using the latest available GGMs.In this study we use the satellite-only combined GOCE/ GRACE model GOCO02S (Goiginger et al. 2011), which is provided with a spectral resolution complete to degree 250 of spherical harmonics, to compute the gravity field quantities (the free-air gravity anomaly and geoid-to-quasigeoid correction values).The GOCO02S gravity field quantities are then compared with the EGM08 results in order to assess the expected improvement of gravity field at medium wavelengths.The study area is situated in Himalayas and Tibet, where the investigated gravity field quantities reach extreme values.

THEORy
The computation of gravity field quantities from the GGM coefficients utilizes a spectral representation of gravity field by means of spherical harmonics.The disturbing potential T (i.e., the difference between the Earth's gravity potential W and the normal gravity potential U) at an arbitrary point (r, Ω) is computed using the following expression (e.g., Heiskanen and Moritz 1967) where R = 6371 × 10 3 m is the Earth's mean radius (which approximates the geocentric radii of the geoid surface), Y n, m are the (fully normalized) spherical harmonic functions of degree n and order m, T n, m are the (fully normalized) coefficients of the disturbing potential, and nmax is the upper summation index of spherical harmonics.The 3-D position is defined in the system of spherical coordinates (r, Ω), where r is the spherical radius and Ω = (φ, λ) denotes the spherical direction with the spherical latitude φ and longitude λ.
The geoid height N is defined by the well-known Bruns' (1878) formula in the following form where 0 c is the normal gravity calculated at the reference ellipsoid using, for instance, Somigliana's (1929) formula.The disturbing potential T in Eq. ( 2) is computed at the geoid surface.The geocentric radius of the geoid surface is denoted as r g .
By analogy with Eq. ( 2), the quasigeoid height ζ (i.e., the height anomaly) is defined as (Molodensky et al. 1960) where Q c is the normal gravity at the telluroid.The computation of the normal gravity value Q c can be done according to the expression given in Heiskanen and Moritz [1967, Eq. (2-123)].The disturbing potential T is in this case calculated at the surface point.The geocentric radius of the surface point is denoted as r P .
The geoid-to-quasigeoid correction is defined approximately as a function of the simple planar Bouguer gravity anomaly Δg SB and the topographic height H of computation point.With reference to Heiskanen and Moritz (1967, Chapters 8-12 and 8-13), we write where c is the integral mean of the normal gravity along the normal between the reference ellipsoid and the telluroid.
The simple planar Bouguer gravity anomaly Δg SB in Eq. ( 4) is given by (ibid.) where G = 6.674 × 10 -11 m 3 kg -1 s -2 is Newton's gravitational constant, o t t is the mean topographic mass density, and Δg is the (free-air) gravity anomaly.We note that either the normal or orthometric heights can be used to compute both the geoid-to-quasigeoid correction [in Eq. ( 3)] and the Bouguer gravity anomaly [in Eq. ( 4)].
A more accurate expression for computing the geoidto-quasigeoid correction in the spectral domain was given by Sjöberg (2006); see also Sjöberg and Bagherbandi (2012) The last constituent on the right-hand side of Eq. ( 6), i.e., Vbias t c , defines the topographic bias.This bias represents the error in the analytical downward continuation of the external gravitational potential inside the topographic masses.For the adopted constant topographic mass density distribution, the term Vbias t is computed according to the following spectral expression (Sjöberg 2007) The computation of the topographic bias in the spectral domain is numerically stable when using a low-degree series.Ågren (2004) and Sjöberg (2007), however, demonstrated that the computation became unstable at a very-high degree spherical harmonic terms of a power series of topographic heights.The terms : , , ,...
/ in Eq. ( 7) define the spherical height functions : , , ,... i.e., where P n is the Legendre polynomial of degree n with the argument t equal to cosine of the spherical distance ψ between the spherical directions X and Xl; i.e., t = cosψ.

RESUlTS
The numerical analysis was conducted at the study area bounded by the parallels of 22 and 60 arc-deg northern latitudes and the meridians of 60 and 120 arc-deg eastern longitudes comprising the Himalayas and Tibet Plateau.All computations were realized on a 1 × 1 arc-deg geographical grid.
The GOCO02S coefficients were used to compute the gravity field quantities with a spectral resolution complete to degree 250 of spherical harmonics.The same spectral resolution was used to compute the gravity field quantities using the EGM08 coefficients.The results from both models were compared.The topographic bias was computed in spectral domain using the coefficients of the global topographic model DTM2006.0 complete to degree/order 250.The normal gravity field quantities were computed according to the GRS-80 parameters (Moritz 1980).The average density of the upper continental crust 2670 kg m -3 (see Hinze 2003) was adopted as the topographic mass density.We note here that authors of EGM08 and DTM2006.0 recommended using the full series expansion (up to the maximum available degree and order).However, for the comparison with the GOCO02S results we used the same maximum degree of expansion for all three models.
The GOCO02S and EGM08 coefficients were used to generate (free-air) gravity anomalies.The regional map of the GOCO02S gravity anomalies is shown in Fig. 1.The GOCO02S gravity anomalies vary between -182 and 217 mGal with the mean of -10 mGal and the standard deviation is 39 mGal.The differences between the GOCO02S and EGM08 gravity anomalies within the study area are shown in Fig. 2.These differences vary from -51 to 60 mGal; the mean and RMS of differences are 0 and 7 mGal respectively.
The DTM2006.0 coefficients were used to generate the topographic bias according to Eq. ( 7).The result is shown in The geoid-to-quasigeoid correction was computed according to Eq. ( 6) using the GOCO02S and DTM2006.0 coefficients.The result is shown in Fig. 4. It varies between -0.15 and 3.62 m with the mean of 0.32 m and the standard deviation is 0.80 m.The same computation was done using the EGM08 and DTM2006.0 coefficients.The differences between the values of the geoid-to-quasigeoid correction computed using the GOCO02S and EGM08 coefficients (complete to spherical harmonic degree 250) are shown in

dISCUSSIOn
The practical computation of the geoid-to-quasigeoid correction according to Eq. ( 6) was done in two numerical steps.First, the disturbing potential values at the surface and geoid points were generated from the GOCO02S (as well as EGM08) coefficients.The topographic bias term was then calculated using the DTM2006.0 coefficients.As seen from the results in Figs. 3 and 4, the topographic bias represents a significant contribution to the geoid-to-quasigeoid  correction.According to Sjöberg and Bagherbandi (2012), this contribution is about 90%, while only about 10% is attributed to the disturbing potential term.The topographic bias in Eq. ( 7) is defined as a function of variable topography while adopting a constant topographic density.This term is thus strongly correlated with the height of computation point as well with the surrounding terrain configuration (especially in the close proximity of computation point).The consideration of surrounding terrain in the topographic bias is the principal difference compared to the approximate definition in Eq. ( 4), where only the topographic height of computation point is taken into account in the functional model.For a more accurate computation, the constant topographic density in Eq. ( 7) should be replaced by a more realistic mathematical model of the actual topographic density distribution.Since current knowledge about the actual density distribution is still restricted by the lack of geological (rock density) data especially in the chosen study area, consideration of a more realistic density model is not applicable.Another substantial difference between the approximate and accurate computational models in Eqs. ( 4) and ( 6) is associated to a different definition of gravity values.In the approximate formula [Eq.( 4)], the free-air gravity value at computation surface point is assumed.Note that the Bouguer gravity reduction is here formally considered as the topographic term [with analogy to the topographic bias in Eq. ( 6)].Conversely, the difference between the disturbing potential values on the surface and geoid points in the first two constituents on the right-hand side of Eq. ( 6) equals the integral mean of the gravity disturbance along the plumbline within the topography (multiplied by the topographical height).Hence, the approximate definition in Eq. ( 4) uses the free-air gravity anomaly at the surface point, while the accurate model utilizes implicitly the mean value of the gravity disturbance along the plumbline within the topography.For more details about these theoretical aspects of a rigorous definition of the geoid-to-quasigeoid correction we refer readers to Tenzer et al. (2006).
The inaccuracies in computed values of the geoid-toquasigeoid correction are mainly due to three sources of errors.The largest contribution to the total error budget is more likely due to uncertainties within the actual topographic mass density distribution.The approximation error due to using the constant topo-density value introduces the same relative error in the computed topographic bias, provided that the density errors propagate proportionally to the errors in the topographic bias.For the maximum value of the topographic bias of 3.7 m (see Fig. 3), the relative density error of 10% corresponds to the maximum error in computed topographic bias of about 37 cm.It is obvious that the largest errors due to density uncertainties are mainly found in mountainous regions with variable geological structure.The inaccuracies in computed correction values caused by the errors of terrain model can formally be attributed to the height error of computation point and the height errors within surrounding terrain configuration.If the height errors of surrounding terrain have a prevailing random character, they cancel out by averaging.The most significant then becomes the errors in heights of computation points.To quantify this error, we can use the approximate formula in Eq. (4).When substituting from Eq. (5) to Eq. ( 4), it is clear that the largest inaccuracy in the computed value of the geoidto-quasigeoid correction is due to the height uncertainty in the computed Bouguer reduction term.An error analysis according to the Gauss error propagation law reveals that the errors in computed correction are N fgrelated to the height error For the maximum topographic heights (< 9 km), the height error of 10 m causes the error of about 2 cm in the geoid-to-quasigeoid correction.For the elevation of 1 km, the same height error of 10 m causes an inaccuracy in the computed correction value of only about 2 mm.In our numerical realization, the errors in computed correction are due to the DTM2006.0 omission and commission errors.The third main source of error is due to the uncertainties in input gravity data.In our numerical scheme these errors are due to the omission and commission errors of used GGMs.Here we investigated the errors as a result of the EGM08 commission errors at medium wavelengths.This has been done by comparing the GOCO02S and EGM08 results.This comparison revealed large differences.As seen in Fig. 2, the differences in the computed (free-air) gravity anomalies are within ±60 mGal.The corresponding differences in the geoidal heights shown in Fig. 4 are within ±0.2 m.We note that the differences in the computed values of the geoid-to-quasigeoid correction very closely agree with the geoid/quasigeoid height differences between GOCO02S and EGM08.The differences in computed values of the gravity anomalies and geoid/quasigeoid are obviously highly spatially correlated but the gravity differences have a prevailing higher-frequency pattern.As also seen in Figs. 2 and 4, the spatial structure of these Fig. 6.Approximation errors in the computed values of the geoid-to-quasigeoid separation using Eq. ( 4).Values are in meters.
differences indicates that they are attributed to a mediumwavelength GGM spectrum.This is confirmed based on the analysis of the degree and cumulative geoidal height differences between GOCO02S and EGM08 (see Fig. 6).The cumulative geoid height differences N f up to degree N were computed from The significant degree geoid height differences are seen approximately above spherical harmonic degree 100.It was demonstrated in recent studies that the GOCE gravity gradiometry data improved the accuracy of the GGM coefficients at this part of the gravity spectrum (see e.g., Pail et al. 2010;Goiginger et al. 2011).The low degrees of GOCO02S are primarily determined by GRACE, whereas the GOCE gravity gradiometry observables start to significantly contribute at degree about 100.Beyond degree 150, the combined model is dominated mainly by GOCE.
In geodetic literature (e.g., Heiskanen and Moritz 1967), the geoid-to-quasigeoid correction is conventionally defined approximately as the difference between Helmert's (1890) orthometric and Molodensky's (1945) normal heights.The resulting formula for computing this correction is then defined as a function of the planar simple Bouguer gravity anomaly and the topographic height of computation point [cf.Eq. ( 4)].However, this approximation yields large errors in the computed correction values, particularly in mountainous regions.Santos et al. (2006) demonstrated that the inaccuracies revealed via numerical examples at the study area situated in the Rocky Mountains using Helmert's orthometric heights can reach centimeters and locally even decimeters.The computation thus typically requires the application of a more accurate formula [which is given in Eq. ( 6)].The inaccuracy caused by using the approximate expression in Eq. ( 4) for computing the geoid-to-quasigeoid correction is shown in Fig. 7.The approximation errors [i.e., differences between the solutions obtained using accurate and approximate formulas in Eqs. ( 6) and ( 4) respectively] are between 1.56 and -0.67 m; the mean and RMS of approximation errors are -0.07 and 0.14 m respectively.The largest positive differences were found at computation points of which elevations to exceed about 5 km.The values of the geoidto-quasigeoid correction computed according to Eq. ( 6) are thus significantly larger than the corresponding values computed using the approximate formula in Eq. ( 4).Bretterbauer (1986) assumed that the relation between the geoid-toquasigeoid correction and elevation is linear.The non-linear increase of this correction with elevation was demonstrated, for instance, by Flury and Rummel (2009).The non-linear character between these two quantities is seen also in Fig. 8.The correlation between these two quantities is 0.95.

SUMMARy And COnClUSIOnS
We have applied accurate expression for computing the geoid-to-quasigeoid correction at the study area where this correction reaches global maxima.The expression which utilizes the spherical harmonic representations of the Earth's gravity field and topography was used to calculate this correction.The topographic bias was computed using the DTM2006.0 coefficients.The correction term related to the gravity field (disturbing potential values at surface and geoid points) was evaluated using GOCO02S and compared with the results obtained from EGM08.All values were computed with a spectral resolution complete to degree 250 of spherical harmonics (which is the maximum resolution of GOCO02S).The constant value of topographic density 2670 kg m -3 was adopted.The computed correction values were compared with the corresponding results obtained after using the approximate formula.The results were then compared and analyzed.
Our results revealed that the geoid-quasigeoid correction (computed with a spectral resolution up to degree/order 250) vary between -0.15 and 3.62 m.This correction is highly correlated with the topography with the maxima at computation points with the largest elevations.These values differ significantly from the values computed using the approximate formula.The geoid-to-quasigeoid correction computed approximately reaches the maxima up to about 2 m.This value is much smaller than the maximum of 3.62 m obtained when using the accurate formula (computed for the maximum degree 250 of spherical harmonics).Sjöberg and Bagherbandi (2012) estimated that the maximum value of this geoid-to-quasigeoid correction is 5.47 m when using the EGM08 and DTM2006.0 coefficients with a spectral resolution complete to degree 2160 of spherical harmonics.The errors due to using the approximate formula [in Eq. ( 4)] thus can even exceed the value of the correction itself.
The GOCE gravity gradiometry data improved the accuracy of GGMs especially at medium wavelengths (within a frequency bound from about 100 to 250 of spherical har-Fig.7. Scatter plot of the geoid-to-quasigeoid separation with the heights of observation points.monics).Our results have shown that this improvement is significant; the differences found between the computed GO-CO02S and EGM08 gravity field quantities reach ±60 mGal (for the gravity anomalies) and ±0.2 m (for geoid-to-quasigeoid correction and consequently also the geoid/quasigeoid heights).These results indicate that the EGM08 commission errors at medium wavelengths cause almost the same errors in computed correction values as the topo-density uncertainties.

Fig. 3 .
Fig. 3.The topographic bias within the study area is everywhere positive and reaches the maxima of 3.7 m.The mean value is 0.4 m and the standard deviation is 0.8 m.The geoid-to-quasigeoid correction was computed according to Eq. (6) using the GOCO02S and DTM2006.0 coefficients.The result is shown in Fig.4.It varies between

Fig. 1 .
Fig. 1.GOCO02S gravity anomalies generated on a 1 × 1 arc-deg geographical grid with a spectral resolution complete to spherical harmonic degree 250.Values are in mGal.

Fig. 2 .
Fig. 2. Differences between the GOCO02S and EGM08 gravity anomalies generated on a 1 × 1 arc-deg geographical grid with a spectral resolution complete to spherical harmonic degree 250.Values are in mGal.

Fig. 5 .
Fig. 5.These differences are between -0.20 and 0.13 m; the mean and RMS of differences are 0.00 and 0.01 m respectively.

Fig. 5 .
Fig. 5. Differences between the GOCO02S and EGM08 values of the geoid-to-quasigeoid separation computed on a 1 × 1 arc-deg geographical grid.Units are in m.The topographic bias for both GGM solutions was computed using DTM2006.0.