1 Introduction

Geoscientists frequently use three sources of information to investigate the Earth’s interior structure. The first set is understood as direct evidence from rock samples by drilling projects, the second set includes the records of seismic waves, which are generated, for example, by earthquakes, explosions, volcanoes and other natural or anthropic sources, and the third set of information is the gravity field models generated through modern satellite gravity missions such as Challenging Mini-satellite Payload (CHAMP), Gravity Recovery and Climate Experiment (GRACE) and Gravity field and steady state Ocean Circulation Explorer (GOCE), which can provide global and homogeneous coverage of data. A huge improvement can also be obtained in the accuracy and spatial resolution of these models by combining them with airborne and ground-based gravity data as well as satellite altimetry observations over the oceans (see e.g. Hamayun 2014). Other important sources for studying Earth’s interior are meteorites and its magnetic field.

The primary/uppermost interface of the Earth’s interior is the Mohorovičić discontinuity (or Moho), the boundary between the Earth’s crust and upper mantle. Typically today, one uses both gravimetric and seismic methods to investigate the depth of this surface. In case of the seismic method, Moho is identified by the travel time of the seismic wave reflected at the Moho boundary, whereas in the gravimetric method, gravity data are used under a certain isostatic hypothesis, which assumes isostatic equilibrium of the crust on the dense underlying mantle (e.g., Abrehdary 2016). However, due to the high cost seismic data acquisition and lack of global coverage, its application for the recovery of the Moho interface is limited. Therefore, in the absence or low coverage with seismic data, gravimetric or combined gravimetric-seismic methods may improve the result. Classical isostatic approaches, as Airy/Heiskanen (A/H), Pratt/Hayford (P/H), and Vening Meinesz (VM) models, which mainly depends on the excess/deficit of topographic (solid) mass over continents/sea, respectively, are typically not able to model realistically the actual Moho geometry, because the isostatic mass balance depends on global mass loading and effective elastic thickness, rigidity, rheology of the lithosphere and viscosity of the asthenosphere (Watts 2001; Oruç et al. 2019). Moreover, the glacial isostatic adjustment (GIA), present-day glacial melting, plate motion, mantle convection and other geodynamic processes contribute to the overall isostatic balance. The classical isostatic models assume local (A/H and P/H) or regional (VM) isostatic balance, while Moritz (1990) generalized the VM model to a more realistic global compensation with a spherical sea level approximation, as expressed by Sjöberg (2009) in mathematical form, who presented a new solution according to Moritz theory under the name of Vening Meinesz-Moritz (VMM) model.

Many studies on Moho depth determination in Fennoscandia have been published. For example, Bungum et al. (1980) estimated the crustal thicknesses for 11 Fennoscandian seismograph stations equipped with three-component long period instruments, using the so-called spectral ratio technique of Phinney, and they demonstrated that the largest Moho depths are found for stations located in the north-east areas of Scandinavia and Finland, with a local maximum in the Bothnian Bay. Luosto (1997) investigated the structure of the Earth’s crust in the region as revealed from refraction and wide-angle reflection studies and he showed that the thickest values are found in the middle of the shield and in Central Finland. Bagherbandi and Sjöberg (2012) studied the non-isostatic effects (NIEs) on crustal thickness and concluded that after correcting the gravimetric–isostatic model for the NIEs, the final result for the Moho can be considered as bias-free. Silvennoinen et al. (2014) determined a Moho boundary map for the northern Fennoscandian Shield based on combined controlled-source seismic and receiver function data, and they observed large Moho depths in the northern part of the study area, deepening towards the northeastern corner. Bagherbandi et al. (2015) presented another gravimetric–isostatic crustal thickness model for Fennoscandia, showing how the main geological unit structures attribute/contribute in isostatic balance by affecting the Moho geometry.

