Geoid model determination for the Hellenic area “Hellas Geoid 2023”

Abstract

The latest geoid model "HELLAS GEOID 2023" (HG2023) derived by the Hellenic Military Geographical Service is the most comprehensive model for the entire Hellenic area. Long-term gravity data, orthometric and geometric heights, seabed topography and a high-resolution digital terrain model were implemented in the calculations. Data evaluation and their accuracy estimation were of major importance to ascertain compatibility among data sources. In this study, data from neighboring countries were used as they were essential for the completeness of the model developed. The technique remove–compute–restore was adopted for the separation of the gravity signal in the heterogeneous data utilizing EIGEN 6C4 (full degree and order 2190), since the global geodynamic model was found to fit best in the Greek region. Mean gravity residual contribution to the local geoid model was calculated using Stokes' theorem in the frequency domain and implementing the Fourier transform using the Wang and Core modification. The final surface of the resulting gravimetric geoid was adapted to the existing height system of the State utilizing the Least Squares Collocation method by fitting points of known orthometric and geometric heights, distributed throughout Greece. The external accuracy of the estimated geoid model was estimated at the level of 6 cm.

1 Introduction

The Geoid is essentially the real shape of the Earth, without topographic and atmospheric masses. In geodesy it is defined as the equipotential surface of the Earth’s gravity field which coincides with the sea surface in the absence of disturbing factors like tsunamis, ocean currents, salinities, wind, etc. and it extends through the continents (Vaníček and Krakowski 1986). Height determination and vertical control implementing a precise geoid model constitutes one of the most challenging research subjects of geodesy, and it attracts more attention since 1980s, related with the wide spread and the intensive use of GNSS techniques in surveying. Colorado experiment demonstrated that a 2 cm precision geoid model is achievable (Wang et al. 2021), using numerous geoid computations by different groups around the world at a rough area with a large topographic gradient.

The release of EGM2008 (Pavlis et al. 2012) defined a new era concerning the resolution and accuracy of geoid modelling because of its high degree and order due to numerous and heterogeneous data that contributed to computations. Since then, under the light of the new satellite dedicated gravity missions and terrestrial and shipborne gravity campaigns, other high resolution global models emerged such as EIGEN 6C4 (Förste et al. 2014), XGM2019e (Zingerle et al. 2020), SGG-UGM-2 (Wei et al. 2020). Besides their universal usage the global geopotential models constitute the base for local or regional models.

Greece has one of the longest coastlines in the world with a total coastal length of approximately 21,000 km and counts almost 30,000 islands, islets and rocks with 7,000 of them having area of more than 100 m² (HMGS, 2020). In this particular geomorphology, collecting the necessary density of high-precision ground and marine gravity data is a challenging task.

During the last decades, various studies (Balodimos, 1972; Arabelos et al. 1982; Tziavos, 1984; HMGS, 1983; Andritsanos et al. 1999; Andritsanos, 2000; Mylona, 2005; Vergos, 2006; Tziavos et al. 2009; Grigoriadis, 2009; Daras et al. 2010; Papadopoulos et al. 2019) conducted which are related to geoid model determination in the Greek area. This was mainly realized by adapting a GGM to the Greek classic levelling datum and altering the high frequency signal of it using local observations for a better fit. Exploitation of GPS data boosted geoid determination (Andritsanos et al. 1999), while investigation on the effects of topography and bathymetry resulted in defining Greek Geoid 2010 (Tziavos et al. 2009). Research in local geophysical dynamics led to the estimation of geopotential constants (e.g. W₀) to unify a local vertical datum connecting isolated height frames of islands to an equipotential surface defined on the mainland (Kotsakis et al. 2010).

In 2021, HMGS revised the State's first and second order gravity networks and created an updated database of gravity observations. Employing these data and a modern Digital Terrain Model (DTM) of high resolution (5 m), HMGS computed a new gravimetric geoid model for the Greek region. This surface is adjusted to points with known orthometric and geometric heights, and the resulting final surface is the first version of Hellenic Geoid 2022 v.1 (abbreviate as HG2022, Paraskevas et al. 2022). HG2023, which is described and discussed in this study, is an updated approach and constitutes the best fitted geoid surface to the orthometric elevations of the State’s triangulation network.

2 Data

2.1 Digital terrain model and seafloor topography

The DTM that has been used in this study is produced by the “Hellenic Cadastre” (a legal entity under the public law) by aerial photos of the Large Scale Orthophotos (LSO) project as a part of a wider project “Data and Information Technology Infrastructure for a contemporary cadastre” co-funded by the EU. Its nominal resolution is 5 m and the time period of the aerial photos was 2007–2009. HMGS completed and corrected the DTM using orthometric heights from the triangulation network. In addition, a higher resolution DTM (2 m) also derived from aerial photos of 2015–2016 by Hellenic Cadastre also, was used locally, where it was available, around several points in order to increase final precision. The overall combined DTM is validated at known elevation points and results are presented in Table 1.

Table 1 Statistics of the residuals derived by stations with known orthometric height (Fig. 2) and the produced 5-meter DEM (values are in meters)

The SRTM3arc_v4.1 (Farr et al. 2007) was used to supplement the aforementioned DΤM in the neighbor countries. The digital seafloor topography was retrieved from EMODnet Digital Bathymetry (2020), with a resolution of about 115 m. This model was checked and supplemented where necessary with data from other studies (Paraskevas et al. 2021; Mintourakis, 2020) and digitized contour lines from the Hellenic Navy Hydrographic Service (HNHS) maps, next to coast regions. Using the above datasets, a joint Digital Terrain Model and Seafloor Topography was created utilizing the coastline for the Greek area in vector format derived by HMGS. The joint model was downsampled to 10 m, 50 m and 100 m in order to be utilized in several steps in this study (Fig. 1).

