1. Introduction
The International Association of Geodesy (IAG) establishes a coordinated approach for monitoring global changes inside the Global Geodetic Observing System (GGOS). The structure of GGOS/IAG is able to provide control information and products for prediction and control of several geo-phenomena like Earth´s orientation in the space, mass transport in the Earth´s System, geohazards monitoring, sea level changes, local variability and forecasting (^{Plag and Pearlman, 2009}). A unified reference for heights is fundamental for activities related to the modelling of regional and global changes as established in the GGOS theme 1 - Unified Height System (^{Ihde, Sánchez and Sideris, 2010}); (^{Kutterer and Neilan, 2015}).
An essential aspect which must be considered, is that the discrepancy of each Vertical Datum (VD) to the referred global equipotential surface propagates to each functional which depends on the heights (e.g. gravity anomalies). This is the called VD indirect effect (^{Gatti, Reguzzoni and Venuti, 2012}). To avoid this problem, it is necessary to refer all national VDs to the same global equipotential surface.
Considering the referred aspects, the IAG fixed the foundations of the IHRS by the Resolution 1/2015 (IAG, 2015), where the global equipotential reference surface is fixed by the geopotential value W _{ 0 } = 62 636 853.4 m^{2}s^{-2}, and the primary vertical coordinates of points P _{ i } on the Earth´s surface are given by the geopotential numbers expressed as C _{ Pi } =W _{ 0 } - W _{ Pi } .
Several strategies for connecting VDs in the geopotential space started to be developed in the last three decades, as a consequence of satellite based geodetic approaches and modern gravimetry. More details can be found in: (^{Heck and Rummel, 1990}); (^{Rummel, 2000}); (^{Lehmann, 2000}); (^{Amjadiparvar, Rangelova and Sideris, 2016}). Nowadays, among the main activities of the international geodetic community are the search of conventions and standards for realizing the IHRS, as well as the determination of the offsets of each national VD related to the new conventional global equipotential surface.
The SIRGAS WG III (Geocentric Reference System for the Americas project - Working Group III - Vertical Datum), has among its tasks the unification of the National Vertical Reference Networks (NVRNs) of member countries in terms of physical heights and a common VD. Since 1997 the SIRGAS WG III establishes several activities like workshops and schools on Vertical Reference Systems in South America, Central America and Caribe. Always, these activities are directed for contributing to the development and diffusion of IAG resolutions, standards related to the connection of NVRNs, and their link to a common global reference, now the IHRS (^{de Freitas, 2015}).
The contents in this paper are inserted in the purposes of SIRGAS WG III, by introducing alternatives for determining the offset of the Ecuadorian Vertical Datum (EVD), related to the IHRS in the geopotential space. The free Geodetic Boundary Value Problem (GBVP) solution based on heterogeneous data, allow determining the disturbing potential in the EVD, and in its 4° x 4° surrounding region. The size of the EVD adjacent region was chosen with the purpose of considering a representative area containing gravimetric information, however, in future studies, variations in the size of the EVD adjacent region will be tested. These variations will be tested mainly with the aim of avoiding the effect of omission errors of global models in mountainous regions.
The heterogeneity of data sources considering distribution, reference systems, tide systems, reference epochs, precisions and different spectral resolutions, were the main undergoing aspects. In addition to geodetic point-data, we also considered global databases available in Global Geopotential Models (GGMs) and Digital Elevation Models (DEMs). The Least Squares Collocation (LSC) method and Remove-Restore (RR) spectral decomposition technique seem to be the most adequate ways for facing of the referred data basis with different spectral characteristics (^{Barzaghi, 2016}). We applied also the Residual Terrain Model (RTM) technique for modelling the short wavelengths of the gravity field (^{Hirt, Featherstone and Marti, 2010}) (^{Tziavos and Sideris, 2013}).
In 2010, IGM-EC (Ecuadorian Military Geographic Institute) carried out the preliminary adjustment of the Ecuadorian Fundamental Vertical Network (EFVN), adopting as its origin the Mean Sea Level (MSL) determination for the observation period 1988-2009, and materializing the vertical reference in the BM03 bench mark, which was located 6.2707 m above the MSL. Thus, the vertical offset computation on the Equatorial Vertical Datum by means of the solution of the fixed Boundary Value Problem will be performed considering the BM03 as the computation point.
2. Theoretical Foundations
The geopotential in a point P is given in terms of the normal potential U _{ P } and by the disturbing potential T _{ P } by:
The normal potential is obtained from a Geodetic Reference System (GRS), based on the reference level ellipsoid (in our case associated to the GRS80) (^{Hofmann-Wellenhof and Moritz, 2006}), and considering the gradient of the normal potential in the form:
where h _{ P } is the ellipsoidal height. Then, the local geopotential offset (e.g. in the EVD related to the IHRS) can be expressed by:
where the global reference value of IHRS is given by W _{ 0 } and the geopotential value associated with the local VD is W _{ 0i } . The gradient of the normal potential ?U _{ 0 } /?h is approximately the normal gravity value (?γ). In this approach, the main problem is to determine the T _{ 0i } (local disturbing potential) value. The disturbing potential is determined based on some strategy for solving the GBVP.
The adopted strategy for computing the disturbing potential T in this work, concerns with free-air surface gravity anomalies or Molodensky’s gravity anomalies, which do not imply (like in the classical solution of the GBVP by the Stoke's integral) arbitrary reductions of known gravity values at Earth’s surface to the geoid (for details see ^{Hofmann-Wellenhof and Moritz, 2006}, pp. 290-293). The gravity can be obtained from several sources like direct observations on the continent and on the oceans, as well indirectly from aerial gravimetry, satellite gravity, altimetry missions, and by modelling the direct and indirect effects of topographic mass. In the present work, continental and oceanic areas are involved in the GBVP solution.
Considering the usual reference surface related to heights with physical meaning presented in the Figure 1, we adopted the determination of the disturbing potential based on the Molodensky’s gravity anomalies (∆g _{ M } ), including the atmospheric (δg _{ A } ) and the Honkasalo (∆g _{ H } ) corrections (see eq. 4). This kind of anomaly is given by:
The functional γ_{ ∑ } is the normal gravity value at the telluroid (see Figure 1), obtained from γ. The normal gravity value in the level ellipsoid is given e.g. by the Clairaut’s formula and depending of the latitude ϕ Erro! Fonte de referência não encontrada. (^{Hofmann-Wellenhof and Moritz, 2006}):
In this equation γ_{ e } is the normal gravity at the equator and the other involved parameters are:
In the previous equations ω is the angular velocity of the Earth, ? = (a - b)/a is the geometrical flattening of the level ellipsoid, with a and b being respectively the major and minor semi-axis. Other equivalent possibility for computing the normal gravity is to consider the Somigliana’s formula:
where γ_{ p } is the normal gravity at the pole. All the parameters must be related to a reference ellipsoid, for this work were considered those related to the GRS1980.
The normal gravity above the level ellipsoid is obtained by considering the normal gravity gradient as proposed by Bruns (details in ^{Hofmann-Wellenhof and Moritz, 2006}, p. 81-82). Considering the ellipsoidal height h, in a rigorous form, the normal gravity above the level ellipsoid (γ_{ h } ) is given by Erro! Fonte de referência não encontrada.; and considering a simplification given by the mean value for the normal gravity gradient, the normal gravity above the level ellipsoid is given in mGal by Erro! Fonte de referência não encontrada.:
According to ^{Heikkinen (1979}), the ^{Honkasalo term (Honkasalo, 1964}), which removes the average part of the tidal force, must be include when gravity data are referred to the International Gravity Standardization Net (IGSN71). However, this term has been deemed inappropriate because of resulting errors in the GBVP solution. Therefore, following the recommendation of the International Association of Geodesy (IAG), the Honkasalo term has to be removed from IGSN71 gravity values (^{Uotila, 1980}). In the equation Erro! Fonte de referência não encontrada., the Honkasalo term is removed by adding the latitudinal-dependent Honkasalo correction Δg _{ H } (12) in mGal (^{Hinze et al., 2005}):
2.1 Computation of the disturbing potential
The Fredholm integral over the telluroid for computing the disturbing potential at the Earth’s surface based on the Molodensky’s theory (^{Hofmann-Wellenhof and Moritz 2006}, p. 311-313) is given by:
where n means the normal direction and is the euclidean distance from the integration element to the computation point. It presumes an interactive solution. This approach can be substituted by decomposition in the series T _{ P } = T _{ 0 } +T _{ 1 } +T _{ 2 } +... in the general form:
Which can be decomposed as:
where R is the mean Earth’s radius and the ensemble of terms g _{ 0 } , g _{ 1 } , g _{ 2 } … is known as Molodensky’s series (^{Hofmann-Wellenhof and Moritz (2006}, p. 311-313). The terms are associated with local corrections and criteria for linearization. Usually, the computational methods for the Molodensky’s solution involve only the two first corrections given by (^{Gemael, 2012}, pp. 270):
where, h and h _{ P } are the heights of the integrating and computing points respectively and is the gravity mean value.
In Erro! Fonte de referência não encontrada. ψ is the angle of the geocentric distance between the computation point and the point where the gravity anomaly is known. The Stokes’s function is given by:
The spherical distance can be computed based on the coordinates of the referred points (ϕ, λ) and (ϕ’, λ’) respectively as (^{Gemael, 2012}, p.146):
The Molodensky’s solution is applied for a geocentric GRS. For an arbitrary GRS, it is rewritten as:
where, ΔGM is the difference between the geocentric gravitational constant in the global and the local systems. Sometimes the term T _{ zero } = ∆GM/R is called zero-degree term.
In practices, the computation of the Molodensky’s integral is adequate to be solved in the frequency domain based on the Fast Fourier Transform (FFT) technique or by LSC approach (^{Sansò and Sideris, 2013}). These approaches allow the spectral decomposition of the disturbing potential in the so-called Remove-Restore (RR) process (^{Forsberg and Tscherning, 1981}). The solution by spectral decomposition allows using related gravimetric information coming from different sources. It is the case where it is possible to integrate computed gravity anomalies, predict values from modern GGMs, and values from terrain effects like in the RTM technique based on DEMs (^{Forsberg, 1984}); (^{Hirt, Featherstone and Marti, 2010}). In this sense it is possible to consider the RR process by the spectral decomposition as:
The gravity anomaly corresponding to the RTM effect can be obtained from the residual high frequency topography, given by a high degree DEM that is subtracted by a mean long wavelength surface corresponding to the maximum harmonic degree of the used GGM. The corresponding RTM residual gravity anomaly in a planar approximation (^{Forsberg, 2008}) is given by:
where, G is the universal gravitational constant, ρ is the standard density of the crust, and h is the topographic height given by a DEM. The integral equation Erro! Fonte de referência não encontrada. is solved by using only the small residual Molodensky’s gravity anomalies Erro! Fonte de referência não encontrada.. Then it is possible to compute T _{ RES } , and the completely disturbing potential is obtained in a composition process by:
Analogously, the height anomaly is compute as:
The equations Erro! Fonte de referência não encontrada. and Erro! Fonte de referência não encontrada.are related by the Bruns formula Erro! Fonte de referência não encontrada., and from Erro! Fonte de referência não encontrada. is possible to compute the geopotential (W _{ P } ) by using the equation (1).
3. EVD Characteristics and Used Data Sets for GBVP solution
The IGM-EC in cooperation with the US Inter American Geodetic Survey (IAGS) established the EVD in 1948. The EVD is associated with the La Libertad’s tide gauge ( Figure 2). The reference Mean Sea Surface (MSS) level was materialized by the INOCAR (Ecuadorian Oceanographic Institute of the Naval Forces), in cooperation with the IGM-EC (Ecuadorian Military Geographic Institute) by the fundamental bench mark (BM) called BM03 (λ = -80°54'19.4667"; φ = -2°13'10.1178") (IGM-EC, 2013). The surrounding datum area involves the ocean-continent transition and has an irregular bathymetry/topography in the range from -4698 m to 1552 m ( Figure 2).
In the proposed method, to model the disturbing potential in the EVD, we consider in situ data (terrestrial and ocean) in a 4° x 4° region centered in the EVD, and global information coming from GGMs and DEMs. We also consider Sea Surface Topography (SST) and ocean gravity anomalies information coming from satellite altimetry.
Because the records of the gravity dataset have heterogeneous characteristics, it is necessary to perform an outliers filtering, the adopted statistical criteria for detecting and filtering outliers was the three-sigma rule. We considered as outliers all the records whose differences ( ) between computed and modeled gravity anomalies are larger than three times its standard deviation.
Terrestrial gravity: Terrestrial gravity data coming from gravimetric surveys carried out within the period 1960 - 2015. This gravity information was surveyed by the IGM-EC and by oil companies for geodetic and oil prospection purposes respectively. A set of 1986 terrestrial gravity data records were used in this work (Figure 2, brown points). For all gravity records the latitude, longitude and levelled height are known. The terrestrial gravity records have heterogeneous characteristics, principally in terms of observations epochs, observation techniques, and data processing; moreover, the gravity observation precisions are unavailable. For these reasons is necessary to use a statistical method to remove outliers. The details of this process are described in the methodology section.
Ocean gravity: The ocean gravity data set, was provided by the Bureau Gravimetric International (BGI) data base and contains 8549 stations at the 4°x4° surrounding datum area ( Figure 2, blue points). This gravity records belong to 21 gravity surveys carried out within the period 1957 - 1987 for multiple purposes. For each gravity record it is also known the latitude, longitude, free-air gravity anomaly and Bouguer gravity anomaly. As in the case of the terrestrial gravity, the precisions for the ocean gravity records are unavailable and an outlier filtering, based on a statistical analysis was performed (details in the methodology section).
DTU15 marine gravity anomalies: The BGI ocean gravity data is complemented by marine gravity anomalies from the Global Marine Gravity Model (GMGM) DTU15 derived from satellite altimetry (^{Andersen and Knudsen, 2016}). The DTU15 is based on five years of altimetry information from Criosat-2 and Jason-2 missions. The DTU15, originally with a spatial resolution of 1 arc minute was resampled to 5 arc minutes (Figure 2, green points) aiming to reduce the computational cost, mainly for the RTM process (residual terrain effect on gravity anomalies computations). This resampling was performed without apparent prejudice to the quality of the final solution in view of the small variability of such anomalies in the ocean areas.
Global Models: Besides the in-situ information, GGMs, DEMs and the GMGM were used in the RR technique for modelling the gravity field. The GGMs and GMGM allow modelling the long and intermediate wavelengths, while DEMs allow modelling the short ones. In this sense, the GGMs GO-CONS-GCF-2-DIR-R5 (DIR-R5) (^{Bruinsma et al., 2013}) with harmonic development until the degree and order 300, and the EIGEN-6C4 (^{Förste et al., 2014}) developed until the degree and order 2190, each one associated with the GMGM DTU15, were used aiming to produce two independent solutions. The DIR-R5 model was generated through the combination of information from the GOCE, GRACE and LAGEOS missions. Comparing the DIR-R5 with previous versions of the GOCE direct approach, one can appreciate that its spectral behavior shows improvements while the formal errors decrease. On the other hand, the EIGEN-6C4 model has been generated through the combination of information from the LAGEOS, GRACE, GOCE and EGM2008 missions, gravity information from the DTU global model and terrestrial data. The EIGEN-6C4 presents improvements in terms of spectral behavior and formal errors when compared to the EGM2008. Because the DIR-R5 is a satellite-only model, it does not have the influence of the local vertical datum, and therefore it can be used with its maximum degree of expansion. In the case of the EIGEN-6C4, it’s necessary to truncate its maximum degree expansion, in order to avoid local vertical datum effects.
We considered two DEMs in the computations: The SRTM15_PLUS (^{Olson et al., 2014}) involving bathymetry/topography until the wavelength of 15 arc second (approx. 450m) and the SRTM1 with resolution of 1 arc second (approx. 30m) for the continental area (^{Werner, 2001}). The two DEMs are merged to cover the whole study area (oceanic and continental part). The omitted high frequencies in the GGMs were modeled by using the RTM technique based on the effects of short wavelengths residual bathymetry/topography in the gravity field.
4. Methods
The EVD offset in the geopotential space was estimated by the free GBVP solution. For this purpose, surface gravity anomalies ( ) (terrestrial and oceanic) were computed according to Erro! Fonte de referência não encontrada.. Due to the non-homogeneous spatial distribution of oceanic gravimetry records, marine gravity anomalies from DTU15 model were used to complement the BGI data base. According to a statistical criterion, those gravity records considered as outliers are eliminated from the data set. The RTM technique was used to modelling the high frequencies on the gravity anomalies (
The residual height anomaly (ζ_{ RES } ) was estimated on the EVD location as a function of the residual gravity anomalies ( ), for later compute the height anomaly (ζ) by the restoration of ζ_{ GGM } and ζ_{ RTM } according the equation Erro! Fonte de referência não encontrada.. The disturbing potential (T _{ P } ) is compute from ζ_{ P } through the Bruns’s formula Erro! Fonte de referência não encontrada., the geopotential value (W _{ 0 } ^{ i } ) is computed as function of the disturbing potential from (1), and the EVD offset is computed according to Erro! Fonte de referência não encontrada.. Details about these computations are described in the following subsections.
4.1 Gravity anomalies computation and outlier elimination
Gravity anomalies are computed from direct gravity observations by Erro! Fonte de referência não encontrada., considering the atmospheric and the Honkasalo corrections. These computed anomalies and those coming from the DTU15 are used to compute residual gravity anomalies (in each data point in Figure 2), and later estimate the geopotential value at the EVD through the free GVBP solution. Due to the data heterogeneity and the lack of related metadata, in a first step it was necessary to establish a preliminary check of data consistency. In this sense, we compared the computed continental and ocean gravity anomalies with gravity anomalies coming from global models (GGMs and GMGM). Thus, to filter the BGI data set (ocean gravity anomalies), we used as reference the gravity anomalies from the GMGM DTU15. From this analysis, 66.91% of the marine gravity anomalies were eliminated from the BGI gravity data set. Thus, the standard deviation for the residual gravity anomalies changes from 7.37 mGal to 4.40 mGal, and the correlation coefficient from 0.9888 to 0.9962. The marine gravity outliers belonged to some specific ship tracks, hence we presumed that the high percentage of records considered as outliers is due to errors in data handling and processing. On the other hand, to filter all the gravity records (in situ ocean and continental gravity anomalies and DTU15 gravity anomalies), we used as reference the gravity anomalies from the GGM EGM2008 (^{Pavlis et al., 2008}), also in its maximum degree (n= 2190) and order (m=2159). This analysis concludes that 4.53% of the gravity records are outliers and should be deleted from the gravity data set. A correlation coefficient of 0.9640 and a standard deviation of 16.74 mGal were computed for the residual gravity anomalies before the outlier filtering. In the case of the filtered data set a correlation coefficient of 0.9820 and a standard deviation of 11.71 mGal were computed for the residual gravity anomalies. The gravity anomalies from the model grids in its maximum resolution were interpolated at the gravity records locations by kriging.
4.2 Residual gravity anomalies
Residual gravity anomalies were computed for all the available gravity records according to Erro! Fonte de referência não encontrada.. Geopotential modelling from gravity anomalies is best performed using residual quantities according to the remove-restore method (^{Forsberg and Tscherning, 1997}). The residual gravity anomalies are computed considering the gravity anomalies from the GGMs DIR-R5 (^{Bruinsma et al., 2013}) in its maximum degree and order (n=m=300) and from the EIGEN6C4 with maximum degree (nmax = 1000).
Using the RTM technique, as described by Tziavos and Sideris (^{2013}), masses above a smoothed reference surface (lower resolution DEM) are removed and the voids are filled to compute the short wavelengths of the gravity field related to the effects of the residual topographic masses. For these purposes, two DEMs were merged to estimate the effects of the residual topography: the SRTM15_PLUS with spatial resolution 15" arc (continental and oceanic coverage - bathymetry); and the SRTM1 with spatial resolution 1" arc (only continental coverage).
The spatial resolution of the smoothed reference surface (required for the RTM computations) must be equivalent to that of the GGM used to represent the low frequencies of the gravity field. In this work, this lower resolution DEM was generated by applying a low-pass filtering to the merged DEM. For these purposes, the low pass filtering is performed by a moving average with a moving window which size was chosen experimentally by a statistical analysis of the residual gravity anomalies computed by Erro! Fonte de referência não encontrada. and obtained for different versions of the smoothed DEM (different window size). Thus, the chosen window size corresponds to the one that generates the lowest Root Mean Square Error (RMSE) for the residual gravity anomalies ( ) (^{Tziavos, Vergos and Grigoriadis, 2009}); (^{Carrion et al., 2015}). The chosen window size for each GGMs and the corresponding RMSE is showed in the Table 1.
The effects of the residual topography on the gravity anomalies were estimated from RTM by using the TC module of the GRAVSOFT (Geodetic Gravity Field Modelling Programs) software (^{Forsberg and Tscherning, 2003}), while the gravity anomalies from GGMs were obtained from the ICGEM (International Center for Global Gravity Field Models) calculation service (^{Barthelmes and Köhler, 2016}) through the SPGG (Single-Point Global Earth Models Generator) software (^{Nicacio, 2017}). The residual gravity anomalies computed for the optimal RTM solution (Table 1) were used to estimate the disturbing potential in the EVD by the GBVP solution in the free form. The residual gravity anomalies and the corresponding histograms are shown in the Figure 3.
It should be noted that the residual gravity anomalies are obtained with the purpose of applying the remove-restore technique, so they were calculated with a satellite-only GGM and with a combined GGM with truncated series expansion. Thus, the residual gravity anomalies values can be high.
nmax | RTM solution | Moving average ratio (min) | RMSE (mGal) | |
---|---|---|---|---|
DIR-R5 | 300 | RTM17 | 17 | 25.85 |
EIGEN6C4 | 1000 | RTM5 | 5 | 14.91 |
According to ^{Torge and Müller (2012}), the spatial resolution of the GGMs is related to their maximum degree expansion, and can be calculated according to the expression:
The output DEMs (coarse elevation grid), generated by applying the moving average to the original DEM (detailed elevation grid), and according to the ratios of the Table 1, have a pixel size of 34' for the DIR-R5 solution and 10' for the EIGEN6C4 solution. This pixel sizes approaches the spatial resolution of the GGM calculated according to Erro! Fonte de referência não encontrada., 36' and 10.8' approximately.
4.3 Disturbing potential estimationFigure 5(a) Figure 5 (b)
The residual height anomaly (ζ_{RES}) at the EVD was estimated by the Fast Collocation (LSC) method (^{Bottoni and Barzaghi, 1993}) as a function of the residual gravity anomalies computed as described on the previous section. The procedure followed for this calculation can be seen in the Figure 4.
According to ^{Forsberg (1984}), the statistical properties of the gravity field variations are described by the covariance function, thus, to predict the residual gravity anomaly by Fast Collocation, it is essential to know the covariance function of the residual gravity anomalies used as input to model the gravity field.
The empirical covariance functions (ECF) were generated in the space domain with the EMPCOV program ( ^{Tscherning, 2009}) as a function of the residual gravity anomalies and fitted to the Tscherning/Rapp degree-variance model. The ECF parameters are shown in Table 2.
DIR-R5 (nmax=300) | EIGEN6C4 (nmax=1000) | |
---|---|---|
Optimal ratio (°) | 0.067 ≈ 4' | 0.067 ≈ 4' |
Signal standard deviation | 12.786 | 5.841 |
Noise standard deviation (mGal) | 5.063 | 3.830 |
The optimal ratio (also denoted sampling interval size) for the empirical covariance function corresponds with the residual gravity anomalies grid spacing (4'). The signal standard deviation is used as input for the computation of the Analytical Covariance Function (ACF), while the noise standard deviation is used for the Fast Collocation computation.
The variance for the EIGEN6C4 solution is higher than the variance for the DIR-R5 solution; this is mainly due to DIR-R5 omission errors.
To model the empirical residual gravity anomaly covariance function, is required to determine an analytical covariance function that fits the empirical values. The analytical covariance functions were generated with the COVFIT program (^{Tscherning and Knudsen, 2009}) and its parameters are showed in the Table 3.
DIR-R5 | EIGEN-6C4 | |
---|---|---|
RMSE Analytical Covariance function adjustment (mGal) | 1.0516 | 0.7193 |
Depth of Bjerhammar sphere (m) | -16684.81 | -6849.49 |
Scale factor for error variances | 0.2559 | 0.3127 |
The analytical covariance functions, on the application of the LSC method, were generated considering variations for the parameters: Bjerhammar sphere depth and scale factor for error variances (AA), and using gravity error degree variances for the GGM coefficients. These parameters were modified according to the degree of adherence of the analytical covariance functions with the empirical covariance functions. An iterative least squares procedure allows estimating these parameters by fitting the global model to the empirical covariance values for the local area. The empirical and analytical covariance functions are shown in ^{Figure 6}.
In general, adhesion between the analytical and empirical covariance functions is recorded, however, a lower adhesion was recorded for the solution given by DIR-R5 with a maximum degree expansion.
5. Results
5.1 Residual height anomaly estimation on the EVD
The estimation of ζ_{ RES } through the Fast Collocation method (ζ_{ LSC } ) (Table 4) was performed by using the generated analytical covariance functions. The height anomaly on the EVD (ζ_{ EVD } ) was computed by the restoration of the long and short wavelengths of the gravity field. This computation was accomplished by the restoration of the residual topography effects on the height anomalies (ζ_{ RTM } ) and the height anomalies from the GGMs (ζ_{ GGM } ). The zero-degree term (ζ_{ 0 } ) Erro! Fonte de referência não encontrada.was considering on the ζ_{ GGM } computation to be consistent with the GRS80 ellipsoid (^{Hofmann-Wellenhof and Moritz, 2006}):
where GM_{GGM} is the geocentric universal gravitational constant related with the GGM and GM_{GRS} is the geocentric universal gravitational constant related with the geodetic reference system (in this work, the GRS80). Table 4 shows the computed height anomalies ζ_{EVD} Erro! Fonte de referência não encontrada.and the corresponding values for ζ_{ RTM } , ζ_{ LSC } and ζ_{ GGM } .
ζ_{ LSC } (m) | ζ_{ RTM } (m) | ζ_{ GGM } (m) | ζ_{ LSC+RTM+GGM } (m) | |
---|---|---|---|---|
DIR-R5 | -0.217 | 0.3738 | 10.7305 | 10.8873 |
EIGEN6C4 | 0.027 | 0.0315 | 10.6858 | 10.7443 |
In the case of the solution given for EIGEN6C4, small quantities associated with the effects of the residual topography and the modeling of the residual height anomalies by LSC are registered.
A difference of 14 cm is registered between the two solutions. Due to the considerable magnitude of this difference, these results will be validated in future studies, for which alternative solutions will be tested considering other GGMs, degree expansion and variations in the size of the EVD adjacent zone.
5.2 Ecuadorian Vertical Datum offset
The geopotential on the EVD (W _{0} ^{i}) was estimated to obtain the Local Vertical Datum (LVD) discrepancy related to the IHRS Global Vertical Datum (W _{0}). The disturbing potential (T _{P}) was computed by the Bruns’s formula Erro! Fonte de referência não encontrada. as a function of the height anomaly (ζ_{P}) that was estimated by the free GBVP solution. Knowing the disturbing potential and the normal potential on the computation point (EVD), the geopotential is computed according to (1). Because W _{P} must be calculated referring to the mean tide system (following the IAG recommendations for the establishment of the IHRF), the two quantities on the right side of (1) should be also referred to this tide system. The normal potential U _{P} is computed considering the GRS80 ellipsoidal parameters and referred to the zero-tide system. Thus, is necessary to compute U _{ P } on the mean tide system. According to ^{Ihde et al. (2008}), the effect of the permanent tide on the potential (W _{2}) can be computed as:
where φ is the latitude of the calculation point.
Considering W _{2}, U _{P} can be expressed in the mean tide system, and thus the geopotential W _{EVD} can be computed referred also to the mean tide system. After estimating the potential value W _{P} (1), it was compared with the global vertical reference value W _{0} to calculate the vertical offset (βH _{0}) on the EVD as:
Table 5 presents the quantities involved in the vertical offset (βH _{ 0 } ) computation.
DIR-R5 (n=300) | EIGEN6C4 (n=1000) | |
---|---|---|
W _{ 0 } (IAG) (m^{2}/s^{2}) | 62636853.400 | |
U _{ 0 } (GRS80) (m^{2}/s^{2}) | 62636860.850 | |
T _{ P } (m^{2}/s^{2}) | 106.1527 | 104.7584 |
W _{ P } (m^{2}/s^{2}) | 62636850.5043 | 62636849.1101 |
U _{ P } (m^{2}/s^{2}) | 62636745.3195 | |
βH _{ 0 } (m) | 0.2873 ± 0.0274 | 0.4303 ± 0.0901 |
dW(m^{2}/s^{2}) | 2.8014 ± 0.2674 | 4.1956 ± 0.8789 |
The free solution of the GVBP, when using the GGM DIR-R5 with expansion degree 300, generates an offset of approximately 29 cm. In the case of the computed offset with the GGM EIGEN6C4 with expansion degree 1000 generates an offset of approximately 43 cm. The uncertainty in the offset computation for the EIGEN6C4 solution (± 0.0901 m) is considerably greater than the uncertainty for the DIR-R5 solution (± 0.0274 m).
In future studies, the election of the maximum degree expansion of the GGMs should consider an exhaustive analysis in terms of minimizing the GGMs commission errors. The choice of the optimal maximum degree expansion allows to reduce the noise from GGMs commission errors that are transmitted to residual functionals (e.g., gravity anomalies and gravity disturbances) computed by the remove-restore technique.
For the free solution of the GBVP, Molodensky or surface gravity anomalies were calculated. For the computation of these gravity anomalies it is required to have normal heights to reduce the normal gravity to telluroid. Because normal heights are not available for Ecuador, these were approximated by levelling heights or by heights derived from height anomalies from GGMs and ellipsoidal altitudes (GNSS), this approximation implies the loss of rigor in the application of the method.
6. Summary and conclusions
The GBVP solutions in the free form based on the Molodensky’s theory and fast collocation approach were performed considering the GGMs DIR-R5 and EIGEN6C4 for the representation of the long wavelengths of the gravitational field. The geopotential value (W _{ 0i } ) computation were carried out on the reference level corresponding to the MSL determined for the period 1988 - 2009 and materialized on the fundamental bench mark (BM03 - origin of the EVRN). With this approach, variations on the maximum degree of spherical harmonics expansion of the GGMs were tested. The vertical discrepancies related to W _{ 0 } were computed in terms of geopotential values and considering the zero-degree term (ζ_{ 0 } ) on the computation of height anomalies from the GGMs (restoration stage). The zero-degree term represents the difference of reference between the GGM and the geodetic reference system used in the offset representation (in this work, the GRS80).
Although the Molodensky gravity anomalies computation is performed without considering hypotheses of densities and gravity reductions, the heights used to reduce the theoretical gravity from the ellipsoid to the telluroid are linked to a local height system. Future studies related to the linkage of the EVD to the IHRS will take into account this aspect through the fixed solution of the GBVP.