The primary aim of this study is to determine a new Moho depth model (named HVMD19) for Fennoscandia, based on gravimetric–altimetric and topographic information using the VMM isostatic model and a seismic crust model. The seismic model is used to determine the mean Moho depth and Moho density contrast and for correcting the gravity data for NIEs. For this purpose, the data files of the recent XGM2019e gravimetric-altimetric model (Zingerle et al. 2019), Earth2014 Earth topographic/bathymetric model (Hirt and Rexer 2015) are used along with CRUST1.0 (Laske et al. 2013) crustal seismic model to a regional Moho depth model over the area with a resolution of 1° × 1°. As the observed gravity data set suffers disturbing signals from the topography and density variations related to bathymetry, ice, sediments, and other crustal components, we remove these signals by applying so-called stripping gravity corrections (see Tenzer et al. 2015). However, there are also other, mainly deeply located, mass inhomogeneities, which to large extent can be removed as gravimetric so-called NIEs, as a procedure to filter out as much as possible of the remaining signals in the gravity data caused by disturbing effects in the Earth’s interior by comparing the gravimetric and seismic data in the frequency domain. The NIEs could be attributed to different mechanisms, such as crustal thickening/thinning, thermal expansion of mass of the mantle (Kaban et al. 2004), Delayed Glacial Isostatic Adjustment (DGIA), and plate flexure (Watts 2001, p. 114). By the effect of DGIA we mean the expected integrated effect of future GIA. However, there are also other causes of NIEs that makes the Moho topography inconsistent with an isostatic model (Bagherbandi and Sjöberg 2012). Hence, a filtering technique is used for removing the remaining NIEs not possible to remove by the above stripping technique. Note that the Moho geometry can also be determined without an isostatic assumption (e.g. Oldenburg 1974). However, also in this case density heterogeneities in the crust and upper mantle should be stripped, therefore leading to the same problem.

In Fennoscandia, the crust is currently rising, because of the DGIA as a response of the viscous mantle to melting of the Late Pleistocene ice sheet. As a result, there is a negative anomaly in the present-day static gravity field as isostatic equilibrium has not yet been reached (e.g., Bjerhammar et al. 1980; Root et al. 2015). In this context, our second aim of the paper is to correct the specific NIEs of the gravity field for the DGIA directly prior to the general NIE. In other words, the DGIA effect will be determined and removed separately. Then two questions arise: (1) Can the DGIA gravity signal be isolated and removed from the observed gravity signal, and (2) does this specific correction significantly improve the Moho depth determination when considering that much of this effect is already part of the NIE? Question (1) is dealt with in Sect. 3.3 as part of the description of all additive corrections to gravity as described in Chap. 3, and question (2) is answered in the Chap. 5, where the results of our computations are presented and discussed. However, first we present the basic formulas for the VMM solution in the next chapter. In Chap. 4 the data are described, and Chap. 6 concludes the paper.

2 Altimetry–Gravimetry Moho Recovery

2.1 The Basic Equation of the VMM Method

The isostatic gravity disturbance \(\delta g_{I}\) is the sum of the refined gravity disturbance \(\delta g_{B}^{R}\), reduced for the gravitational contributions of the topography and known anomalous crustal density structures, and the isostatic compensation attraction \(A_{C}\) (Sjöberg 2009, 2013), which in case of isostatic equilibrium (i.e. total compensation) vanishes:

$$\delta g_{I} (P) = \delta g_{B}^{R} (P) + A_{C} (P) = 0.$$
(1)

This formula can be re-formulated as a non-linear integral equation for the VMM problem of isostasy:

$$R\iint_{\sigma } {K\left( {\psi ,s} \right)d\sigma = - \left[ {\delta g_{B}^{R} (P) + A_{{C_{0} }} (P)} \right]}/G\Delta \rho = f(P),$$
(2)

where G is the Newton's gravitational constant, R is the radius of the mean Earth sphere, \(\Delta \rho\) is the Moho density contrast, \(\sigma\) is the unit sphere, K is the integral kernel function of the spherical distance \(\psi\) and the parameter \(s\left( { = 1 - D/R} \right)\) (s is a simple function of the D (Moho depth), which is the unknown of the integral equation). The spectral representation of K reads:

$$K\left( {\psi ,s} \right) = \sum\limits_{n = 0}^{\infty } {\frac{n + 1}{{n + 3}}\left( {1 - s^{n + 3} } \right)} P_{n} (cos\psi ),$$
(3)

where \(P_{n} (cos\psi )\) is the Legendre polynomial of degree n and \(A_{{C_{0} }} \left( P \right)\) is the nominal compensation attraction, which can be derived from an a priori crustal seismic model.