Fig. 1

Produced joint digital elevation model with sea floor topography

¹

Fig. 2

Stations with known Orthometric Elevation. Black tringles represent triangulation points and red dots gravity measurements

2.2 Control points

In total, over 3000 benchmarks of Greece’s national triangulation network (more than 10% of the network), were used (Fig. 3). The majority of the above-mentioned benchmarks evenly distributed over the continental area of Greece were used as control points at various stages of the geoid computation process. It was assumed that the average accuracy of the orthometric heights was 5–10 cm (e.g. Ampatzidis et al. 2018) despite the fact that it was nominally reported (Takos, 1989) as 3 cm.

The geometric heights of the benchmarks were determined by GNSS observations from 2001 to 2022 in several campaigns (Anagnostou 2007; Papadopoulos et al. 2020). The estimated average accuracy of the geometric heights is about 2 cm.

Fig. 3

GNSS/Levelling Measurements. Green squares indicate GNSS observations over 12h, blue circles at least for 6h, while black triangles for at least 2h

2.3 Global geopotential model (GGM)

Geoid heights derived by GNSS/Levelling were compared with four global geopotential models developed to a full degree and order 2190 (Table 3), and it was concluded that the EIGEN6C4 model provided the best fit against the EGM2008, SGG-UGM-2 and XGM2019e models (Tables 3 and 4 and Fig. 4).

These models were checked by GNSS/Levelling observations at 3358 points and their geoid undulations. The determination of the GGM geoid undulations was carried out through the general formula (Rapp 1997) as in Eq. 1.

$$N=\zeta +\frac{{\varDelta g}^{FA}-0.1119H}{\stackrel{-}{\gamma }}{\rm H}+{{\rm N}}_{0}$$

(1)

where ζ and Δg^FA denote the height anomaly and free-air gravity anomaly, respectively, that are computed from the corresponding series expansions (up to n_max=2190) based on the spherical harmonic coefficients of each model and the GRS80 parameters. The value of 0.1119H is the attraction of an infinite Bouguer plate of height H assuming a constant density of 2.67 gcm^− 3 (Heiskanen and Moritz 1967), $\stackrel{-}{\gamma }$ is the mean normal gravity between each point on the ellipsoid and the corresponding point on the telluroid. The term N₀ represents the contribution of the zero-degree harmonic to the GGM-based geoid undulations and is computed using the Eq. 2 (Heiskanen and Moritz 1967):

$${N}_{0}=\frac{GM-{GM}_{0}}{R\gamma }-\frac{{W}_{0}-{U}_{0}}{\gamma }$$

(2)

where R and γ are regarded as mean values of the earth radius and normal gravity acceleration computed for the GRS80 ellipsoid (Moritz, 1992; Hofmann-Wellenhof & Moritz 2006, Table 2.5) and are summarized at Table 2.

Table 2 Adopted Initial values for Earth constants (Eq. 2).

By IERS conventions (Sanchez et al. 2015; 2016), Wo = 62636853.4 m²/s² and Earth's geocentric gravitational constant (GM) of the adopted GGMs is GM = 0.3986004415 × 10¹⁵ m³/s². Using the above values for the zero degree term in Eq. (2) yields N ₀= − 17.68 cm.

Table 3 GGMs used for the evaluation tests at the 3358 GNSS/levelling benchmarks

The computation of the GGM geoid undulations from Eq. (1) was performed in tide-free system with respect to GRS80 as a fixed reference ellipsoid. SGG-UGM-2 and XGM2019e were transformed to tide-free system using the conversion of the coefficient C₂₀ implemented in the Colorado Experiment (Wang et al. 2021):

$$C_{{20}}^{{TF}} = C_{{20}}^{{ZT}} + 3.11080x10^{{ - 8}} \times 0.3/\sqrt 5$$

(3)

Orthometric heights refer to the Hellenic Vertical Datum (HVD) so they were defined in the mean-tide system (${H}^{MT})$. Since ellipsoidal heights were defined in the tide-free (h^TF) system, the orthometric heights should be transformed to the tide-free one (${H}^{TF}$) following Ihde et al. (2008).

$${H}^{ZT}={H}^{MT}+0.0994-0.29541{\text{sin}}^{2}\varphi -0.00042{\text{sin}}^{4}\varphi$$

(4)

$${H}^{TF}={H}^{ZT}-0.06034+0.17901{\text{sin}}^{2}\varphi +0.00182{\text{sin}}^{4}\varphi$$

(5)

where φ is the geocentric latitude.

The analysis of the GGM geoid heights was carried out implementing the in-house Matlab software GRAVSynth (Papanikolaou and Papadopoulos, 2015), slightly modified in order to adopt the tide-free system by the transformation term given in Eq. 3 and the GRS80 Earth constants instead of the WGS84 ones. For each of the above four GGMs the geoidal undulation was computed at all points and afterwards compared to their measured values. Residual statistics and the systematic difference (bias) were exported as a report after Eq. 6.

$$\Delta N^{i} = {\text{ }}N^{i} _{{GNSS/LEV}} - {\text{ }}N^{i} _{{GGM}} - {\text{ }}bias$$

(6)

