On the connection of the Ecuadorian Vertical Datum to the IHRS

2.1 Additional remarks about HRS

According to IAG Resolution #1 – 2015, the definition of the IHRS is given in terms of geopotential parameters, the vertical coordinates being geopotential numbers referred to an equipotential surface of the earth’s gravity field realized by a global conventional value W ₀, computed as detailed in Sánchez et al. (2016). The establishment of the IHRS meets the requirements of GGOS, which according to its terms of reference must: (1) support a precise (at centimetric level) combination of physical and geometric heights on a global scale; (2) enable the unification of existing LVDs; and (3) guarantee vertical coordinates with global consistency (the same accuracy everywhere) and long-term stability (the same order of accuracy at all times) (Sánchez and Sideris 2017). The realization of the IHRS must be carried out in terms of geopotential values and be established according to well documented conventions. The coordinates in the IHRF are the geopotential numbers: the transformation of the geopotential numbers into physical heights and the geometric link to the reference surface are a matter of the IHRF realization and not of definition. As already pointed out, in order to establish the IHRS, it is necessary to adopt a conventional W ₀, estimate the local values W _0i , and determine their differences according to equation (1) (Ihde et al. 2010).

2.2 Fixed GBVP

Potential boundary value problems (BVPs) are applied in physical geodesy to the determination of the gravitational potential V outside the earth masses (which is a harmonic function). If the densities within the planet and its boundary were known, the potential of the Earth’s gravitational field could be found through an integration over the total Earth’s volume υ (Torge and Müller 2012).

where G is the universal gravitational constant, ρ is the density of topographic masses, and l the distance to the mass element and the attracted point (Hofmann-Wellenhof and Moritz 2006).

However, the densities inside the earth are in general poorly known and this formula (equation (2)) cannot give an accurate estimate of the gravitational potential. So, usually, the external gravitational potential is obtained by means of the potential theory and considering the terrestrial surface as boundary. Several GBVPs, which are based on the BVPs of the potential theory, have been developed. As an example, one can consider the GBVP problem that is used in classical geoid determination, which is based on the knowledge of gravity anomalies Δg _geoid = g _geoid – γ _ellipsoid (clearly dependent on reductions for the effect of the topographic masses above the geoid). Another very well-known GBVP is the Molodensky’s problem (Molodensky et al. 1962) which allows finding the solution without topographic reductions for the Earth’s crust gravity effect.

The Molodensky’s theory is the basis of the modern vision for solving the GBVP by applying the boundary condition on the physical surface of the earth related to the normal derivative of the so-called disturbing potential T (difference between the true earth’s gravity potential (W) and modeled theoretical ones (U)).

Since the publication of the basic works of Stokes (1849) and Molodensky et al. (1962), the determination of the external gravity field of the earth from gravity observations performed on the earth’s surface has been related to the solution of a free BVP of potential theory (Heck 2011).

The advent of accurate satellite positioning techniques (GNSS) and the consequent knowledge of earth’s surface geometry makes it possible to replace the free GBVP by the fixed one. Thus, the fixed GBVP solution is based on the determination of the geopotential W (x, y, z) in the space outside the terrestrial surface (Heck 2011). In the fixed solution, gravity disturbances (δg) are used instead of gravity anomalies (Δg), and are the input for the disturbing potential computation.

The terms in equation (4) are obtained by solving the Hotine’s integral as (Hotine 1969) follows:

where R is the mean earth’s radius and σ is the unit sphere. The term u _i is associated with local corrections and criteria for linearization. Usually, the computation methods involve only the first two corrections given by equations (6) and (7) (Hofmann-Wellenhof and Moritz 2006).

where h and h _P are the heights of the integrating and computing points, respectively (referred to the ITRF2008, epoch 2016.43), and l 0 is given by

The integral kernel H(Ψ) in equation (5) has the following form (Hofmann-Wellenhof and Moritz 2006):

where Ψ is the angular spherical distance between the computation and integration point.

2.3 Remove–compute–restore technique

According to Forsberg (1997), the remove–compute–restore technique allows a frequency analysis of the data. The long wavelength component of the data is modeled via a GGM. The gravimetric signal that dominates the short wavelengths of the gravity field, derived from the gravitational effect of the topographic masses, can be used to smooth the gravitational field before performing the modeling processes. The long and the short wavelength components are removed from the observed Q quantity by subtracting the corresponding terms (Q _RTM, Q _GGM) according to the below expression:

In this approach, the effect of topography on gravity (Q _RTM) is obtained as the effect of the topographic masses in between the actual digital terrain model (DTM) and a mean long wavelength DTM surface corresponding to the maximum degree of the GGM used to compute the gravity effect of the long and very long wavelengths (Q _GGM). This residual terrain effect in planar approximation (Forsberg and Tscherning 2008) can be expressed as follows:

where G is the universal gravitational constant, ρ is the standard density of the topographic masses, h _P is the height of the computation point given by a DTM, z is the height relative to the reference surface, h _ref is the mean long wavelength of DTM, [x _P, y _P, z _P] are the computation point coordinates, and [x _Q, y _Q, z _Q] are the integration point coordinates (Forsberg and Tscherning 1997).

On the other hand, the Q _GGM quantity, corresponding to the long and medium wavelength of the gravity field, is computed from a GGM which is given as a truncated spherical harmonic series. The residual Q _r component is then used in the “compute” step which can be based on, for e.g., LSC.

2.4 LSC and least squares fast collocation (LSFC)

The LSFC method, developed by Bottoni and Barzaghi (1993), is a variant of LSC, aimed at reducing computing time and required computational resources. The Fast Collocation method involves imposing constraints on input data. In case homogeneous and gridded data are considered, the covariance matrix has a particular structure that allows a fast computation of the collocation solution. If a two-dimensional data grid is considered and a covariance function that is only dependent on the plane distance d _PQ between points P and Q

one can prove that the covariance matrix C computed for this grid using the covariance function of equation (12) is a symmetric block Toeplitz matrix and that each block has a Toeplitz structure. Given this highly regular structure, efficient algorithms for storing the matrix and for solving the collocation formula can be found.

2.5 Study region

The study has been carried out in a region centered on the EVD (Figure 1), more precisely in an area centered on BM03, i.e., the reference benchmark linked to the tide gauge named La Libertad. The study region is delimited by a square of 4° and having as centroid the EVD. The region between latitudes 0°13′10.1172″S–4°13′10.1172″S, and longitudes 78°54′19.4652″W−82°54′19.4652″W, is characterized by an irregular topography/bathymetry which ranges approximately from −4,698 to 1,552 m, involving the ocean – continent transition zone. The region shows high seismic activity, mainly due to the subduction of the Nazca Plate under the South American Plate; additionally, deformations of continental geological structures generate superficial earthquakes (Yepes et al. 1994). Historically, earthquakes of great magnitude have been recorded in Ecuador. The country is located in the Pacific “Circle of Fire,” the area with the largest seismic activity on the planet (Gonzales et al. 1988). Climatic phenomena are periodically related to the influence of ocean currents, among which we have the equatorial currents, Humboldt current, and the El Niño phenomena. These phenomena generally produce variations in temperature and, consequently, variations in the density of oceanic waters due to thermal expansion/contraction.

Figure 1

EVD location (La Libertad).

As said, the EVD discrepancy computation in relation to the IHRS was performed considering the BM03 benchmark (Figure 1) as the computation point. The coordinates of BM03 refer to the observation made by the Geographical Army Institute of Ecuador (IGM-EC) prior to the earthquake that occurred on 16th April 2016 and are in the SIRGAS 95 reference epoch 1995.4 (IGM-EC 2013). BM03 was affected by the aforementioned earthquake, the epicenter of which was located near the EVD, with magnitude 7.8 on the Richter seismological scale.

The computations developed in this study were performed considering the pre-seismic positions, and therefore, the displacement in the VD due to the seismic event must be considered for future computations and applications. We found the existence of several MSL determinations carried out by Ecuadorian Navy Oceanographic Institute (INOCAR) at the La Libertad tide gauge station since its installation. Also, according to INOCAR reports, BM48 is the origin point benchmark of the Ecuadorian Vertical Reference Network (EVRN). However, currently the EVRN origin point considered by the IGM-EC refers to the MSL value defined for the tide gauge observation period 1988–2009 and materialized in the BM03 benchmark.

2.6 Gravity dataset

Due to the non-availability of a unique and consistent geodetic database containing the information required for the development of this work, it was necessary to acquire gravity data from different sources. The diversity of data sources also implies the existence of records with heterogeneous characteristics, mainly in terms of observation methods and epochs, spatial distribution, associated precisions, and equipment used.

The main features of the geodetic data used in this work are detailed in the following. In total, 4,808 points with terrestrial gravimetric information (Figure 2b) were compiled within the region of 4° × 4° centered on the EDV point.

Figure 2

(a) Ocean gravity data from BGI database (black), and from SGGSA (red). (b) Terrestrial gravity data. (c) Airborne gravity records. (d) Gravity anomalies from WGM2012 model. (e) Gravity anomalies from DTU15 model.