2.2 The VMM Solution for the Moho Depth

Here, we use the method presented by Sjöberg (2009) for computing the Moho depth. To do so, the function \(f(P)\) in Eq. (2) is rewritten as a harmonic series:

$$f\left( P \right) = \sum\limits_{n = 0}^{\infty } {\sum\limits_{m = - n}^{n} {f_{nm} } } Y_{nm} \left( P \right)$$
(4)

with

$$f_{nm} = \frac{1}{4}\iint\limits_{\sigma } {f(Q)Y_{nm} (Q)d\sigma }$$
(5)

and

$$Y_{nm} (Q) = \left\{ \begin{gathered} \overline{P}_{nm} (\cos \theta )\cos (\left| m \right|\lambda )\quad m \ge 0 \hfill \\ \overline{P}_{nm} (\cos \theta )\sin (\left| m \right|\lambda )\quad m < 0 \hfill \\ \end{gathered} \right.,$$
(6)

where \(Y_{nm} (Q)\) is the fully-normalized spherical harmonic of degree n and order m and θ and λ are co-latitude and longitude.

By rewriting also the left-hand side of Eq. (2) as a series in spherical harmonics, Sjöberg (2009) could, after linearization, identify each harmonic of D as a harmonic in f, and by summing up he received the following first-order solution for the Moho depth:

$$D_{1} \left( P \right) = \sum\limits_{n = 0}^{\infty } {\left( {2 - \frac{1}{n + 1}} \right)} \sum\limits_{m = - n}^{n} {f_{nm} Y_{nm} \left( P \right)} ,$$
(7)

and it is used as input to Eq. (2) for the second-order formula:

$$D\left( P \right) = D_{1} \left( P \right) + \frac{{D_{1}^{2} \left( P \right)}}{R} - \frac{1}{32\pi R}\iint_{\sigma } {\frac{{D_{1}^{2} \left( Q \right) - D_{1}^{2} \left( P \right)}}{{\sin^{3} \left( {\psi /2} \right)}}d}\sigma_{Q} .$$
(8)

3 Additive Corrections to Gravity Disturbance

3.1 Stripping Gravity Corrections (For Crust Density Variations)

The main input data set used for VMM Moho depth determination is the the refined Bouguer gravity disturbances, i.e. free-air gravity disturbances corrected for topography, bathymetry, ice thickness and sediment basins (or stripping gravity corrections). To achieve this, we apply a uniform mathematical formalism of computing the gravity corrections \(\delta g^{TBIS}\) of the density variations within the Earth’s crust by Tenzer et al. (2015):

$$\delta g^{TBIS} = \delta g^{T} + \delta g^{B} + \delta g^{I} + \delta g^{S} ,$$
(9)

where \(\delta g^{T}\) is the topographic gravity correction, and \(\delta g^{B}\), \(\delta g^{I}\) and \(\delta g^{S}\) are the stripping gravity corrections due to the ocean/bathymetry, ice and sediment density variations, respectively. Then the refined Bouguer gravity disturbance \(\delta g_{B}^{R}\) is realized from the XGM2019e free-air gravity disturbance \(\delta g_{FA}^{{{\text{XGM2019e}}}}\) by:

$$\delta g_{B}^{R} = \delta g_{FA}^{{{\text{XGM2019e}}}} - \delta g^{TBIS} .$$
(10)

3.2 NIEs

In reality the crust is not in complete isostatic equilibrium, and the observed gravity data are not only generated by the topographic/isostatic masses but also from masses in the deep Earth interior, that leads to non-isostatic effects (NIEs; Bagherbandi and Sjöberg 2012, 2013). According to Sjöberg (2009) the major parts of the long-wavelengths of the gravity field are due to density variations in the Earth’s mantle and core/mantle topography variations, and these parts thus belongs to the NIEs. As already mentioned, the NIEs could be the contribution of different factors, such as crustal thickening/thinning, thermal expansion of mass of the mantle, DGIA, plate flexure, and effect of other phenomena, implying that this contribution to gravity will lead to systematic errors/NIEs of the computed Moho topography. Hence, also the refined Bouguer gravity disturbance, used for Moho depth determination, must be corrected/reduced for the NIEs.

Assuming that the seismic Moho model CRUST1.0 is known and correct, the gravity effect of the NIEs can be estimated from its difference to the erroneous preliminary VMM gravity effect:

$$\delta g^{NIE} = \frac{GM}{{R^{2} }}\sum\limits_{n = 0}^{{n_{\max } }} {\left( {n + 1} \right)} \sum\limits_{m = - n}^{n} {c_{nm}^{NIE} Y_{nm} \left( P \right)} ,$$
(11a)

where

$$c_{nm}^{NIE} = c_{nm}^{{{\text{CRUST1}}{.0}}} - c_{nm}^{VMM} .$$
(11b)

Here \(\begin{array}{*{20}c} {c_{nm}^{NIE} ,} & {c_{nm}^{VMM} ,} & {c_{nm}^{{{\text{CRUST1}}{.0}}} } \\ \end{array}\) are the spherical harmonic coefficients of the gravity disturbances of the NIE, VMM and CRUST1.0, respectively.

As we already mentioned, this study emphasizes on investigating the DGIA effect on the gravity disturbance and Moho depth determination in Fennoscandia. To achieve this goal, we divide NIEs into NIE and NIE1, i.e. without and with special removal of the DGIA effect.

3.3 DGIA

Delayed GIA (DGIA) expresses the delayed adjustment process of the Earth to an equilibrium state when former ice sheet loads have vanished. The ongoing adjustment of the Earth's body to the redistribution of ice and water masses is evident in various phenomena, which have been studied to infer the extent and amount of the former ice masses, to reconstruct the sea level during a glacial cycle and to constrain rheological properties of the Earth's interior. Here we aim at answering the question whether the effect of the gravimetric DGIA correction is significant for Moho determination in Fennoscandia.

3.3.1 The DGIA Effect on Gravity Disturbance

Several studies have been performed to strip the GIA related gravity signal from the observed gravity field. However, the conclusions of these studies differ significantly. For example, following an idea by Bjerhammar et al. (1980), who used the spherical harmonic spectral window 10–22, Sjöberg et al. (1991, 1994) stripped the GIA signal from other gravity sources by a band-limit filtering of the harmonic series of the gravity potential both to degrees 10°–22° and 10°–100° in attempts to outscore non-GIA related components of the gravity field. This approach obviously can only be approximate, as the GIA signal will always correlate also with other harmonics as emphasized, e.g. by Kakkuri and Wang (1998) and Root et al. (2015). Sjöberg et al. (1991, 1994) estimated the GIA related gravity signal by using harmonic window 10–100 to about − 30 mGal in the central uplift region, and by also using regression analysis in 1994 they got − 28 mGal, which agrees well also with the minimum received by Root et al. (2015, Fig. 6a), although admittedly the location of the minimum is somewhat shifted. Here, we will use the spectral window 10–23, which yields the highest correlation (0.92) with the recent NKG2016LU GIA model (Vestøl et al. 2019), in close agreement with Bjerhammar et al. (1980) and Joud (2018).

Usually the DGIA effect is part of the general NIE, but here we aim at specifically determine and remove its effect prior to estimating the remaining NIE. The purpose for this is twofold: firstly, we want to directly estimate the DGIA effects on gravity and Moho depth, and, secondly, we want to investigate how efficient the NIE procedure is to remove the DGIA effect.

The first task will be answered as follows. By considering the spectral window that provides the optimum correlation between the GIA and gravity (cf. Sect. 3.3.1), one may estimate the DGIA effect on gravity as:

$$\delta g^{DGIA} = \gamma \sum\limits_{n = 10}^{23} {\left( {n + 1} \right)\sum\limits_{m = - n}^{n} {A_{nm} Y_{nm} } } ,$$
(12)

where \(\gamma\) is normal gravity, \(Y_{nm}\) and \(A_{nm}\) are spherical harmonics and coefficients of the gravitational potential, and the resulting effect on the Moho depth is approximately given by (Sjöberg 2009, Eq. 28):

$$D^{GIA} \approx - \frac{{\delta g^{DGIA} }}{2\pi G\Delta \rho }.$$
(13)

Hence, if \(\delta g^{DGIA}\) is within − 30 mGal, \(D^{DGIA}\) will hardly exceed − 1.5 km, which will be further analyzed in Sect. 5.4.

4 Description of the Data

In order to determine the Moho depth model for Fennoscandia, named HVMD19, bounded by latitudes (55° N–71° N) and longitudes (5° W–31° E), we use different data sources. In the following, we further utilize the XGM2019e gravitational model (Zingerle et al. 2019) on a 1° × 1° arc-deg spherical grid with a spectral resolution complete to a spherical harmonic 180° to produce the free-air gravity disturbance along with the Earth2014 topography model and seismological model of CRUST1.0. XGM2019e is a combined global gravity field model represented through spheroidal harmonics up to d/o 5399, corresponding to a spatial resolution of 2 min. The model is formed by merging terrestrial, satellite-only, and satellite altimetry gravity data. We then correct the gravity disturbances for the density variation of the oceans, ice sheets and sediment basins (i.e. stripping gravity corrections) according to Tenzer et al. (2015), the NIEs and also DGIA, (see Eqs. 9–12). As already mentioned, Fennoscandia is an interesting area for GIA research, and we therefore compute the NIE in two ways: Case (1) by Eqs. (11a) and (11b), and Case (2) by first determining \(\delta g^{DGIA}\) as in Fig. 1, removing it from the NIE calculated in case 1 and computing a new NIE. Finally these NIEs with and without specific removal of the DGIA effect are added to the refined Bouguer gravity disturbances. CRUST1.0 is used to determine the mean Moho density contrast and in the determination of the NIEs. Also, CRUST1.0 model is utilized to obtain the mean Moho depth for each block. In oerder to validate our Moho depth estimates, the CRUST19 model (Szwillus et al. 2019) is first considered, which displays global maps of Moho depth and average P wave velocity of the crystalline crust with a resolution of 1° × 1° using a nonstationary kriging algorithm, including the interpolation error and error covariance and then the GRAD09 model is chosen, which is compiled from more than 250 data sets of individual seismic profiles, 3-D models obtained by body and surface waves, receiver function results and maps of seismic and/or gravity data compilations in the European Plate (Grad et al. 2009). As already explained in Sect. 3.3.1, we use the spectral representation of Eq. (12) in estimating the DGIA signal on gravity. The statistics of the estimated refined Bouguer gravity disturbance, NIE and DGIA are given in Table 1.

Fig. 1
figure 1

The DGIA effect on gravity by Eq. (12) with \(n_{\max } = 23.\) (Unit: mGal)

Full size image
Table 1 Statistics of the free-air and refined Bouguer gravity disturbance estimated from XGM2019e model
Full size table

Using Eq. (12) for estimating the DGIA effect one arrives at Fig. 1, which is in fair agreement with Root et al. (2015). From Table 1, one can see that the minimum contribution to gravity is − 19.6 mGal.

Figure 2 maps the free-air gravity disturbance computed on a 1° × 1° grid using XGM2019e coefficients complete to 180° of spherical harmonics, showing large positive disturbances in the southern Norwegian high mountain massive, and locally in the north Norwegian sea, whereas the largest negative anomalies are seen in the Tromsö region in Norway, as well as in the Baltic Sea and Gulf of Bothnia, the latter most likely due to the DGIA effect (cf. Root et al. 2015).

Fig. 2
figure 2

The free-air gravity disturbance. (Unit: mGal)

Full size image

Figure 3 images the Refined Bouguer gravity disturbance (reduced for gravity of topography, density heterogeneities related to bathymetry, ice, sediments, and other crustal components by applying stripping corrections) and corrected for NIEs (disregarding special correction for the DGIA effect). As one would expect, these corrections have been able to change the small free-air gravity disturbances in Fig. 2 to large variations. Particularly notable are the large negative disturbances in the Baltic Sea and eastern Finland.

Fig. 3
figure 3

The estimated refined Bouguer gravity disturbance disregarding special correction for the DGIA effect. (Unit: mGal)

Full size image

5 Results for MD and Comparisons

This section presents and interprets our (preliminary and final) Moho depth results before and after correcting for DGIA effects, named PHVMD19 and HVMD19, and it provides comparisons with the CRUST19 and GRAD09 models. Finally, the result of the gravimetric Moho depth determination is interpreted in the context of geological and tectonic setting of the study area.

5.1 The Moho Models PHVMD19 and HVMD19

Here, the refined Bouguer gravity disturbances described in Sect. 4 are used to estimate the Moho depth. For this purpose, as already showed in Sect. 2.2, we assume that the Moho density contrast is known, e.g., from seismic observations, and the Moho is assumed to be in isostatic equilibrium.

The results of the estimated model HVMD19 and its difference with PHVMD19 are shown in Fig. 4a, b, and Table 2 gives some complimentary statistics.

Fig. 4
figure 4

a The estimated HVMD19. (Unit: km). b Difference between PHVMD19 and HVMD19 model. (Unit: km)

Full size image
Table 2 Statistics of the Moho depths estimated from XGM2019e model
Full size table

The figures display small differences (within 1.4 km), but one can see that the maximum Moho depths in the Gulf of Bothnia and in the Central Finland decreases from 58 to 56 km due to the special correction for DGIA effect. Small Moho depths are generally observed in the Norwegian Sea, increasing to largest values in the Gulf of Bothnia (56.5 km) and Central Finland (56 km). Table 2 shows that HVMD19 is somewhat smoother than PHVMD19.

5.2 Comparison with Other Moho Models

In Table 3, PHVMD19 and HVMD19 are compared with the Moho depths from CRUST19 and GRAD09 seismic models. By simple comparison we find that the RMS differences with these models are 3.8/4.2 km and 3.7/4.0 km when the above strategy for removing the DGIA effect is/is not applied, respectively, and the mean value differences are 1.2/1.4 km and 0.98/1.4 km, respectively. Hence, the specific correction for the DGIA effect on gravity disturbance is small, but significant in the central uplift region.

Table 3 Statistics of the differences between the PHVMD19/HVMD19 and Moho depths from the seismic models of CRUST19 and GRAD09
Full size table

Figure 5a, b illustrates the difference between the estimated Moho depths from the preliminary and final HVMD19 and that from the CRUST19 model. As can be seen from this comparison, our result closely agree with CRUST19 in most areas (and differences are mainly within ± 6 km). Taking a closer look at the change from Fig. 5a, b, one can observe that the largest discrepancies, for example those in the Gulf of Bothnia and Baltic Sea, have decreased, which can be attributed to the special correction for the DGIA effect. However, the difference between our model and CRUST19 reaches 10 km along the shoreline of Norway, which to some extent could be explained by the possible use of incorrect crust and upper mantle densities and mean Moho depth in the gravimetric–isostatic model (see Bagherbandi et al. 2015; Abrehdary and Sjoberg 2019). Another explanation could be that (in contrast to most of the inner parts of the Fennoscandia shield) active intraplate tectonic motions occur along the Norwegian coastline (e.g., Ekman 1985), which may cause non-isostatic effects for the Moho constituents. It remains to see how they can be corrected. Furthermore, one can see negative differences within − 8.5 km on the Norwegian Sea, which may be caused by poor interpolation of sparse CRUST19 data.

Fig. 5
figure 5

a Difference between PHVMD19 and CRUST19 model. (Unit: km). b Difference between HVMD19 and CRUST19 model. (Unit: km)

Full size image

Another comparison is performed between the PHVMD19/ HVMD19 and the GRAD09 model in Fig. 6a, b. As in Fig. 5a, b one can see large positive differences in parts of the Fennoscandian shield, especially in most of Norway and parts of central and northern Sweden (i.e. Scandes), demonstrating that the topography of the southern Scandes is not compensated by a mantle thermal anomaly at depth (see Pascal and Olesen 2009). Also, there are significant negative differences in Central Finland and in the Baltic region. This misfit between HVMD19 and the seismic models may be explained by low quality of the latter in this region and/or that the refined Bouguer gravity disturbances in these areas are over-compensated. The other important anomalies are related to the offshore area, where the GRAD09 model estimate an accuracy of about 3 km. These discrepancies should be attributed to the gravity stripping corrections.

Fig. 6
figure 6

a Difference between PHVMD19 and GRAD09 models. (Unit: km). b Difference between HVMD19 and GRAD09 models. (Unit: km)

Full size image

5.3 Geological Setting of Fennoscandia Region

In order to validate HVMD19, we also compare it with the map of geological provinces in Fig. 7, which helps in understanding the relationship between tectonic history and crustal structure. As can be seen from Fig. 4a, the Moho depth is 55–56 km in Central Finland, with a thinning to the east approximately along the Archean–Proterozoic border. The Moho depth in Finnish and Swedish Lapland is at an intermediate level of about 46–50 km, while elsewhere the depth of the Svecofennian crust varies from 40 to 45 km with a thinning towards the western edge of the Svecofennian domain in SW and under the Caledonides and into the Norwegian Sea. An almost E–W oriented zone of thinned (or normal depth) crust is recognized along the Rapakivi formations in Southern Finland, and another area of crustal thinning is seen in the Bay of Bothnia and north-central Sweden. Moving from the shield to the Palaeozoic zone in SW and south, the Moho depth thins to 30–35 km, which is the typical Moho depth in Central Europe. The thinnest Moho can be seen in the NW in the transition from the continental to the oceanic crust.

Fig. 7
figure 7

Map of the geological provinces in Fennoscandia

Full size image

5.4 Moho Depth vs. GIA Effect

The DGIA is a non-isostatic gravity signal, which will falsely affect the Moho depth estimate if not considered. Admittely, it is part of the general NIE considered in PHVMD19, but in HVMD19 the DGIA effect is specifically accounted for prior to the NIE determination and removal. With reference to Table 4 and Figs. 4a, b and 8 one can see that the remaing DGIA effect in PHVMD19 (which is removed in HVMD19) is as much as − 1.4 km in the Gulf of Bothnia, while the large positive effect in NW Norwegian Sea has some other origin. As already noted in Sect. 5.3, applying the specific DGIA effect also reduces the difference between the gravimetric–isostatic determined Moho depth and the seismic ones. This implies that a specific correction for the DGIA yields a small but significant in the Moho depth in the center of the uplift region.

Table 4 Statistics of the remaining DGIA effect in PHVMD19, removed in HVMD19
Full size table
Fig. 8
figure 8

The remaining DGIA effect in PHVD19. (Unit: km)

Full size image

6 Conclusions

In this study we determined a new gravimetric–isostatic Moho depth model for Fennoscandia using the VMM isostatic hypothesis from the refined Bouguer gravity disturbance and Earth2014 topographic data over 1° × 1° blocks on land and ocean. To that end, the refined Bouguer gravity disturbance was primarily reduced for gravity of topography, density heterogeneities related to bathymetry, ice, sediments, and other crustal components using stripping gravity corrections, and it was further corrected for NIEs of nuisance gravity signals from mass anomalies below the crust due to crustal thickening/thinning, thermal expansion of the mantle, Delayed Glacial Isostatic Adjustment (DGIA) and plate flexure. As Fennoscandia is a important area for GIA investigation, we specifically studied the DGIA effects on gravity and Moho depth to figure out if we can improve the stripping of the GIA related gravity signal from the observed gravity field by a specific correction prior to applying the general NIE. To fulfill this, we used the spectral HW 10–23 of the gravity field to calculate DGIA effect, as it yields the maximum correlation, of the order of 0.92, with the newest land uplift model NKG2016LU of the Nordic region. Also we validated the Moho depths estimated from the PHVMD19 and HVMD19 with the seismic based CRUST19 and independent GRAD09 model, showing that the RMS difference of HVMD19/PHVMD19 and the seismic CRUST19 and GRAD09 models are 3.8/4.2 km and 3.7/4.0 km when the above strategy for removing the DGIA effect is/is not applied, respectively, and the mean value differences are 1.2/1.4 km and 0.98/1.4 km, respectively. Hence, the specific correction for the DGIA effect slightly improves the agreement of the gravimetric–isostatic model and seismic models to about 10%.

Finally, we believe that the significant differences between gravimetric and seismic Moho geometries along the Norwegian coastline may be caused by remaining non-isostatic effects in the gravimetric models. However, this issue needs further investigation.