OFFSET EVALUATION OF THE ECUADORIAN VERTICAL DATUM RELATED TO THE IHRS

Avaliação da Discrepância do Datum Vertical Equatoriano em Relação ao IHRS

Abstract

Considering the definition of the International Height Reference System (IHRS) in the geopotential space (Resolution 1/2015, International Association of Geodesy - IAG), among the present main objectives of the international geodetic community is the materialization of IHRS around the world. One fundamental task for this is the offset determination of each national vertical datum related to the IHRS. In this manuscript we establish the relationship between the Ecuadorian Vertical Datum (EVD) and the IHRS in the geopotential space following the foundations of the Resolution 1/2015 IAG. Gravity data, heights from the Ecuadorian Fundamental Vertical Network, Global Geopotential Models and Digital Elevation Models were used in the computations. Based on the Least Squares Collocation method, empirical covariance functions and spectral decomposition techniques, we realized the modelling of the geopotential in the study region (4° x 4° centered in the La Libertad tide gauge, Ecuador). Based on the referred approaches, we solved the free Geodetic Boundary Value Problem for determining the discrepancy of the EVD related to the IHRS. An offset of approximately 29 cm ± 3 cm was estimated for the W ₀ - W ₀ ⁱ relation when the GO_CONS_GCF_2_DIR_R5 model was used in the modeling of the medium and long wavelengths of the terrestrial gravity field, and approximately 43 cm ± 3 cm when the EIGEN6C4 model was used.

Keywords
International Height Reference System (IHRS); Ecuadorian Vertical Datum; Global Geopotential Models; Geodetic Boundary Value Problem; Least Squares Collocation

Resumo

A partir da definição do International Height Reference System (IHRS) no espaço do geopotencial (Resolução 1/2015 - International Association of Geodesy - IAG), a comunidade geodésica internacional tem entre seus principais objetivos a materialização do IHRS em todo o globo. Para tanto, um dos aspectos fundamentais é a determinação da discrepância de cada datum vertical nacional em relação ao IHRS. No presente trabalho estabelecemos a relação entre o Datum Vertical Equatoriano (DVE) e o IHRS, no espaço do geopotencial, seguindo aos pressupostos da Resolução 1/2015 da IAG. Dados gravimétricos, altitudes oriundas da Rede de Nivelamento Fundamental de Equador, Modelos Globais de Geopotencial e Modelos Digitais de Altitude foram utilizados nos cálculos. A modelagem do geopotencial na região de estudo (4° x 4° com centro no marégrafo La Libertad - Equador) foi realizada mediante o método de Colocação por Mínimos Quadrados, funções de covariância empíricas e técnicas de decomposição espectral. Assim, foi efetivada a solução livre do Problema de Valor de Contorno da Geodesia no DVE de forma a determinar a sua discrepância com o IHRS. Um offset de aproximadamente 29 cm ±3 cm foi estimado para a relação W ₀ - W ₀ ⁱ quando utilizado o modelo GO_CONS_GCF_2_DIR_R5 na modelagem dos médios e longos comprimentos de onda do campo de gravidade terrestre, e de aproximadamente 43 cm ±3 cm quando utilizado o EIGEN6C4.

Palavras-chave
International Height Reference System (IHRS); Datum Vertical do Equador; Modelos Globais de Geopotencial; Problema do Valor de Contorno da Geodesia; Colocação por Mínimos Quadrados

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 ₀ = 62 636 853.4 m²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 ₀ - 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:

$$W_P=U_P+T_P.$$ (1)

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:

$$U_P=U_0+\frac{\partial U_0}{\partial h}{h}_P,$$ (2)

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:

$$\delta W_{0i}=W_0-W_{0i}=W_0-\left(U_{0+}\frac{\partial U_0}{\partial h}{h}_{0i}+T_{0i}\right),$$ (3)

where the global reference value of IHRS is given by W ₀ and the geopotential value associated with the local VD is W _0i . The gradient of the normal potential ?U ₀ /?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:

$${∆g}_M=g_P-{\gamma }_{\mathrm{\Sigma }}+{\delta g}_A+{\mathrm{\Delta }g}_H.$$ (4)

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):