where Nⁱ_GNSS/LEV represents the GNNS-geoid undulation (N = h-H) of each point, Nⁱ_GGM is the corresponding calculated geoid undulation and bias is the systematic difference. The mean differences (bias) mirror in a way the average difference of the reference surface defined by each GGM compared to the Hellenic Vertical Datum. Statistics for every dataset are presented in Table 3, while statistics of the GGM comparison to the measured N are presented in Table 4 taking into account the bias.

Table 4 Comparison of the statistics of each GGM with measured geoid heights (GNSS/Levelling)

Fig. 4

Height-dependent variations of the differences N_GNSS/LEV - N_GGM at the 3358 GPS/leveling benchmarks: a N_GNSS/LEV - N_EGM08 b N_GNSS/LEV - N_{EIGEN 6C4} c N_GNSS/LEV - N_SGG−UGM, d N_GNSS/LEV - N_{XGM2019E_2159}.

The height variations in the GGM residuals presented in Fig. 4 over the test network reveal a linear dependency between orthometric heights and the geoid residuals (see e.g. Kotsakis et al. 2010). A slightly better performance is found at the EIGEN 6C4 model (Table 4 and Fig. 4). In Fig. 4, the EIGEN 6C4 model achieves a higher concentration of residuals near zero resulting in a lower declination of the associated line compared to the other models. Furthermore, the line’s constant factor exhibits a lower absolute value at low altitudes remaining closer to zero, indicating a weaker correlation with elevations than the other models.

2.4 Terrestrial data

HMGS’s archive consists of more than 35000 raw gravity observations which were collected since the 1950s. All gravity measurements were reduced to a single reference system and constitute a dataset which for the needs of this work was named "HMGS 2021" (Paraskevas et al. 2022). All the available datasets refer to the National Gravity Station located at the HMGS facilities according to the 2021 adjustment of the 1st and 2nd order gravimetric network. After a quality check procedure, many observations were rejected as outliers or repeated ones at the same station. The largest measured gravity variations over time at one station are up to three hundred µGal. The above observations were also checked for horizontal and vertical displacement. Horizontal positions were corrected using their initial description where it was necessary. Orthometric heights were compared to the 5m DTM, and large value residuals led to either a height correction using the 2m DTM or rejection. In addition, many of the old stations were re-surveyed using GNSS receivers in RTK mode (precision 3–4 cm). Gravity values of old stations accompanied with raw observations were recalculated, otherwise, the absolute gravity value of each station was reduced to HMGS 2021 datum. Moreover, positioning in gravity surveys after 2007 was obtained using GNSS receivers (in static or RTK mode).

Various procedures took place in order to detect localized gross errors in the gravity data including visual inspection and the application of the Least Squares Collocation method (Vergos et al. 2005). As an enhancement to the visual inspection an additional test was performed implementing a MATLAB script, to check each measurement by means of the following criterion:

$$For\,Distance\, \left(i, j\right)\,<2 \,Km$$

$$\text{And}\,H\,\left(i\right)-H\,\left(j\right)>{\text{1\,m}}$$

(7)

$$Gradient=\frac{\left|{gFA}_{i}-{gFA}_{j}\right|}{\left|H\left(i\right)-H\left(j\right)\right|}<1 \,mGal/m$$

where i, j every pair of points in a 2 km neighborhood, gFA are free air anomalies (calculated using Eq. 16), H is the Orthometric Elevation referred to HVD (Hellenic Vertical Datum). At the stations where the Gradient value (Eq. 7) was larger than 1 mGal/m, all measurements at their neighborhood were checked one by one for gross errors. The above criterion is based on empirical trials which are given in order to derive the most reliable results. The chosen threshold of 1mGal/m is an empirical value which returns a reasonable amount of points to check. Free air gravity anomalies are dependent on topography, and generally their spatial variation should be small. Using the above algorithm we applied a local inspection for gross errors in the dataset that could affect our results.

Finally, the differences between the observed and EIGEN 6C4 gravity values were calculated, and their statistics are listed in Table 5. After rejecting faulty values by the adopted criterion (Eq. 7), the remaining final absolute gravity values form the aforementioned database as GroupA, consisting of 21279 gravity points, distributed all over Greece (Fig. 5).

2.5 Marine data

In addition, iso-gravity contours and free-air marine gravity anomaly values from earlier works (Morelli et al. 1975a; Morelli et al. 1975b; Casten and Makris, 2001) are kept in HMGS’s archive. Marine gravity measurements were, also, retrieved from other sources (BGI, 2022). For the above data, raw observations were not available in order to be safely reduced to the same datum as the terrestrial data. Moreover, it was not possible to validate their accuracy for an optimal combination to create a new data set. For these reasons, more recent marine gravity measurements were sought out that would be used as check data. Datasets consist of marine gravity measurements (Paraskevas et al. 2021) called Group B and a grid of free air gravity anomalies (Mintourakis, 2014), called Group C. The Group B dataset was aligned to the chosen gravity reference system as there were available raw gravity measurements. In order to reduce Group C observations to the common datum, they were compared to gridded Free Air Anomalies of Group B and the calculated bias (1st order correction) was added to original data.

2.6 Neighbor countries

As far as it concerns the adjacent countries the only available gravity data were from open sources (BGI, 2022). The lack of dense raw observations and metadata led us to search for another source of data. For this study, it was decided to calculate Free Air anomalies using the selected GGM and the spectral enhancement method (Hirt et al. 2011). This method reduces the frequency gap between the global geopotential models and the external gravity field of the Earth. In the present study, we combine the EIGEN 6C4 model up to its maximum degree with components of higher frequencies. These components are actually the terrain effects of the so-called residual terrain model (RTM) (Forsberg, 1984). This method has been successfully applied in many studies (Forsberg, 1985; Hirt et al. 2010; Sampietro et al. 2017; Zaki et al. 2018) to improve the performance of a GGM. The effects of the Residual Terrain Model (RTM) correspond to the topographic effects with respect to a reference surface. These effects can be linked to a Bouguer plate taken between the real topographic surface and the reference surface from which classical terrain corrections are subtracted (Forsberg, 1984):