In addition to the observed gravity, the records provide information about latitude, longitude, and leveled heights. Terrestrial gravity data records have been collected for geodetic purposes by the IGM-EC (Ecuadorian Military Geographic Institute) since the 1960s. In this group, gravity records were mostly observed on the principal highways mainly along with leveling lines. Some data refer to gravimetric densification campaigns. The ocean gravity data (3,209 records in black in Figure 2a) were made available by the Bureau Gravimétrique International (BGI) (2017). A smaller portion of the records located near the coast (321 records in red Figure 2a) was made available by the Subcommittee of Gravimetry and Geoid for South America (SGGSA) and refer to surveys for exploratory geophysics. The ocean gravity records of the BGI database, which are within the EVD contiguous region (4° × 4°), belong to 21 surveys carried out between 1957 and 1987. For each BGI record, the latitude, longitude, observed gravity, Free-Air anomaly, and Bouguer anomaly are known in the GRS67 (Geodetic Reference System 1967).

The airborne gravity records (3,443 points) contained in the 4° × 4° region in the vicinity of the EVD (Figure 2c) belong to the gravity data base of the SGGSA. The airborne gravity records are reduced to the terrain level according to WGS84 and EGM96 geoid, referred to the IGSN71 (International Gravity Standarization Net 1971) gravity datum and have information related to latitude, longitude, and height derived from a DEM. Gravity anomalies of the DTU15 model (Andersen and Knudsen 2016) derived from satellite altimetry (Figure 2e) and located in the oceanic part of the study region are also used as a source of complementary gravity information for modeling the earth’s gravity field. Finally, gravity values were also obtained from the World Gravity Map2012 (WGM2012) (Bonvalot et al. 2012). This is a set of maps of gravity anomalies and digital grids computed on a global scale based on gravity and height models. The WGM2012 model consists of a set of three maps of gravity anomalies (surface free-air anomalies, Bouguer anomalies, and isostatic anomalies) derived from the GGM EGM2008 and fill-in gravity information by the DEM ETOPO1. WGM2012 is the first set of maps of gravity anomalies that considers the contribution of most masses present on the surface of the planet (atmosphere, oceans, inland or continental seas, lakes, ice caps, etc.). The gravity anomalies of the WGM2012 model (Figure 2d) are obtained by rigorous computation, consistent with geodetic and geophysical definitions. The model provides homogeneous gravity information of the terrestrial gravity field on a regional and global scale. Since the gravimetric records come from different sources, it is necessary to consider their heterogeneity. According to Moritz (1980), LSC is a method that allows determining the anomalous gravity field by combining geodetic observations of different kinds (of different natures or with different spectral resolutions and precisions), through the efficient estimation of the gravitational field, by using the statistical characteristics of heterogeneous data in the form of covariance functions.

2.7 Gravity disturbance computation

Gravity disturbances were used in the computation in order to avoid local VD effects on geopotential modeling associated with reductions related to simplifying assumptions of the crust structure, in addition to the consideration of heights in local references. The processing performed on different data sources is detailed in the following sub-sections.

2.8 Outlier filtering

Elimination of the ocean and terrestrial gravity outliers was based on the comparison of the observed and modeled gravity disturbances. The EIGEN6C4 GGM (considering its maximum degree in spherical harmonics expansion, n = 2,190), which in the working area has an accuracy of approximately 5 cm (determined by comparison with GNSS/lev records), was used as the basis for outliers filtering. The effect of residual terrain modeling (RTM) (Forsberg 1984) on EIGEN6C4 reduced the gravity disturbances and was considered to improve outlier detection. Because the EIGEN6C4 information reproduced the long and very long wavelengths of the gravity field, to be comparable with the observed gravity disturbances, they must also include the contribution of the short wavelengths of the gravity field obtained by RTM. The reference DTM surface for the RTM computation was determined by smoothing the maximum resolution DEM (SRTM1 – spatial resolution of 1 arc sec (Farr et al. 2007)) by applying a low-pass filter. The low-pass filtering was performed by computing the moving average having a size which is defined through statistical analysis of the root mean square (RMS) of the residuals obtained for different versions of the smoothed DEM (Tziavos et al. 2009; Carrion et al. 2015). Once the best cap size (2′) is defined, outliers in the gravity dataset were selected by applying the three-sigma criterion (3σ) based on the δg _res statistics.

where δg _obs is the gravity disturbances obtained from the observed gravity records, δg _EIGEN6C4 is the gravity disturbances derived from the EIGEN6C4, and δg _RTM is the contribution of the residual topography on the gravity disturbances.

The elimination of outliers was performed independently for each data subset. Gravity disturbances used to compute residual gravity disturbances (δg _res) are presented in Figure 3.

