Regional Gravity Model of Greece Based on Satellite, Marine and Terrestrial Data

Greece is a country structured by land, several islands and sea. A regional gravity model of such a country demands the involvement of several types of gravity data in order to cover all its territory. In this paper, we present the development of a regional combined gravity model of Greece and its surrounding area, by integrating terrestrial, marine (shipborne and altimetry-derived gravity data) and satellite data (GOCE and GRACE data). These kinds of data, especially the terrestrial and marine ones, have passed quality and validation control since they were collected from different organizations, which means that they probably have been acquired and processed with different parameters and formulas. Following that, their integration was accomplished with the application of the least-squares collocation (LSC). Therefore, a newly combined regional gravity model of Greece and its surrounding areas has been developed. This gravity model can be valuable for regional geological and geophysical studies of Greece, since it provides homogeneously the distribution of Complete Bouguer anomaly (CBA) all over Greece, with updated gravity data and in good agreement with the initial gravity datasets that have been used.


Introduction
Greece's territory presented study area are characterized from the existence of the mainland but also the sea area (Ionian and Aegean Sea), where there are also numerous islands. Greece is located in the eastern Mediterranean Sea, a country of south-eastern Europe, the southernmost of the Balkan countries. When it comes to gravity coverage of this territory, everyone can understand that the terrestrial gravity measurements are not adequate for the whole study area. Even the terrestrial measurements on the islands of Greece are not enough, since the Ionian and Aegean Seas cover a big part of the investigation area. Therefore, marine measurements also have to be taken into consideration. But unfortunately, we may still have areas with data gaps. These areas could be at the land, due to steep terrain or unapproachable areas, or even sea areas with no shipborne gravity data. Therefore, we can use the satellite gravity data existing in order to fill these data gaps. These types of data, especially generated from the GRACE/GOCEbased satellite-only global geopotential models (GGMs), are available for free and are widely used for geological purposes (Alvarez et al., 2015;Braitenberg, 2015;Choi et al., 2006;Eshagh & Pitoňák, 2019;Fuchs et al., 2016;Pal & Majumdar, 2015;Vaish & Pal, 2015). The point is that all these different inhomogeneous types of data have to be integrated in one common regional gravity model, valuable for regional geological and geophysical applications and modelling (Belay et al., 2021;Sadiq et al., 2010;Sobh et al., 2019). Mo re specifically, in this paper, by Regional gravity model, we refer to a combined Complete Bouguer Anomaly (CBA) map of Greece.
Several authors have published so far gravity maps of Greece, mostly based on terrestrial and marine data. Some of them acquired the majority of the data, especially the terrestrial ones while others compiled together smaller and bigger existing gravity databases. Below we mention the most important ones. Lagios et al., (1988Lagios et al., ( , 1996 created the CBA of Greece based on the database of the Section of Geophysics & Geothermics of the University of Athens, including terrestrial 33,000 points. Makris et al. (2013) compiled the gravity map of Greece from which he produced a 3-D gravity model constrained also by seismic data. Grigoriadis (2009) and Grigoriadis et al. (2016) focused primarily on the study of the Earth's gravity field in the Hellenic region from a geodetic and geophysical perspective. They developed a new gravity database of high accuracy and resolution combining heterogeneous data, such as geopotential models (EGM96 and EIGEN-GL04C), terrestrial/marine gravity (several existing databases such as Lagios et al., 1996), satellite altimetry-derived gravity data (DNSC08) and height data. Most of the previous gravity maps include the ''Morelli'' marine gravity data (Behrend et al., 1996;Morelli et al., 1975a, b), covering not only the Greek seas but also most of the Mediterranean Sea. Finally, Natsiopoulos et al. (2023) present their improved geoid model for the Hellenic area, using GOCE SGG data.
For the purposes of the terrain corrections required to calculate the CBA of the gravity data, we proceeded to the production of a Digital Elevation Model (DEM) of Greece, with cell grid 0.005°. It has been developed based on the combination of the European DEM (EU-DEM), v.1.1 (Copernicus, 2021) with Gebco bathymetry data (GEBCO Compilation Group, 2021). The maximum depth at the sea is 5123 m while the maximum elevation on the land is 3046 m.

Data
As we have already mentioned, a regional gravity model is constructed based on different types of gravity measurements that need to be homogenized. In this paper, we will integrate terrestrial, marine and satellite data for the production of the combined regional gravity model of Greece.

Terrestrial Data
Regarding the terrestrial gravity data there are several datasets as mentioned before. For that reason, we checked and considered all the available options of existing datasets and Bouguer maps that cover all the Greek area (Grigoriadis, 2009;Grigoriadis et al., 2016;Lagios et al., 1988Lagios et al., , 1996Makris et al., 2013). After checking and considering all the alternative existing gravity datasets, we ended up using the most updated and validated terrestrial gravity dataset, the one developed by the Hellenic Military Geographic Service. Its data are provided as a Free Air Anomaly map (Hellenic Military Geographical Service, 2021) with grid spacing 0.008°9 0.008°, a development based on gravity data originated from diachronic gravity measurements all over Greece. According to the technical manual of the map (Kagiadakis et al., 2021), the initial amount of gravity points collected from 1950s to nowadays was 26.832, including though several duplicate points and stations with reduced horizontal gravimetric or elevation accuracy. The quality control procedure followed by HMGS (Kagiadakis et al., 2021) includes initially the review of their gravimetric network of I and II class. Afterwards, the processing of the available data has been carried out in order to homogenize them and calculate their accuracy. The loops were processed from the beginning, based on the prementioned review of the network in order to absorb the geodynamic alterations of the Earth's crust due to the diachronic changes among the measurements. Finally, only the gravity points with standard deviation less than 0.5 mGal were accepted, based on Kagiadakis et al. (2021). For that reason, the Free Air Anomaly map and therefore the dataset used in this paper contain 18,358 points, distributed all over Greece. The used coordinate system is WGS84, while the Somigliana (1930) normal gravity (g o ) equation has been used, with the parameters for WGS84 reference ellipsoid, as being described in Blakely (1995) (Eq. 1). The latitude corrections were carried out using the Gravity Formula of 1984 (Eq. 2).
The Simple Bouguer Anomaly (g SBA ) has been initially calculated (Eq. 3). The assumed constant density for the Bouguer correction has been set up to 2.67 g/cm 3 , generally used by several other researchers for the Greek area (Dilalos & Alexopoulos, 2017Dilalos, 2018). The coordinate system is the Hellenic Geodetic Reference System (EGSA'87). Beyond that, the terrain corrections have been calculated using the 'Gravity and Terrain Correction' extension of Oasis Montaj (Geosoft, 2010), based on the DEM depicted in Fig. 1, for a radius up to 167 km. The terrain corrections with Oasis Montaj extension are calculated using a combination of the method of flat-topped square prism by Nagy (1966) and the method square prism by Kane (1962). Beyond that, three more gravity datasets, acquired from the authors during the previous years, have been embedded to the general (HMGS) dataset. More specifically, a total of 1478 measurements from Athens basin (Dilalos, 2018;Dilalos & Alexopoulos, 2017, 2019a, b, 2020Dilalos et al., 2019a, b), Pineios river delta of Thessaly (Dilalos et al., 2022) and southeastern basin of Thessaly (unpublished results) have been added. Therefore, the final terrestrial gravity dataset used in this paper contains 19,836 gravity points, covering all the mainland and islands of Greece. Finally, the CBA (g CBA ) has been calculated (Eq. 4), including the calculated Terrain Corrections.
In the above equations, g FA is the Free-Air Gravity values of the prementioned HMGS dataset in mGal, q is the selected Bouguer density of the rock slab in g/cm 3 , g curv is the Curvature correction applied, based on LaFehr's formula (LaFehr, 1991) and), g TC is the calculated Terrain Correction in mGal.
The Terrestrial CBA Map is illustrated in Fig. 2, with values ranging from -251.00 mGal to 156.80 mGal. The lower values are located mostly across the mountainous area of Pindos ridge (western Greece) but also at the area of eastern Macedonia. On the other hand, the higher values are located on the islands of Aegean along with their surrounding mainland (Attiki-Evoia) but also part of Crete.

Marine Data
For the gravity information of the Greek sea areas two types of gravity data have been used: shipborne available data from several previous researches and altimetry-derived gravity data, mostly for the areas with no availability of shipborne gravity measurements. Therefore, a marine CBA map around Greece (Fig. 3) has been produced by the integration of the gravity data provided by shipborne gravimetric observations and satellite altimetry-derived gravity data. Shipborne gravimetric observations can provide a significant enhancement of the short-wavelength signals especially in shallow water and coastal zones where the altimetric data accuracy is reduced (Zaki et al., 2022).
In Fig. 3, we can observe the shipborne measurements used in this paper illustrated with grey dotted lines. These lines are covered with quite good density, the western and southern area of Greece (Ionian Sea), the Gulf of Corinth, the sea surrounding Crete Island, the southern Cyclades and part of the Aegean Sea. On the other hand, the gravity information of the remaining areas where there are no shipborne available data, namely the greater part of the Aegean Sea and Northern Ionian Sea, has been obtained from altimetry-derived gravity data. The dataset used, called ''NTUA-FAAv1'' (Mintourakis, 2020), has been produced from these data, providing the Free-Air Anomaly map of these areas. The specific data have been processed with the ''Iterative subarc-levelling technique'' (Mintourakis, 2020), which seems to offer a slight improvement to the gravity modelling close to the coastline, regarding mainly the short-wavelength gravity signals.
At this point we have to mention that all the shipborne gravity data have been retrieved from international organizations, such as Ifremer, National Oceanic and Atmospheric Administration (NOAA/ NCEI), Marine Geoscience Data System (MGDS), Lamont-Doherty Earth Observatory (LDEO) and Bureau Gravimetrique International (BGI)/IAG International Gravity Field Service. In the end, the total amount of gathered marine gravity data was almost 3,675,250 points. Taking into account that all these were acquired during different projects, ships, periods, guidelines and organizations along with the fact that may have been processed with different formulas for data reduction, we had to run an extensive quality control. All these shipborne data were projected in a common GIS environment categorized by their source. The first step was thorough optical check, trying to locate adjacent stations, from different sources, that could have very different values. Some of them could be referred to Potsdam gravity system and others to IGSN'71. Beyond that, different reduction formulas may have been used and for that reason we had to search each dataset's specifications and parameters. This kind of problematic data could be corrected by re-calculating the gravity values. The purpose was to build a common homogenous marine gravity dataset, determined by the same type characteristics, such as data reduction formulas, reference system and coordinate system. Afterwards, gravity maps of these marine data had been constructed in order to locate problematic individual gravity points and biases. These could be revealed from severe contouring problems caused for example by wrong digits of the gravity value. Therefore, after this timeconsuming process we ended up rejecting a great amount of the marine gravity points due to errors, problems, doubtful data, or even unknown processing steps. Sometimes we even had duplicate datasets, retrieved from the different organization/databanks but with different values or even coordinates. For that reason, all the gravity points that we used, were referred to the IGSN'71 datum (Morelli et al., 1974).
The final marine gravity dataset after all the quality control, is comprised of almost 2,989,700 marine gravity points, based on marine gravity data (Armijo et al., 2002;Ewing, 2015;Foucher, 1992;Gutscher, 2002;Hooft et al., 2015;Le Pichon & Chamot-Rooke, 1995;Mascle, 1998;Rawson, 2014;Sofianos, 2015;Taylor, 2014) but also data retrieved from the databases of Bureau Gravimetrique International/BGI (Doi: 10.18168/BGI) and https:// campagnes.flotteoceanographique.fr/search. The rest of the data reduction process, beyond the initial correction applied from the organizations (e.g. Eotvos correction) for the production of the Sea CBA map ( Fig. 3), was based on the formulas used for the terrestrial data (see Sec 2.1). Its minimum Bouguer gravity value equals to -155.42 mGal and its maximum is 279.59 mGal. The greater gravity anomalies are located in the west of the Ionian Islands (Ionian Sea), where we have great depths (Fig. 1) due to the existence of the Hellenic trench. Additionally, we have a cluster of great gravity values at the area of southern Aegean Sea.

Satellite-Only GGMs
The satellite data used in this paper were basically originated from the satellites that were launched a few years ago by European Space Agency (ESA) and National Aeronautics and Space Administration (NASA), known as GOCE and GRACE correspondingly. Both of these satellites and their satellite-only GGMS have been used by several researchers globally for several geological purposes (Alvarez et al., 2015;Braitenberg, 2015;Choi et al., 2006;Fuchs et al., 2016;Novák et al., 2019;Pal & Majumdar, 2015;Vaish & Pal, 2015).
GOCE is ESA's satellite mission that combines gravity gradiometry and GPS tracking to determine the Earth's mean gravity field. The GOCE gravity gradient data are delivered from the on-board gradiometer. The gravity gradients are derived from the differences among the accelerations of three pairs of accelerometers (Rummel et al., 2011). The current effective spatial resolution of GOCE data is around 80-100 km at the Earth's surface (van der Meijde et al., 2015). The mean height above the ellipsoid is almost 270 km. This contributes to the improvement of gravity data in the medium wavelength range, between 75 and 200 km spatial resolution and has its largest contribution over regions where sparse or lowquality terrestrial data are available. The GOCE anomalies are composed of the low and the middle frequencies.
On the other hand, the GRACE mission launched in March 2002, consisting of two identical satellites that map earth's gravity field, by making accurate measurements of the distance between the two spacecrafts (Tapley et al., 2004). Their continuously provision of data ended somewhere in October 2017. Nowadays, the GRACE-FO satellites (launched on 22 May 2018) continue their scope (Boergens et al., 2020). Based on the fact that GRACE data are provided in monthly datasets, we can also have some indications for mass changes and crustal deformation of the area (Li et al., 2020).
The satellite-only GGMs that have been used in this paper are GOCO06S (Kvas et al., 2021) and ITU_GGC16 (Akyilmaz et al, 2016). Their characteristics and accuracy statistics can be found in Table 1. Kvas et al. (2021) state that the GOCO06s model provides a consistent combination of spaceborne gravity observations from several satellite missions and measurement techniques. The long to medium spatial wavelengths are mainly determined by GRACE due to the high sensitivity of the intersatellite ranging observation, while the kinematic LEO orbits mainly contribute up to degree and order 150. Finally, the medium to short wavelengths of the solution, beginning from degree 120, are mostly provided by the GOCE gradiometer observations. ITU_GGC16 is a static global gravity model up to degree order 280 computed from the combination of ITU_GRACE16 (up to d/o 180) and GO_CONS_GCF_2_TIM_R5 (up to d/o 280). The integration is performed at the normal equation level along with variance component estimation.
The prementioned satellite models were downloaded by the International Centre for Global Earth Models, known as ICGEM. It provides a helpful web interface to calculate and download the computed values of gravity field functionals of the Earth's global gravity field, on freely selectable grids, with respect to a user-defined calculation grid (Ince et al., 2019). Taking advantage of this online unique Table 1 Characteristics and statistics for the used, satellite-only, Bouguer anomaly models regarding the study area

Model
Year service, we have calculated and downloaded the required data of the gravity anomaly distribution for the entire coverage of Greece, with a selected spacing of 0.05°(36,381 points) and mean sphere radius equal to 6,378,137 m. High-level processing facilities carry out the initial corrections on the satellite raw data up to level 2 product (Dahle et al., 2013;Frommknecht et al., 2011;Gruber et al., 2007) that are being distributed to the public. Beyond these initial corrections, the rest of the necessary corrections needed for the calculation of the satellite CBA maps (Figs. 4, 5), were based on the formulas and parameters used for the terrestrial data (see Sec 2.1). In Fig. 4 we can see the CBA maps based on the GOCO06S satellite model (Kvas et al., 2021) and in Fig. 5 the one for satellite maps are quite smooth, due to their type of data and their spatial accuracy, which is not as good as the one from the terrestrial and marine gravity data.

Combined Gravity Model
Due to the reduction of the gravity field with orbit altitude, the satellite-only GGMs will never be able to reach very high spatial resolutions of only a few kilometers. However, precise knowledge of the Earth's gravity field structure with very high resolution is essential not only for a range of geo-science disciplines, such as solid Earth geophysics for lithospheric modelling and geological interpretation, exploration geophysics but also for geodesy (e.g., surveying, inertial navigation) and civil engineering (e.g., construction, modelling of water flow for hydro-engineering). For this reason, gravity models derived only from satellite data are complemented by combined gravity models, which contain much higher resolution gravity information retrieved from terrestrial gravity measurements. Several global combined GGMs have been presented in the literature so far (e.g., Fecher et al., 2017;2005). The combined models have commonly higher resolution than the satellite-only ones, which could be down to 10 km or better.

Least-Squares Collocation
For the integration of the terrestrial, marine and satellite datasets presented in this paper, we used the procedure of least-squares collocation (LSC), which has been used successfully by several authors in the past for gravity field modeling (Belay et al., 2021;Kamto et al., 2021;Ramouz et al., 2019;Saadat et al., 2019;Sadiq et al., 2010;Sobh et al., 2019;Yildiz, 2021).
As Moritz (1978) states, LSC is practically a mathematical technique for determining the Earth's figure and gravitational field by a combination of heterogeneous data of different kinds. LSC is a useful procedure that provides information about the errors along with the estimated regional gravity field in forms of variances or full variance-covariance matrix. The quality of the inversion of gravity data is highly depended on the covariance function of the observations, which is the observed gravity signal (Knudsen, 1987). Based on Moritz (1980), LSC is basically a statistical estimation technique combining least square adjustment and least square prediction into a linear regression model. But as Sadiq et al. (2010) also state, LSC practically minimizes the errors providing an outcome with a smoother solution related to the one of the spherical harmonics, due to the downward continuation (Moritz, 1980). The LSC has been primarily introduced by Krarup (1969) and then was thoroughly analyzed by Moritz (1973Moritz ( , 1978 and Grafarend (1976). As Moritz (1978) states, we can try to estimate T (disturbing potential) by a linear combination f of applicable base functions u 1 , u 2 , …., u q : where practically P represents the given space point, for which the above functions are being considered, and b k the appropriate coefficients. If we have a number of errorless values of the disturbing potential (T), at n space points P i (Eq. 6), we can assume that the calculation f to T is provided with their interpolation (Eq. 7): We can generalize the above interpolation equations considering n values of linear functionals L 1 T, L 2 T, …., L n T, such as gravity anomalies given by Eq. (8), where xyz represents coordinates.
This procedure of fitting an analytical calculation to n given functionals is called collocation. On the other hand, the rms interpolation error (m P ) of the basic interpolation function u i (P) is: where M is an appropriate average and e P is the local interpolation error that is defined as: Taking into account that we have the minimum error at the collocation procedure (m P = min), we will have the following equation: where L Q k defines that the linear functional L K is appliable to variable Q. Now, if we put: We will then have the following LSC formula At this point we must take into account the fact that the observed functionals l i have quite often measuring errors, considered as noise (n i ) that have to be included in the collocation Eq. (13). Therefore, if the errorless part of the functional is denoted as s i we have: We will then have the signal covariance matrix [C ij ] stated by Eq. (15). If the corresponding error covariance matrix is [D ij ], it will result in the Eq. (16) representing the total covariance matrix of the functionals: As it is clarified in Moritz (1978), a best statistically linear estimate is provided, while the measuring errors do not affect the signal.

Application of Least-Squares Collocation
The LSC in the present work has been executed in the MATLAB environment, using the griddataLSC script (MATLAB Central File Exchange, 2020) that carries out the data interpolation, empirical covariance estimation and function fitting. The particular script offers the choice of 6 covariance functions, such as exponential (Shaw et al., 1969), Reilly (Collier, 2002), Gaussian (Kearsley, 1977), logarithmic (Forsberg, 1987 and 2nd order Markov (Kasper, 1971).
The authors carried the LSC with all provided functions of griddataLSC script, through the MATLAB environment, but ended up using the results of 2nd order Markov, as the most satisfactory ones. This was based on the production of the relevant maps, necessary to run visual checks but also compare their characteristics (e.g., minimum and maximum values). Therefore, the use of Gaussian and log3D functions reduced the minimum values of the results a lot (-132.78 and -154.90 mGal correspondingly) or the maximum ones (264.23 and 296.30 mGal correspondingly). On the other hand, Reilly function provided very sharp transitions of the values and reduced maximum values (277.24 mGal). Exponential and triangulation functions provided similar results to the 2nd order Markov. Beyond that, the covariance function parameters (C 0 and D) that have been determined were lower for the 2nd order Markov. The applied LSC procedure with the 2nd order Markov function calculated the combined regional CBA map of Fig. 6.
There are several studies that have also preferred the application of the 2nd Markov function for gravity modelling (e.g., Kamto et al., 2021;Märdla et al., 2017;Oja et al., 2019). This type of covariance function seems to be better than other ones, based on the fact that it encounters for both direction and distance among the data points, which is important if the spatial variability is anisotropic. This will provide a more accurate modelling, especially for large variations of the gravity field, like expected in Greece. Beyond that, during the 2nd order Markov function there is a spatial autocorrelation of nearby data points for the short-range correlation while for the long-range one we have large-scale patterns Vol. 180, (2023) Regional Gravity Model of Greece (Oliver & Webster, 2015). Kasper (1971) also proposes a second-order Markov process as a statistical model for gravity anomalies in a local region.

Results and Discussion
The combined regional gravity map of Greece that has been developed (Fig. 6), in terms of CBA map, is the result of LSC with 2nd order Markov function, integrating the terrestrial, marine and satellite gravity data of the Greek territory. The minimum observed gravity value is -178.17 mGal and the maximum is 290.23 mGal (Table 2). Generally, the lower gravity values are located across the Pindos range, which is located at the western mainland but also at the area of eastern Macedonia. The values of the lower gravity anomaly at the combined model have been reduced compared to the ones of the terrestrial data (Fig. 2) for the same area (down to -251.01 mGal). This is probably caused due to fact that this area was also determined by satellite data (Figs. 4,5) with quite smaller gravity anomaly values, down to -128.09 mGal the minimum. On the other hand, the higher gravity values (up to 290.23 mGal) are located at the Ionian Sea (Fig. 6) and at the southern Aegean Sea, between the Cyclades and Crete Island. These two areas, are determined with similar ranges of gravity values at both marine (Fig. 3) and satellite (Figs. 4,5) data, which are ranging up to 279.60 and 302.50 mGal correspondingly (Table 2). Beyond these datasets, there is also a small impact from the terrestrial data of the Aegean islands (Cyclades), that reach in the terrestrial CBA map (Fig. 2) values up to 156.80 mGal. Therefore, the results of the regional combined gravity model (Fig. 6) seem to be satisfactory, without problems.
The relatively sharp transition from higher gravity anomaly to lower at the southeast of the area (Ionian Sea) is probably determined from the existence of the Hellenic trench, extending along this area. The higher gravity values at the Aegean Sea, between Cyclades and Crete are due to the existence of the Aegean microplate there.
The authors believe that the new complete gravity anomaly (CBA) map of Greece, presented in Fig. 6, has two advantages related to the existing ones. At first, the individual gravity datasets used (terrestrial, marine and satellite ones) are the updated and revised, compared to the older ones and especially the satellite ones. Beyond that, the compilation has been conducted with the use of LSC, which is considered to be an efficient procedure to combine such heterogeneous data. In Table 3 a summary of the statistics regarding the existing Bouguer maps, presented in the introduction, along with the ones of the new complete Vol. 180, (2023) Regional Gravity Model of Greece 2819 gravity map (Fig. 6) are presented. At this point, we have to mention that given the fact that we could not retrieve the original statistics for some of the existing maps, we provide the calculated ones based on their digitization by the authors. Based on Table 3, the new CBA map has the lowest minimum value, while its maximum one is close to the one of Makris et al. (2013). Mean values are quite similar for all except for Grigoriadis et al. (2016), which has the lowest one. Finally, regarding the standard deviation of the new map is lower than the one of Makris et al. (2013) and higher than the ones of Grigoriadis et al. (2016) and Lagios et al. (1988). The newly developed combined gravity model of Greece, in terms of CBA (Fig. 6), has been compared to the input gravity datasets, that is the terrestrial HMGS terrestrial data, the marine dataset (comprised of the shipborne and altimetry-derived gravity ones)  Table 4). The STD differences of 10.47, 15.36, 11.70 and 10.60 mGal were calculated between the combined gravity model and the terrestrial (Fig. 7), marine (Fig. 8) and satellite datasets correspondingly.
More specifically, in Fig. 7, the CBA differences between the LSC and terrestrial-only model values have been illustrated. The LSC model seems to have higher CBA values (positive difference values) at areas of Greece (up to 72.80 mGal) that could be considered mountainous, such as across the Pindos range and Crete. This could be explained by the sharp topography that affects the terrestrial data. On the other hand, lower CBA values (negative difference values) have been revealed across flat areas at the LSC CBA map (down to -39.20 mGal). These differences could be due to the resolution difference between the terrestrial-only and LSC models. In Fig. 8, the CBA differences between the LSC and marine-only model values have been illustrated. At this point we can observe even higher CBA values (positive difference values) from the LSC CBA model (up to 153.50 mGal), compared to the one of marine-only measurements. It is quite obvious that these values are located across the Hellenic trench, that seem to affect them. On the other hand, lower CBA values (negative difference values) have been revealed mainly at the Aegean Sea (down to -54.60 mGal), between Cyclades and Crete, probably due to the existence of the Aegean microplate there.
Beyond the comparison between the initial gravity datasets and the produced LSC CBA model that have been discussed in the previous paragraphs, an additional comparison is also presented in Fig. 9. At  Makris et al. (2013). The CBA values of LSC model are higher (positive differences) mainly across the Hellenic trench (up to 71.25 mGal), that probably due to the marine datasets that have been used for both models. Accordingly, the LSC CBA model seems to have lower values (negative differences down to -75.96 mGal) at the Aegean area. This is similar to the comparison presented in Fig. 9, between the LSC and marine-only CBA values. But there are also lower CBA values defined in LSC model at the greater part of Greece terrestrial land. This could be due to the evolvement of quality of the satellite data that have been used between the two CBA models.

Conclusions
A new combined regional model of Greece and its surrounding area have been developed based on the combination of terrestrial, marine and satellite datasets. For the appropriate combination of these datasets the Least-square collocation has been applied due to their inhomogeneity. It is a regional gravity model of Complete Bouguer Anomaly values that can be valuable for regional geological and geophysical studies of Greece, since it provides homogeneously the distribution all over Greece, with updated gravity data.
The collocation procedure has to be done with great care, after running the appropriate checks of all the datasets. It is a relative time-consuming process, which is though essential for the homogeneity of the final CBA model.
The comparison of the combined regional gravity model that has been developed with the initial input datasets revealed that they are in good agreement with relatively low STD differences, ranging from 10.50 to 15.40 mGal, depending on the type of data. Based on the comparison maps between the initial datasets and the LSC model we found out that there are some geological structures that seem to affect them, such as the Hellenic trench, the Pindos range and the Aegean microplate. At these areas the most significant differences have been observed, revealing that it is important for a CBA model to include all types of gravity datasets (terrestrial, marine and satellite).
Finally, from the comparison of the new LSC model with the existing ones, the necessity of the continuous update of the gravity datasets have been revealed in order to maintain the best possible quality of the result.

Data availability
Original data can be retrieved from the corresponding public data repositories, as provided in the text.

Declarations
Conflict of interest The authors have not disclosed any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.