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

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


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, Bulletin of Geodetic Sciences, 24(4): 503-524, Oct-Dec, 2018 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.

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:   =   +   . (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: 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 freeair 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. ( 12)).This kind of anomaly is given by: ∆  =   −  Σ +   + Δ  . (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  (5) (Hofmann-Wellenhof and Moritz, 2006): In this equation  e is the normal gravity at the equator and the other involved parameters are: Bulletin of Geodetic Sciences, 24(4): 503-524, Oct-Dec, 2018 (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: 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 (10); 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 (11): ℎ =  − 0.3086ℎ.
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 (4), the Honkasalo term is removed by adding the latitudinal-dependent Honkasalo correction Δg H (12) in mGal (Hinze et al., 2005): Figure 1: Reference surfaces for heights in the geopotential space.The normal height is the distance Q 0 Q along the normal plumb line from the level ellipsoid to the telluroid , or the distance Q'P from the quasigeoid to the point P. Note that the normal plumb line, on which these distances are determined, is slightly curved due to the non-parallels of the level surfaces.
The set of points Q which accomplish the property U Q =W P form the surface .The distance  from the telluroid to the surface is called height anomaly.The quasigeoid and the geoid are respectively the surfaces determined by the corresponding height anomalies  and geoid heights N over the level ellipsoid.

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, 2002, pp. 270): where, h and h P are the heights of the integrating and computing points respectively and ̅ is the gravity mean value.
In (19)  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: Geodetic Sciences, 24(4): 503-524 Oct-Dec, 2018 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: Δg   = Δg   − Δg   − Δg   . ( 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 ( 13) is solved by using only the small residual Molodensky's gravity anomalies (22).Then it is possible to compute T RES , and the completely disturbing potential is obtained in a composition process by:  =   +   +   . (24) Analogously, the height anomaly is compute as: The equations ( 24) and ( 25) are related by the Bruns formula (26), and from ( 24) is possible to compute the geopotential (W P ) by using the equation ( 1). ( Bulletin of Geodetic Sciences, 24(4): 503-524, Oct-Dec, 2018

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, freeair 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.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.

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 (4).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 ( RTM M g  ) omitted by the GGMs, which are used to compute the residual gravity anomalies according to the equation ( 22).
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 ( 25) . The disturbing potential (T P ) is compute from  P through the Bruns's formula (26), the geopotential value (W 0 i ) is computed as function of the disturbing potential from (1), and the EVD offset is computed according to (3).Details about these computations are described in the following subsections.

Gravity anomalies computation and outlier elimination
Gravity anomalies are computed from direct gravity observations by (4), 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.

Residual gravity anomalies
Residual gravity anomalies were computed for all the available gravity records according to ( 22).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 coveragebathymetry); 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 ( 22) 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 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.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 (27), 36' and 10.8' approximately.

Disturbing potential estimation
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.The Fast Collocation method is applied by the FASTCOLC program (Tscherning and Barzaghi, 1991) of the GRAVSOFT software, and for this it is required that the available gravimetric information must be interpolated in a regular grid.The generated grids, with spacing of 4 arc minutes, corresponding to each of the proposed solutions, were generated by the weighted average interpolation method using the GEOGRID program of the GRAVSOFT software (Forsberg and Tscherning, 2003).In 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.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.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.

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 ) ( 28 ) 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 ( 28 ) and the corresponding values for  RTM ,  LSC and  GGM .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.

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 (26) 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: 0 =  0 −    0 . ( Table 5 presents the quantities involved in the vertical offset (βH 0 ) computation.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.

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.

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

FigureFigure 5 :
Figure 5 (a) and Figure 5 (b) the residual gravity anomalies grids are shown (for each solution) and their frequency distribution (histograms).

Table 1 :
Moving average ratio.RTM optimal solution

Table 3 :
Parameters of the analytical covariance functions

Table 4 :
Height anomaly on the EVD

Table 5 :
Vertical offset at the EVD