$${\varDelta g}_{RTM}\approx 2\pi G\rho \left({\rm H}-{{\rm H}}_{ref}\right) - TC$$

(8)

where H and ${{\rm H}}_{ref}$ are the heights of the detailed DTM (or the station’s height where available) and the reference surface, and TC corresponds to the classical terrain correction, calculated through Eqs. 17, 18 and 19. The RTM approach requires a smooth, long-wavelength reference surface to be subtracted from the detailed elevation model. If the RTM is used as an augmentation of GGMs beyond their maximum resolution, subtraction of a spherical harmonic reference surface of the same degree is most suitable to high-pass filter the detail elevations (Hirt, 2010; 2013).

ETOPO1 (Ice Surface) (Amante & Eakins, 2009; Pavlis et al. 2012), was also used as a reference DTM in Eq. 9, at the same resolution (0.09 degree) as the one used for the evaluation of EIGEN 6C4. A grid of gravity anomalies with resolution of 0.0016 deg (~ 1600 m) derived from the GGM, was corrected to RTM and compared to the produced Free Air gravity anomaly grid from the terrestrial measurements. The comparison of these two grids showed a shift (1st order correction) between them which was added to that derived from the GGM. Finally, sampling this grid produced point data at its nodes which were included in the final Database as group D.

$$\Delta g^{{Neighbors}} = \Delta g_{{EIGEN6C4}} |_{2}^{{2190}} + \Delta g_{{RTM}} |_{{ETOPO1 - 9000\;m}}^{{SRTMv4.1 - 90\;m}} |\underrightarrow {Static{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} Adjustment{\text{ }}to{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} Adopted{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} DATUM}\Delta g_{{final}}^{{Neighbors}}$$

(9)

2.7 Database

All the aforementioned data merged into an overall database consisted of four separate groups. Group A counts 21279 points of terrestrial measurements observed from 1966 to 2022. Group B comprises 2450 points of marine gravity data from other studies (Paraskevas et al. 2019, 2021). Also, 142869 nodes of a Marine Gravity Data grid with a step of 0.01 degrees forms Group C. Finally, 50373 points of gravity data from neighboring countries in a grid with step of 0.016 degrees, accommodate in Group D. Absolute gravity values of the database points were aligned to the 2021 National Gravity Station value. The database was expressed in geographic coordinates with respect to ITRF2008 (Altamimi et al. 2011)). The Orthometric heights of the database points (if they are available) refer to the HVD which has its starting point (datum) the tide gauge of Piraeus, for the continental Greece, Heraklio for Crete and local tide gauges, at the main islands (Kos, Rhodes, Karpathos etc). Finally, an atmospheric correction was applied to all gravity anomaly values using the formula (Wenzel, 1985):

$$dgatm=0.874-9.9*{10}^{-5}h+3.56*{10}^{-9}{h}^{2},$$

(10)

where dg_atm is in mGal and h is in meters.

Table 5 Statistics of free air gravity anomalies for each group and comparisons to EIGEN 6C4 (D_{X−EIGEN 6C4})

As it can be seen in Table 5, after subtracting the GGM, the terrestrial gravity data

(D_{A−EIGEN 6C4}, D_{D− EIGEN 6C4}) exhibit a larger standard deviation value than the marine ones (D_{B−EIGEN 6C4}, D_{C− EIGEN 6C4}). The standard deviation of the marine gravity data is smaller than the one of the terrestrial data since the gravity field at sea is smoother than the one on land (practically, there is no terrain effect).

2.8 Reduction scheme

Over continental regions, the geoid is mostly located inside the topographic masses, whereas gravity measurements (and levelling to compute gravity anomalies) are taken on the topographic surface. To satisfy the Stokes’ Boundary Value Problem (BVP), the surface anomalies need to be continued downward to the geoid level. In this research Helmert’s second condensation method (Helmert, 1884; Hofmann-Wellenhof and Moritz, 2006) is chosen, which is the method mostly used in practice. “Helmert's second method of condensation consists in condensing all matter external to the geoid into a surface lying in the actual geoid itself” (Lambert, 1930).

According to Martinec et al. (1993), the Stokes-Helmert approach (i.e. using the Helmert’s second condensation method for condensing the topography and Stokes’ integral for obtaining the geoid undulations) could be summarized as follows:

• 1. Computation of the direct topography effect on the observed gravity by replacing the effect of topographical masses on gravity at the Earth’s surface with the effect of the mass layer on the geoid (e.g. the Bouguer plate reduction and the classical terrain correction);

• 2. Calculation of the downward continuation to obtain gravity anomalies (or the vertical gradient of gravity) and reduction of the observed gravity anomalies from the Earth’s surface to the geoid (e.g. the free-air reduction);

• 3. Determination of the co-geoid by computing the reduced gravity anomalies with Stokes’ integral; and

• 4. Computation of the indirect topographical effect on the gravity potential in order to define the geoid.

According to Sideris and Forsberg (1991, Eq. 17), modified by Martinec et al (1993, Eq. 1), steps 1 and 2 could be summarized in Eq. 11:

$${{\varDelta g}^{H}}_{red}={g}_{obs}-\gamma +FA-{A}_{p}+ {A}_{p}^{comp}$$

(11)

where Δg^H_red, is the reduced gravity anomaly at the geoid, g_obs is the observed gravity at the surface of Earth, γ is normal gravity at the ellipsoid, FA is Free Air reduction to the geoid, A_p is the attraction of the topographical masses above the geoid and ${A}_{p}^{comp}$ is the attraction of the condensed topography on the geoid. Sideris and Forsberg (1991, Eq. 22), also show that

$$TC=-{A}_{p}+ {A}_{p}^{comp}$$

(12)

where TC is the classic terrain correction. In other words, mean terrain corrected free air gravity anomalies, termed mean Faye gravity anomalies (Heiskanen and Moritz, 1967) and were used to approximate mean Helmert gravity anomalies at the geoid as in many other studies (e.g. Featherstone et al. 2011, Matsuo and Kuroishi 2020, Claessens and Filmer 2020). Finally Δg_red, used in the Stokes’ integral, are summarized in Eq. 13 (Martinec et al. 1993):

$${\varDelta g}_{red}={g}_{obs}-\gamma +FA+TC+{\delta g}_{ind}$$

(13)

where δg_ind is the secondary indirect effect on gravity. Normal gravity (γ) calculated using the Somigliana’s formula (Somigliana, 1930) as in Eq. 14, and the other quantities are as follows:

$$\gamma =\gamma e*\frac{1+k{sin}^{2}\varphi }{\sqrt{1-{e}^{2}{sin}^{2}\varphi }}$$

(14)

where γe = 9.78032677153489285793 ms^− 2, k = 0.0019318513532606763607 and e² = 0.00669438002290341574957 (Hinze, 2013).

Fig. 5

Gravity data used in this study. Red dots indicate Group A, blue dots Group B, light blue dots Group C and brown dots Group D. Green square are the first and second order national Gravity Network revised in 2021.

Free air reduction:

Free air reduction was calculated using simple free air gradient following Eq. 15 (Heiskanen & Moritz, 1967, 131):

$$FA=-\frac{\partial g}{\partial H}H\approx -\frac{\partial \gamma }{\partial H}H\approx -0.3086H$$

(15)

where H is the orthometric height in meters and thus Free Air gravity Anomalies (gFA) are calculated from the observed absolute gravity values (gobs) using Eq. 16:

$$gFA=gobs-\gamma -FA ,$$

(16)

Terrain correction (TC):

Terrain corrections were calculated using the algorithm proposed by Nagy (1966) and Kane (1962). Near topography and the produced 10m DTM were implemented in TC calculation. Three zones around each point contributed to TC separated by their radial distance. In the near zone, 0 to 10 m from the station, the algorithm sums the effects of four sloping triangular sections, surveyed around the station. In the case where there weren’t field measurements, DTM elevation values at each diagonal corner of the pixel were used. TC₁ is computed by Eq. 17 (Nagy, 1966):

$${TC}_{1}=G\rho \varphi \left(R-\sqrt{{R}^{2}+{H}^{2}}+\frac{{H}^{2}}{\sqrt{{R}^{2}+{H}^{2}}} \right),$$

(17)

where G is the universal constant of gravitational attraction, ρ the terrain density (here 2.67g/cm³), φ the horizontal angle of the triangular section, H the difference between the station elevation and the average elevation of the diagonal corner, and R the specific distance (10 m).

In the intermediate zone (10 m–100 Km), the TC (or terrain effect) is calculated for each point using the flat-topped square prism approach of Nagy (1966). The terrain height is measured on a regular grid (equal spacing) with a total of N×M points extracted from the DTM in use (here 10 m resolution DTM). Then the terrain correction can be computed by summing up these N×M prisms using Eq. 18:

$${TC}_{2}=G\rho \left|\genfrac{}{}{0pt}{}{{z}_{2}}{{z}_{1}}\right.\left|\genfrac{}{}{0pt}{}{{y}_{2}}{{y}_{1}}\right.\left|\genfrac{}{}{0pt}{}{{x}_{2}}{{x}_{1}}\right.xln\left(y+r\right)+yln\left(x+r\right)-z{sin}^{-1}\left(\frac{{z}^{2}+{y}^{2}+yr}{\left(y+r\right)\sqrt{{y}^{2}+{z}^{2}}}\right)\left|\right||$$

(18)

where r is the distance between a unit mass and the station.

In the far zone (100 km–166.7 km), the terrain effect is estimated by the approximation of an annular ring segment to a square prism, as described by Kane (1962). An 100 m DTM was implemented which was produced by downsampling the 10 m one. The gravitational attraction in this zone calculated by Eq. 19 (Kane, 1962):

$${TC}_{3}=2G\rho {{\rm A}}^{2}\frac{{R}_{2}-{R}_{1}\sqrt{{{R}_{1}}^{2}+{H}^{2}}-\sqrt{{{R}_{2}}^{2}+{H}^{2}}}{({{R}_{2}}^{2}-{{R}_{1}}^{2})}$$

(19)

where A is the length of the horizontal side of the prism, R₁ the radius of the inner circle of the annular ring, R₂ the radius of the outer circle of the annular ring and H the height of the annular ring or prism. For marine gravity data TC was calculated only with Eq. 19 using the produced 100 m resolution DTM. Final results were tested in two small onshore areas using TC GRAVSOFT software (Forsberg and Tscherning, 2008) and they agree at the level of 0.4 mgal.

Indirect effect