$$\gamma ={\gamma }_e\left(1+\beta {sin}^2\phi +\beta {\text{'}sin}^{2}2\phi \right).$$ (5)

In this equation γ_e is the normal gravity at the equator and the other involved parameters are:

$$\beta =\frac{5}{2}m-f-\frac{17}{14}mf,$$ (6)

$${\beta }^{\text{'}}=\frac{f^2}{8}-\frac{5mf}{8},$$ (7)

$$m=\frac{a{\omega }^2}{{\gamma }_e}-\frac{3}{2}{m}^2.$$ (8)

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:

$$\gamma =\frac{a{\gamma }_e{cos}^2\phi +b{\gamma }_P{sin}^2\phi }{{\left(a^2{cos}^2\phi +b^2{sin}^2\phi \right)}^{1/2}},$$ (9)

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.:

$${}_h=-\frac{2{\gamma }_e}{a}\left(1+f+m-2f{sin}^2\phi \right)h+\frac{3{\gamma }_e}{a^2}{h}^2,$$ (10)

$${\gamma }_h=\gamma -0.3086h.$$ (11)

Figure 1
Source: de Freitas, 2015

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):

$$\mathrm{\Delta }{g}_H=0.0371\left(1-3{sin}^2\phi \right).$$ (12)

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:

$$T_P=\frac{1}{2\pi }\underset{\mathrm{\Sigma }}{\overset{}{\int }}\left[\frac{\partial }{\partial n}\left(\frac{1}{l}\right)-\frac{1}{\gamma }\frac{\partial \gamma }{\partial n}\right]Td\mathrm{\Sigma },$$ (13)

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 ₀ +T ₁ +T ₂ +... in the general form:

$$T_P=\gamma =\frac{R}{4\pi }{\int }_{\sigma }^{}(g_0+g_1+g_2+\dots \left)S\right(\psi )d\sigma .$$ (14)

Which can be decomposed as:

$$T_0={}_0\gamma =\frac{R}{4\pi }{\int }_{\sigma }^{}{g}_{0}S\left(\psi \right)d\sigma$$ (15)

$$T_1={}_0\gamma =\frac{R}{4\pi }{\int }_{\sigma }^{}{g}_{0}S\left(\psi \right)d\sigma$$ (16)

where R is the mean Earth’s radius and the ensemble of terms g ₀ , g ₁ , g ₂ … 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):

$$g_0=\mathrm{\Delta }{g}_M,$$ (17)

$$g_1=\frac{R^2}{2\pi }{\int }_{\sigma }^{}\frac{h-h_P}{l^3}[g_0+\frac{3\stackrel{-}{g}}{2R}{}_0]d\sigma ,$$ (18)

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:

$$S\left(\psi \right)=\frac{1}{sin\left(\frac{\psi }{2}\right)}-6sin\left(\frac{\psi }{2}\right)+1-5cos\psi -3cos\psi *ln\left[sin\left(\frac{\psi }{2}\right)+{sin}^2\left(\frac{\psi }{2}\right)\right].$$ (19)

The spherical distance can be computed based on the coordinates of the referred points (ϕ, λ) and (ϕ’, λ’) respectively as (Gemael, 2012, p.146):

$$\psi =arc.cos\left(sin\phi *sin{\phi }^{\text{'}}+cos\phi *cos{\phi }^{\text{'}}*cos\mathrm{\Delta }\lambda \right).$$ (20)

The Molodensky’s solution is applied for a geocentric GRS. For an arbitrary GRS, it is rewritten as:

$$T_P=\frac{\mathrm{\Delta }GM}{R}+\frac{R}{4\pi }{\int }_{\sigma }^{}(g_0+g_1+g_2+\dots \left)S\right(\psi )d\sigma ,$$ (21)

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:

$${{\mathrm{\Delta }\mathrm{g}}_M}_{RES}={{\mathrm{\Delta }\mathrm{g}}_M}_{OBS}-{{\mathrm{\Delta }\mathrm{g}}_M}_{GGM}-{{\mathrm{\Delta }\mathrm{g}}_M}_{RTM}.$$ (22)

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:

$${{\mathrm{\Delta }\mathrm{g}}_M}_{RTM}=G\rho \underset{-\infty }{\overset{\infty }{\iint }}\underset{z=h_{REF}\left(x,y\right)}{\overset{z=h\left(x,y\right)}{\int }}\frac{z-h_P}{{\left[{\left(x_Q-x_P\right)}^2+{\left(y_Q-y_P\right)}^2+{\left(z_Q-z_P\right)}^2\right]}^{3/2}}{dx}_Q{dy}_Q{dz}_Q,$$ (23)

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:

$$T=T_{RES}+T_{GGM}+T_{RTM}.$$ (24)

Analogously, the height anomaly is compute as:

$$={}_{RES}+{}_{GGM}+{}_{RTM}.$$ (25)

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).

$$T_P={\zeta }_P{\gamma }_0$$ (26)

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.

Figure 2:
Gravity data distribution in the 4° x 4° delimited region around the EVD

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 ($\Delta g_{M_{RTM}}$ ) omitted by the GGMs, which are used to compute the residual gravity anomalies according to the equation Erro! Fonte de referência não encontrada..

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 ₀ ⁱ ) 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.

Figure 3:
(a) DIR-R5 (nmax=300); (b) EIGEN6C4 (nmax=1000)

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.

Table 1:
Moving average ratio. RTM optimal solution

	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:

$${res}_{GGM}=\frac{{180}^{°}}{n_{max}}\left[°\right],or{res}_{GGM}=\frac{20000km}{n_{max}}\left[km\right].$$ (27)

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.

Figure 4:
Residual height anomaly estimation

Figure 5:
Residual gravity anomalies 4'x4' grid.

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.

Table 2:
Empirical covariance function parameters

	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.

Table 3:
Parameters of the analytical covariance functions

	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.

Figure 6:
Empirical and analytical covariance functions.

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 (ζ₀ ) 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):