Figure 3

Gravity disturbances after elimination of outliers.

The statistics for the sets of gravity disturbances (δg) after outlier elimination, presented in Table 1, show the characteristics of the gravity disturbances for each of the gravity datasets.

2.9 Residual gravity disturbances computation

Residual gravity disturbances (δg _res) were computed considering the contribution of the GO_CONS_GCF_2_DIR_R5 GGM, and the related RTM. The δg _{GOCONS_GCF_2_DIR_R5} gravity disturbances (Bruinsma et al. 2013) were computed for their maximum spherical harmonic expansion (n = 300) and also considering truncation up to n = 200, after evaluating the performance of various combinations of GO_CONS_GCF_2_DIR_R5 (truncated versions) with the EGM2008 GGM. For both solutions, the zero-degree term was included in the computations, in order to consider the difference between the values employed by the GGM and reference ellipsoid for the geocentric gravitational constant GM (Hofmann-Wellenhof and Moritz 2006). The process followed for the computation of residual gravity disturbance values (δg _res) is shown in the flowchart in Figure 4.

Figure 4

Residual gravity disturbances computation.

The maximum expansion degree used to compute the gravity disturbances in the truncated GGM was experimentally defined by analyzing the performance, in terms of height anomalies, of several combinations of GO_CONS_GCF_2_DIR_R5 (varying its maximum degree) with EGM2008 (with an average accuracy of 6 cm, estimated by evaluating geoidal heights in GNSS/level records).

Height anomalies from the combined models were compared with values derived from GNSS/lev records, and the combined GGMs performance was evaluated according to the RMS for the differences ζ _GNSS/lev −ζ _GGM.

In Figure 5, it can be seen that the combined GGM with better performance (RMS = 0.2459 m) corresponds to that established by considering n = 200 as the maximum degree for the GO_CONS_GCF_2_DIR_R5. From this analysis it can be concluded that the satellite-only model GO_CONS_GCF_2_DIR_R5, n = 200, optimizes the performance of the EGM2008 which for the evaluated dataset generates an RMS of 0.3007 m. Because GO_CONS_GCF_2_DIR_R5 is a satellite-only model, the functionals computed from its coefficients are not affected by any LVD offsets. The next step was the computation of the RTM effect (according to equation (11)) and its reduction. To compute the effects of the residual topography on the gravity disturbances ( δ g RTM ), the DEM fusion SRTM1 (1″ of spatial resolution for the continental region) and SRTM15 (15″ of spatial resolution for the ocean region – bathymetry) was used.

Figure 5

GO_CONS_GCF_2_DIR_R5/EGM2008 combinations performance.

The variation in the RMS as a function of the window radius used for the low-pass filter to generate the smoothed DEM, according to the method proposed by Carrion et al. (2015), implied by the GO_CONS_GCF_2_DIR_R5 model can be seen in Figure 6. The optimal RTM solutions (radius = 18 ′ and radius = 21 ′ ) were then selected for getting the final δ g res , computed according to the following expression:

where δg is the gravity disturbances data collected in the study, δg _GGM is the gravity disturbances from the GO_CONS_GCF_2_DIR_R5 GGM, and δg _RTM is the residual topography effect on the gravity disturbances.

Figure 6

GO_CONS_GCF_2_DIR_R5: RMS vs moving average ratio for (a) n = 300 and (b) n = 200.

The statistics for the calculated residual gravity disturbances, shown in Table 2, denote a better adherence of the calculated disturbances with those obtained from the GGM with n = 300, and considering the RTM effect, which is also appreciated in the distribution of frequencies presented in the histograms of Figures 7 and 8.

Figure 7

GO_CONS_GCF_2_DIR_R5 (n = 200): residual gravity disturbances.

Figure 8

GO_CONS_GCF_2_DIR_R5 (n = 300): residual gravity disturbances.

Figures 7 and 8 show the behavior and spatial distribution of the residual gravimetry disturbances in the study region.

2.10 Adjustment by LSFC

The residual height anomaly (ζ _RES) in the EVD point is estimated as a function of the δg _RES and by the LSFC method (Bottoni and Barzaghi 1993). Figure 9 shows the procedure followed for the ζ _RES computation in the EVD. For the application of the LSFC, a 4-min-spacing grid of δg _RES was estimated. The interpolation method used in the grid generation was the inverse distance weighting, which was carried out using the GEOGRID program of the GRAVSOFT package (Forsberg and Tscherning 2008).

Figure 9

Residual height anomaly estimate through the fixed GVBP.