The indirect effect on elevations due to Helmert mass compensation calculated from the joint DTM with seafloor topography. This value is the sum of two terms in the planar approximation (Wichienharoen, 1982):

$${N}_{IND}=-\frac{\pi G\rho }{\gamma }{{\rm H}}_{\left({x}_{p},{y}_{p}\right)}^{2}-\frac{G\rho }{6\gamma }{\iint }_{A}\frac{{{\rm H}}_{\left(x,y\right)}^{3}-{{\rm H}}_{\left({x}_{p},{y}_{p}\right)}^{3}}{{\left[{\left({x}_{p}-x\right)}^{2}+{\left({y}_{p}-y\right)}^{2}\right]}^{\frac{3}{2}}}dxdy ,$$

(20)

where x_p, y_p are the coordinates of the point for which the indirect effect is calculated, x, y are the coordinates of the current point of the integration area A during the calculations, γ the normal gravity on the reference ellipsoid, A is the integration surface with the elevation information and H is the orthometric height.

For the indirect effect of topography calculation, the first term of Eq. 20 was computed point by point from the available elevations, while the second term was calculated using the DTM with a 50 m resolution at integration limits of 1°.

The indirect effect on gravity for this reduction scheme can be given as (Sideris and She 1995; Kuroishi 1995):

$${\delta g}_{ind}\approx \frac{2\pi G\rho {H}^{2}}{R}\approx 0.3086 {N}_{IND}$$

(21)

Statistics of the above processing steps are listed in Table 6.

3 Calculations

3.1 Geoid determination

The selected method for geoid determination in this research is the 2D-FFT approach implementing the "Remove-Compute-Restore" technique (Schwarz et al. 1990). According to this method, the geoid model is computed by aggregating three terms at the restore stage:

$$N = N_{{EIGEN{\text{ }}6C4}} + N_{{IND}} + N_{{res}}$$

(22)

where N_{EIGEN 6C4} is the long to medium-wavelength contribution (calculated from the Global Geodynamic Model in use) as described in Sect. 2.3. N_IND is the indirect effect on elevations due to Helmert mass compensation (calculated from the joint DTΜ with seafloor topography from Eq. 20). The third term in Eq. 22 (N_res) is the contribution of the mean residual gravity anomaly to the geoid, which is calculated by the Stokes’ integral (Heiskanen andMoritz 1967):

$${N}_{res}=\frac{R}{4\pi \gamma }{\iint }_{\sigma }{\varDelta g}_{res}\left(\psi ,\alpha \right)S\left(\psi \right)d\sigma ,$$

(23)

where γ is the normal gravity on the surface of the ellipsoid, R is the mean radius of the earth, σ is the integrating surface and S(ψ) is the Stokes’ function calculated in the form of Legendre polynomials P_n as (Hofmann-Wellenhof and Moritz, 2006):

$$S\left(\psi \right)=\sum _{n=2}^{\infty }\frac{2n+1}{n-1}{P}_{n}\left(cos\psi \right),$$

(24)

Δg_res is the residual gravity anomaly calculated by the Helmert ‘s second method compensation of masses. It is given by the relation (Heiskanen & Moritz 1967; Sideris 1990):

$${\varDelta g}_{res}={\varDelta g}_{red}-{\varDelta g}_{GGM}$$

(25)

The gravity anomalies should be located on a regular grid for the computation of the geoid by means of the FFT technique. Since the gravity data set consists of point data anomalies distributed randomly, they need to be interpolated in order to obtain a regular data grid. Thus, an accurate interpolation procedure should provide a regular grid data with a high reliability level. Helmert’s second method of condensation is not a gravity smoothing terrain reduction and thus it is recommended to use other reduction methods in order to obtain smooth gravity anomalies suitable for gridding (Bajracharya et al. 2001). For this reason, before the interpolation, the short-wavelength effects that appear in the gravity anomaly field, associated with the short-wavelength topography and bathymetry, should be removed. This correction is called the RTM correction (Forsberg and Tscherning, 1997) and should be applied to get smoother gravity anomaly field for gridding. The spatial resolution of the smoothed reference surface must be equivalent to that of the GGM in use to represent the low frequencies of the gravity field. In order to obtain the smoothest gravity anomaly, data were gridded after the remove step of the geoid determination procedure using the following formula:

$${ \varDelta g}_{res\_Final}^{GRID}={\varDelta g}_{res}^{GRID}+2\pi G\rho {\left({\rm H}-{{\rm H}}_{refEIGEN6C4}\right)}^{Grid}$$

$${\varDelta g}_{res}^{GRID}{\underleftarrow{interpolation}\varDelta g}_{red}^{Points}-2\pi G\rho {\left( {\rm H}-{{\rm H}}_{refEIGEN6C4}\right)}^{Points} -{\varDelta g}_{EIGEN6C4}^{Points}$$

(26)

where Δg_red is the reduced free-air gravity anomaly computed using Eq. 13, ρ is the density of the topography (2.67 g/cm³) on land or the density of the topography minus sea water density (2.67–1.027 = 1.643 g/cm³) for marine RTM, Η is the elevation, Η_refEIGEN6C4 denotes the elevation of the reference surface with the same resolution as the GGM EIGEN 6C4 (5’ or about 9.1 km), TC is the terrain correction computed for land and marine points and Δg_EIGEN6C4 is the gravity anomaly computed from the geopotential model EIGEN 6C4. In the final step we restored the term $2\pi G\rho \left({\rm H}-{{\rm H}}_{refEIGEN6C4}\right)$, in order to calculate the geoid (otherwise, quasi geoid will be determined). This final grid could be used in Eq. 23 to calculate N_res (Figs. 6 and 7, Table 7)