$${}_0=\frac{{GM}_{GGM}-{GM}_{GRS}}{r\gamma },$$ (28)

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 .

Table 4:
Height anomaly on the EVD

	ζ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 ₀ ⁱ) was estimated to obtain the Local Vertical Datum (LVD) discrepancy related to the IHRS Global Vertical Datum (W ₀). 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 ₂) can be computed as:

$$W_2\left(\phi \right)=0.9722-2.8841{sin}^2\phi -0.0195{sin}^4\phi \left[m^2{s}^{-2}\right],$$ (29)

where φ is the latitude of the calculation point.

Considering W ₂, 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 ₀ to calculate the vertical offset (βH ₀) on the EVD as:

$$\beta H_0=\frac{W_0-W_P}{{\gamma }_0}.$$ (30)

Table 5 presents the quantities involved in the vertical offset (βH ₀ ) computation.

Table 5:
Vertical offset at the EVD

	DIR-R5 (n=300)	EIGEN6C4 (n=1000)
W0 (IAG) (m2/s2)	62636853.400
U0 (GRS80) (m2/s2)	62636860.850
TP (m2/s2)	106.1527	104.7584
WP (m2/s2)	62636850.5043	62636849.1101
UP (m2/s2)	62636745.3195
βH0 (m)	0.2873 ± 0.0274	0.4303 ± 0.0901
dW(m2/s2)	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 ₀ were computed in terms of geopotential values and considering the zero-degree term (ζ₀ ) 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.

Acknowledgement

The authors would like to acknowledge the SENESCYT, for the José Carrión scholarship, and CNPq (Grant Process No. 306936/2015-1) by the financial support. Also, the Federal University of Parana, the Politecnico Di Milano, and the IGM of Ecuador by the institutional support.

References

• Amjadiparvar, B., Rangelova, E. and Sideris, M. G. (2016) “The GBVP approach for vertical datum unification: recent results in North America”, Journal of Geodesy, 90(1), p. 45-63. doi: 10.1007/s00190-015-0855-8.
  » https://doi.org/10.1007/s00190-015-0855-8.