As explained by Forsberg and Tscherning (2008), GRAVSOFT is a software for regional and local gravity modeling, and consists of a suite of Fortran programs to tackle many different problems of physical geodesy. The programs have been developed since the early 1970’s first at the Geodetic Institute, and later at National Survey and Cadastre of Denmark (now DTU-Space) and the Geophysical Institute (now Geophysics Dept. of the Niels Bohr Institute), University of Copenhagen. For carrying out this work, the GRAVSOFT programs were provided by the Department of Civil and Environmental Engineering of the Politecnico di Milano (POLIMI).

Because of the heterogeneous spatial data distribution, it is not appropriate to determine the grid spacing according to the criterion of the average interpoint distance (Tocho 2006). Therefore, the 4-min-spacing grid was determined experimentally, considering the performance of the LSFC on the residual height anomalies estimation. Two grids of residuals of gravity disturbances were computed, one based on residuals obtained using the GGM GO_CONS_GCF_2_DIR_R5 at full resolution (i.e., d/o 300) and the other based on residuals obtained with the same model at d/o 200. The statistics of these gridded δg _RES are presented in Table 3. These statistics show the characteristics of the gravity data that will be used to estimate the residual height anomaly on the calculation point, using the fast collocation method.

The residual anomaly grids estimated using the two residual gridded gravity disturbances are shown in Figures 10 and 11 for n = 200 and n = 300, respectively.

Figure 10

Residual gravity disturbances (n = 200).

Figure 11

Residual gravity disturbances (n = 300).

The residual gravity disturbances empirical covariance functions (ECF) were estimated for each solution, by using the EMPCOV program (Tscherning 2009) of the GRAVSOFT computational package (Forsberg and Tscherning 2008). The ECF parameters are presented in Table 4.

The empirical and analytical covariance functions for both solutions are presented in Figure 12.

Figure 12

Empirical and analytical covariances functions for (a) n = 300 and (b) n = 200.

Table 5 contains the parameters associated with the analytical covariance functions of the Tscherning and Rapp model (1974), that were estimated using the COVFIT program (Tscherning and Knudsen 2009) of the GRAVSOFT computational package.

The analytical covariance functions, in the application of the collocation method, were generated considering variations for the parameters: depth of the Bjerhammar sphere and scale factor. These parameters were modified as a function of the degree of adherence of the analytical covariance functions with the ECF.

2.11 Ecuadorian vertical offset estimation

According to IAG Resolution #16 (Tscherning 1983), “the indirect effect due to permanent yielding of the earth should be not removed.” Therefore, the zero-tide system is the most adequate tide system, whereby geometry is the mean/zero crust concept (Ihde et al. 2017).

Mäkinen (2017) highlighted the permanent tide system as one of the key points to be considered for the calculations involved in establishing the IHRS, and mentioned the following: “compute everything in ZT system, transfer to MT at the very end, using simplified formulas.” The computations carried out in the following were performed in line with this statement.

The ζ _res value at the EVD point was estimated using the LSFC method, carried out by using the FASTCOLC program of the GRAVSOFT package (Tscherning and Barzaghi 1991). The long and short wavelengths of height anomaly in the same point were then computed and added to the residual height anomaly value to get

The obtained height anomalies in the EVD point using two different d/o expansions of the GO_CONS_GCF_2_DIR_R5 model are listed in Table 6.

About the influence of heights from global models (for the calculation of gravity disturbances) on the accuracy of the values presented in Table 6, we must consider that, in case we use the EIGEN6C4 model, by computations on GNSS/lev points realized on the study area, we have a bias of 0.5 m, which implies a bias in gravity of 0.154 mGal.

As stated above, the values of Table 6 have been computed in the GRS80 and in the ZT system, accounting also for the difference between the GGM and the GRS80 GM values.

Finally, the W _ZT(P) potential values in the ZT system have been obtained as follows (Sánchez et al. 2021):

where γ(P) is the GRS80 normal gravity field value, U(P, GRS80) is the normal potential value, which has been evaluated using the closed formula (2-62) of Heiskanen and Mortiz (equation (25)).

where μ and β are the ellipsoidal harmonic coordinates (equations (26) and (27)), ω is the angular velocity of earth rotation, and E is the eccentricity associated with the considered reference ellipsoid.

where X, Y, and Z are the Cartesian geocentric coordinates for the calculation point.

The term

accounts for the difference between the new standard W ₀ value = 62636853.400 m² s⁻² (IAG resolution in Prague, 2015) and the GRS80 U ₀ value.

According to the IAG recommendations on IHRF, the W _ZT(P) values were then evaluated in the MT system as

where dH _EVD is the Ecuadorian VD offset in relation to W ₀ and in terms of normal height. The results of these computations are summarized in Table 7.