Table 6 Statistics of the processing steps in reducing gravity anomalies to Database (values in mGal)

Table 7 Statistics of produced grids

According to Table 6 it is obvious that reducing the observed gravity to the geoid using Eq. 13 increased the standard deviation, which means that reduced values (Δg_red) are not as smooth as the free air ones (gFA). After subtracting GGM and RTM reduction using Eq. 26Total std is limited to 8 mGal resulting in a smoother surface for interpolation. The produced Δg_res grid appeared to be smooth enough and its standard deviation is 9.2 mGal (Table 7), while after restoring of RTM reduction the final grid is rougher (std ~ 13.7 mGal) containing all the available information for the Stokes’ integral.

To mitigate the truncation error, due to non-integration over the whole earth, the Stokes’ function should be modified. In the context of the present study, the modification proposed by Wong and Gore (1969) was chosen. Wong and Gore modified the Stokes’ kernel by removing the low-degree terms of the Legendre polynomials from the original kernel following the Eq. 27 (Wong and Gore1969):

$${S}_{mod}\left(\psi \right)=S\left(\psi \right)-\sum _{n=2}^{N}\frac{2n+1}{n-1}{P}_{n}cos\psi ,$$

(27)

Modification parameters are the integration’s limit (i.e. the radius of the circular surface σ in the integral of Eq. 23) ψ₀ and the maximum degree of the polynomial N in Eq. 27. To determine the appropriate modification parameters, we employed geoid heights, obtained by GNSS/leveling as reference. The modification parameters were selected by conducting tests for the optimum values that minimize the standard deviation of the difference between the computed and surveyed (by GNSS/leveling) geoid heights at 3258 benchmarks. Table 8 summarizes the statistics (standard deviation) of the above differences from various sets of the two modification parameters. The results show that modification parameters values of ψ₀ = 3° and N = 1400 suit best with the GNSS/leveling geoid heights.

Table 8 Standard deviation of the differences between the computed for various values of modification parameters (ψ₀ and Ν) and GNSS/leveling geoid heights at 3258 benchmarks (values in meters)

The final geoid surface is computed by adding the three grids:

$${N}_{grav}^{Grid}={N}_{EIGEN 6C4}^{Grid}+{N}_{Ind}^{Grid}+{N}_{Res}^{Grid}$$

(28)

Statistics of various calculation steps using 3258 points where GNSS/Levelling values are available and of the calculated grids are listed in Tables 9 and 10, respectively and, the final produced gravimetric geoid surface is presented in Fig. 8.

Fig. 6

Residual gravity after interpolation.

Fig. 7

Final residual gravity anomaly grid after restoring the quantity $2\pi G\rho {\left({\rm H}-{{\rm H}}_{refEIGEN6S}\right)}^{Grid}$ using Eq. 26.

Table 9 Statistics of processing steps in gravimetric geoid determination, calculated at surveyed points by GNSS/Levelling

Table 10 Statistics of grids produced during the processing of geoid determination

3.2 Alignment to the Hellenic vertical datum

The alignment of a local gravimetric geoid to a local vertical datum is realized by Saadon et al. (2022). The main reason of this kind of alignment is the present status of the Hellenic Vertical Datum.

The resulting gravimetric geoid (Fig. 8) is a surface referenced to the datum of the GGM in use as it expresses the low frequencies in the procedure and represents the majority of the signal in the computations. On the contrary, the reference surface of the orthometric heights (vertical datum) in Greece (Hellenic Vertical Datum – HVD) has been determined by the Piraeus tide gauge for the continental area, by the Heraklion tide gauge for Crete and by local tide gauges (Lesvos, Kos, Leros, Megisti etc.) for the island part of Greece. This complexity requires a more sophisticated and comprehensive procedure than a four - parameter transformation which is a pure deterministic approach. In the context of this specific study, Least Squares Collocation was implemented in order to fit the computed gravimetric geoid surface to the local datum by minimizing their difference as in Eq. 29.

$$\varepsilon {\text{ }} = {\text{ }}N^{{GNSS}} {-}{\text{ }}N^{{gravimetric}} {-}{\text{ }}\mu$$

(29)

where N^GNSS is the geoid undulation obtained by GNSS/Levelling (h-H), N^gravimetric is the geoid undulation of the gravimetric geoid as calculated in this study and µ the constant term, absorbing the systematic inconsistencies among the numerous reference surfaces (continental Greece and islands). The covariance relation applied in Eq. 30, fulfills the conditions of its application (positive defined, depends only on distance) in point fitting (Forsberg and Tcherning 2008):

$$C_{{eh}} = {\text{ }}C_{0} \left( {1 + \alpha {\text{ }}s} \right){\text{ }}e^{{ - \alpha {\text{ }}s}}$$

(30)

where C₀ is the variance, s is the distance and α is the parameter that determines the correlation magnitude, while the remaining quantities are determined from the data with a least- squares fit using GRAVSOFT software (Forsberg and Tscherning 2008). The choice of correlation distance depends on the distribution and density of the measured geoid undulations. In this particular study, a correlation distance of 20 km was chosen in order, on the one hand, to respond to the density of the available points and, on the other hand, to eliminate as possible, the correlations among different reference surfaces (due to different vertical datum realizations). The covariance depends only on the relative position of the points and not on their exact coordinates (homogeneity) and on their distance and not on the direction (isotropy) (e.g. Dermanis, 1987).