• Andersen, O. B. and Knudsen, P. (2016) “Deriving the DTU15 Global high resolution marine gravity field from satellite altimetry”.
• BARTHELMES, F. and KÖHLER, W. (2016) “International Centre for Global Earth Models (ICGEM”, in The Geodesists Handbook, Journal of Geodesy 2016, p. 1177-1180. doi: 10.1007/s00190-016-0948-z.
  » https://doi.org/10.1007/s00190-016-0948-z.
• Barzaghi, R. (2016) “The Remove-Restore Method”, in Grafarend, E. (org.) Encyclopedia of Geodesy Cham: Springer International Publishing, p. 1-4. doi: 10.1007/978-3-319-02370-0_19-1.
  » https://doi.org/10.1007/978-3-319-02370-0_19-1
• Bottoni, G. P. and Barzaghi, R. (1993) “Fast collocation”, Bulletin Geodésique, 67(2), p. 119-126. doi: 10.1007/BF01371375.
  » https://doi.org/10.1007/BF01371375.
• Bruinsma, S. et al (2013) “The new ESA satellite-only gravity field model via the direct approach”, Geophys. Res. Lett., 40(14), p. 3607-3612. doi: doi: 10.1002/grl.50716.
  » https://doi.org/doi: 10.1002/grl.50716.
• Carrión, D. et al (2015) “Assessing the GOCE models accuracy in the Mediterranean area”, Newton’s Bulletin, (5), p. 63-82.
• De Freitas, S.R.C. (2015) “SIRGAS-WGIII activities for unifying height systems in Latin America”, Revista Cartográfica Pan-American Institute og Geography and History, 91(91), p. 75-92.
• Forsberg, R. (1984) A Study of Terrain Reductions, Density Anomalies and Geophysical Inversion Methods in Gravity Field Modelling Ohio State University, Department of Geodetic Science and Surveying (A Study of Terrain Reductions, Density Anomalies and Geophysical Inversion Methods in Gravity Field Modelling).
• Forsberg, R. (2008). Terrain Effects in Geoid Computations. In: Lecture notes, international school for the determination and use of the geoid, Como, 15-19 Sept 2008.
• Forsberg, R. and Tscherning, C. (2003) “Geodetic Gravity Field Modelling Programs”. Copenhagen.
• Forsberg, R. and Tscherning, C. C. (1981) “The use of height data in gravity field approximation by collocation”, Journal of Geophysical Research: Solid Earth, 86(B9), p. 7843-7854. doi: 10.1029/JB086iB09p07843.
  » https://doi.org/10.1029/JB086iB09p07843.
• Forsberg, R. and Tscherning, C. C. (1997) “Topographic effects in gravity field modelling for BVP”, in Sansó, F. and Rummel, R. (orgs.) Geodetic Boundary Value Problems in View of the One Centimeter Geoid Berlin, Heidelberg: Springer Berlin Heidelberg, p. 239-272. doi: 10.1007/BFb0011707.
  » https://doi.org/10.1007/BFb0011707.
• Förste, C. et al (2014) “EIGEN-6C4 - The latest combined global gravity field model including GOCE data up to degree and order 1949 of GFZ Potsdam and GRGS Toulouse”, EGU General Assembly, 16, p. 3707. doi: 10.5880/icgem.2015.1.
  » https://doi.org/10.5880/icgem.2015.1
• Gatti, A., Reguzzoni, M. and Venuti, G. (2012) “The height datum problem and the role of satellite gravity models”, Journal of Geodesy, 87(1), p. 15-22. doi: 10.1007/s00190-012-0574-3.
  » https://doi.org/10.1007/s00190-012-0574-3.
• Gemael, C. (2012) "Introdução à geodésia física". Federal University of Parana Curitiba.
• Heck, B. and Rummel, R. (1990) “Strategies for solving the vertical datum problem using terrestrial and satellite geodetic data”, in Sea Surface Topography and the Geoid Springer, p. 116-128.
• Heikkinen, M. (1979) “On the Honkasalo term in tidal corrections to gravimetric observations”, Bulletin Geodésique, 53(3), p. 239-245. doi: 10.1007/BF02523955.
  » https://doi.org/10.1007/BF02523955
• Hinze, W. J. et al (2005) “New standards for reducing gravity data: The North American gravity database”, Geophysics, 70(4), p. J25-J32. doi: 10.1190/1.1988183.
  » https://doi.org/10.1190/1.1988183.
• Hirt, C., Featherstone, W. E. and Marti, U. (2010) “Combining EGM2008 and SRTM/DTM2006.0 residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data”, Journal of Geodesy, 84(9), p. 557-567. doi: 10.1007/s00190-010-0395-1.
  » https://doi.org/10.1007/s00190-010-0395-1.
• Hofmann-Wellenhof, B. and Moritz, H. (2006) Physical geodesy Springer Science & Business Media.
• Honkasalo, T. (1964) “On the tidal gravity correction”, Bulletin of Theoretical and Applied Geophysics, 6, p. 34-36.
• Ecuadorian Military Geographic Institute (2013) Geodetic Control Point Quito.
• Ihde, J., Mäkinen, J. and Sacher, M. (2008) “Conventions for the Definition and Realization of a European Vertical Reference System (EVRS)-EVRS Conventions 2007”, EVRS Conv., p. 1-10.
• Ihde, J., Sánchez, L. and Sideris, M. (2010) Theme 1 : Global Unified Height System Introductory presentation Miami: IAG.
• Kutterer, H. and Neilan, R. (2015) “International Association of Geodesy”, in IAG Travaux Volume 39 Prague, p. 375-379.
• Lehmann, R. (2000) “Altimetry--gravimetry problems with free vertical datum”, Journal of Geodesy, 74(3), p. 327-334. doi: 10.1007/s001900050290.
  » https://doi.org/10.1007/s001900050290.
• Nicacio, E. (2017) “Single Point Global Earth Models Generator”. Curitiba. Available at: www.ciencias geodesicas.ufpr.br/spgg/
  » www.ciencias geodesicas.ufpr.br/spgg/
• Olson, C. J., Becker, J. J. and Sandwell, D. T. (2014) “A new global bathymetry map at 15 arcsecond resolution for resolving seafloor fabric: SRTM15_PLUS”, in AGU Fall Meeting Abstracts, p. 3.
• Pavlis, N. K. et al (2008) “An earth gravitational model to degree 2160: EGM2008”, presented at the 2008 General Assembly of the European Geosciences Union, Vienna, Austria , April 13-18, 84(1), p. 2-4.
• Plag, H. and Pearlman, M. (2009) Global Geodetic Observing System. Meeting the Requirements of a Global Society on a Changing Planet in 2020 Organizado por H. Plag e M. Pearlman. London: Springer.
• Rummel, R. (2000) “Global Integrated Geodetic and Geodynamic Observing System (GIGGOS)”, in Rummel, R. et al. (orgs.) Towards an Integrated Global Geodetic Observing System (IGGOS): IAG Section II Symposium Munich, October 5-9, 1998 Berlin, Heidelberg: Springer Berlin Heidelberg, p. 253-260. doi: 10.1007/978-3-642-59745-9_53.
  » https://doi.org/10.1007/978-3-642-59745-9_53
• Sansò, F. and Sideris, M. G. (2013) “The Forward Modelling of the Gravity Field”, in Sansò, F. and Sideris, M. G. (orgs.) Geoid Determination: Theory and Methods Berlin, Heidelberg: Springer Berlin Heidelberg , p. 3-71. doi: 10.1007/978-3-540-74700-0_1.
  » https://doi.org/10.1007/978-3-540-74700-0_1
• Torge, W. and Müller, J. (2012) Geodesy, De Gruyter Textbook Berlín: De Gruyter.
• Tscherning, C. C. and Barzaghi, R. (1991) “FASTCOLC PROGRAM”. Module of the GRAVSOT software package used to apply the Fast Collocation Method.
• Tscherning, C.C. (2009) “EMPCOV PROGRAM”. Module of the GRAVSOT software package used to determinate Empirical Covariance Functions. Ohio.
• Tscherning, C.C. and Knudsen P. (2009) “COVFIT PROGRAM”. Module of the GRAVSOT software package used to determinate of Analytical Covariance Functions. Ohio.
• Tziavos, I. N. and Sideris, M. G. (2013) “Topographic Reductions in Gravity and Geoid Modeling”, in Sansò, F. and Sideris, M. G. (orgs.) Geoid Determination: Theory and Methods Berlin, Heidelberg: Springer Berlin Heidelberg , p. 337-400. doi: 10.1007/978-3-540-74700-0_8.
  » https://doi.org/10.1007/978-3-540-74700-0_8
• Tziavos, I. N., Vergos, G. S. and Grigoriadis, V. N. (2009) “Investigation of topographic reductions and aliasing effects on gravity and the geoid over Greece based on various digital terrain models”, Surveys in Geophysics, 31(1), p. 23. doi: 10.1007/s10712-009-9085-z.
  » https://doi.org/10.1007/s10712-009-9085-z.
• Uotila, U. A. (1980) “Note to the users of International Gravity Standardization net 1971”, Bulletin geodésique Springer, 54(3), p. 407-408.
• Werner, M. (2001) “Shuttle Radar Topography Mission (SRTM) mission overview”, Frequenz, 55(3-4), p. 75-79.

Publication Dates

• Publication in this collection
  Oct-Dec 2018
• Date of issue
  Dec 2018

History

• Received
  28 Mar 2018
• Accepted
  27 July 2018