One may raise doubts about the prevention of the homogeneity and isotropy of the covariance function, due to presence of different biases among different datasets and sources. As a first comment, in Eq. 29, the constant term µ absorbs part of the systematics, at least as a mean average for all Greece’s consisting vertical systems. In addition, as it is mentioned before, the relatively small correlation length mitigates –as far as possible- the interference of points which are located far away from each other. Finally, we may again underline that due to the mosaic of different reference surfaces in the Hellenic Area we should find a practical way of amalgamating the different vertical datums in terms of a unified geoid computation (Fig. 9).

The final surface of the geoid is calculated using Eq. 31:

$$N_{{HG2023}} = {\text{ }}N_{{gravimetric}} + {\text{ }}\varepsilon _{{grid}}$$

(31)

In total, 3000 points were used for the implementation of Eqs. (30) and (31) in order to align the gravimetric geoid to the Hellenic Vertical Datum (Fig. 10). A crucial process regarding the reliability of the model is its evaluation through cross validation points. As control points 258 stations were chosen, distributed on the mainland and on the islands (Fig. 10). At the control points, we apply the model estimated from the 3000 initial set of points. Table 11 refers to the results of the alignment of the gravimetric geoid to the Hellenic Vertical Datum (upper section) and the external validation of the solution at the 258 points (cross validation points, lower section). The RMS of the differences at the cross-validation points is 6.2 cm which can be considered satisfactory, taking into account the nature of the measurements (old and relatively sparse) and the existence of many islands especially, since -as it has already mentioned- the Hellenic Vertical Datum is not unified and each island realizes its own vertical datum.

Table 11 Statistics of GNSS/LEVELLING points used in LSC and cross - validation points

Fig. 8

Gravimetric Geoid. Contour interval is 1 m

Fig. 9

Correcting Surface according to LSC (see Sect. 3.2, ε_grid in Eq. 31)

Fig. 10

Distribution of 258 control points presented as green squares and 3000 points used in LSC as black triangles

Fig. 11

Hellas Geoid 2023 v1. Contour interval is 1 m

Table 12 presents the statistics of the differences between the geoid undulations derived from 258 GNSS/Levelling benchmarks used as cross-validation points and the selected GGM for Greece, the calculated Gravimetric Geoid Model (Fig. 8) and HG2023 (Fig. 11). Unfortunately, there are no other publicly available local geoid models for Greece in order to compare them to HG2023. The standard deviation of the geoid difference between EIGEN 6C4 and GNSS/Leveling control points is estimated 15.81 cm. The gravimetric geoid model performs better by 10% or 1.6 cm, while HG2023 adjusts better to cross-validation points by 65% or 10.35 cm than the EIGEN 6C4 model resulting in standard deviation value of about 5.5 cm.

Table 12 Statistics of 258 GNSS/LEVELLING points used as cross - validation points

4 Conclusions

Combining heterogenous data in local geoid computations is a challenging issue as they must be transformed to the same gravity datum, consistent in spatial resolution and refer to a common epoch. Nowadays, there are GGMs providing a reliable base for a regional or/and local geoid model. Dense local gravity data contribute the high frequency part of the signal into the model and provide high accuracy of the local geoid. In positioning we consider the dependence of orthometric heights on gravity and elevation datum is a necessity for a complete spatial reference frame.

Modern positioning techniques focus on point survey easily by GNSS observations and a geopotential datum on which physical values can depend (e.g. orthometric heights or/and gravity). Both of them combined together with an appropriate projection, resulting in fulfilling the basic geodetic infrastructure for a State or a wider region. HG2023 stands for the geopotential part of Greece’s new coordinate system towards the adoption of a precise and modern datum. Time variation of the datum is a challenging issue that will be the subject of future work. Greece is an active tectonic region, movements and deformation occur as trends or as abrupt discontinuities and time-depending parameters in datum definitions are of major importance.

This work is the culmination of a long-term effort and systematic GNSS and gravity observations throughout Greece from 2001 till now in the direction of creating a modern coordinate system covering the needs of the Hellenic Army, as well as the country's technical and scientific community. The resulting geoid model is the connection of the National Height System of Greece (HVD) with geometric heights, in the context of performing modern measurements using satellite receivers. The gravimetric surface computations implemented data from the HMGS long-term archive, containing surveys from over 60 years. The most recent measurements led finally to the validation and update of the overall gravimetric data base. The determination of the geometric height through satellite observations (GNSS) contributed to adapting the gravity geoid to the HVD. The resulting final surface is the optimally adapted geoid model to the existing orthometric heights. The external evaluation of HG2023’s accuracy is found at the RMS level of 6.2 cm which could be considered satisfactory, taking into account all the aforementioned problems of the geoid estimation in Greece. HG2023 represents the optimal equipotential surface with a grid step of 0.01° or about 1000 meters, which the HMGS publicly releases (www.hmgs.gr) for use, validation and additional comments. It is acknowledged that theoretically more rigorous approaches do exist for computation and downward continuation of mean Helmert gravity anomalies (e.g. Vaníček and Martinec 1994; Vaníček et al. 1999; Vaníček et al. 2013), while other methods are used worldwide for geoid determination. A future plan will be to implement different methodologies for geoid computations (like LSC or 1-D FFT), or the method handling the topography. However, the gravity dataset should be enriched in regions with poor distribution. Finally, it is also necessary to establish a new National Gravity Station (using absolute gravimeters) and connect it in a global gravity network in order to be more consistent to the global models.

A new version of the local vertical datum comprised of revised orthometric values surveyed by modern means and tied to a single reference is essential. In this case an upgraded geoid model could reach a higher accuracy.