Interpretation of geoid anomalies in the contact zone between the East European Craton and the Palaeozoic Platform—II: Modelling of density in the lithospheric mantle

Summary

We present analysis of lateral variations of density in the upper mantle in the area of contact of the precambrian East European Craton (EEC) and the Palaeozoic Platform (PP) in Poland, obtained by analysis of the gravimetric geoid undulations. A precise 3-D density model of the crust in the study area down to a depth of 50 km, discussed in the first part of this paper (Majdański et al., in press) did not explain all features of the observed geoid. This suggests that these features can be due to density inhomogeneities in the upper mantle. To estimate them, we performed inversion of a residual between the observed geoid and undulations caused by the 3-D density distribution in the crust. Basing on the assumption of local isostatic compensation and Pratt-Hayford isostasy model, the density distribution in the upper mantle was parametrized as a 40-km-thick layer located above the assumed compensation depth of 140 km and subdivided into irregular blocks. The boundaries of the blocks were defined according to boundaries of major tectonic units in the study area and position and shape of the most pronounced anomalies in the residual geoid. A series of sensitivity tests calculated for such density heterogeneities in the upper mantle showed that they can produce geoid undulations of the order of several metres. The density values in each unit were taken as model parameters for the inversion procedure, and inverse problem was solved using global optimization with constraints. The density variations in the upper mantle in the final model correlate well with the surface heat flow. This suggests that these variations can be due to diversity in mantle temperature. The Trans-European Suture Zone (TESZ), which is a major suture separating the EEC from the PP, is not observed as a distinct unit in the mantle. Instead, our study suggests continuation of the lithosphere of the EEC beneath the PP and confirms subdivision of the TESZ into terranes with distinctly different evolution. A simple calculation shows that condition of local isostatic equilibrium is almost achieved at the assumed compensation depth, with the exception of several small-scale units. However, different loads and mass deficiencies in the lithosphere are compensated at different depth levels. This proves that the lithosphere is alternation of mechanically strong and weak layers, both in the crust and upper mantle, and that different loads can be compensated at different depth. However, the total compensation of all lithospheric loads occurs in the upper mantle.

Inverse theory, Gravity anomalies and Earth structure, Europe

1 Introduction

In geodynamic studies, it is important to have knowledge about distribution of density and viscosity in the Earth's mantle. Although Bouguer gravity anomalies have historically been used to constrain these parameters in the crust (e.g. Thompson & Zoback 1979; Hoffmann et al. 2003; Kozlovskaya et al. 2004), geoid anomalies have been considered as sources of information about deeper lithospheric structure (e.g. Doin et al. 1996; Kaban et al. 1999; Ebbing et al. 2006; Coblentz et al. 2007). On a local scale (hundreds of kilometres), geoid anomalies are primarily caused by changes in surface elevation and density distribution in the lithosphere; that is why they can be used, also to study mechanisms of compensation of topographic undulations by density variations in the Earth (isostasy). Initially, it was considered that the local isostatic compensation occurs at the crust-mantle boundary (Airy-Heiskanen isostasy model). However, analysis of the compensation depth of topographic features presented by Martinec (1994a,b) demonstrated that the compensation occurs at a depth of 100–150 km, which better agrees with the Pratt-Hayford model of isostasy. The result obtained by Martinec (1994a,b) suggests that the whole lithosphere (not just the crust) is involved into the compensation.

In our paper, we present analysis of the density distribution and isostatic compensation in the lithosphere for the area adjoining the Trans-European Suture Zone (TESZ), separating the Precambrian East European Craton (EEC) from younger tectonic units of western Europe (Fig. 1). The study area is a region of large contrasts in the crustal structure. As was demonstrated in the first part of our study (Majdański et al., 2009), large variations of crustal thickness and density do not explain all the features of the observed geoid in this area. As the area appears to be in a generally ‘stable’ state (Ziegler et al. 2006), the local isostatic compensation is most probably achieved by variations of density in the upper mantle.

Figure 1.

Geological setting of the main tectonic units of the Trans-European Suture Zone (TESZ) in the Central Europe. Yellow square shows the location of the study area. The geographical location is presented in the corner.

The structure of the upper mantle in the area has not been studied previously by high-resolution passive seismic experiments; that is why precise location of density heterogeneities in the upper mantle cannot be constrained by seismic data, like it was done recently by Ebbing et al. (2006). However, global seismic studies showed that such heterogeneities exist. Namely, velocities of the seismic body waves are higher in the upper mantle of the EEC down to the depths of 250–300 km compared with the values from global models. Models of the upper-mantle velocities (Shapiro & Ritzwoller 2002; Funke et al. 2003; Eaton & Jones 2006) do not indicate any low velocity layer down to 300 km beneath the EEC. In contrast, the typical Variscan region (Fig. 1) is characterized by lithosphere thickness of 90–120 km. Beneath the Caledonides area and the Bohemian Massif (Fig. 1), the thickness of the seismic lithosphere reaches 90–130 km and 120–140 km, respectively, whereas the Carpathians are associated with the lithosphere thickness of 130–150 km (Panza et al. 1980; Du et al. 1998; Bijwaard & Spakman 2000; Plomerova et al. 2002). Zielhuis & Nolet (1994) showed that the TESZ separates regions of relatively high S-wave velocities beneath the EEC from the area of low velocities under the Palaeozoic Platform (PP). Observed difference between S-wave velocities varies from 6–7 per cent at a depth of 80 km to 3–4 per cent at a depth of 140 km.

Geothermal studies indicate that the TESZ is an area of pronounced change in the heat flow, separating the area of the ‘cold’ lithosphere of the EEC with heat flow lower than 40 mW m⁻² from ‘warm’ lithosphere of the Palaeozoic terranes with heat flow in range of 40–70 mW m⁻² (Majorowicz & Plewa 1979). Majorowicz et al. (2003) demonstrated that heat flow variations in the area cannot be explained by changes in radiogenic heat production in the crust only and suggested that significant part of these variations is due to diversity in mantle heat flow. The thickness of the thermal lithosphere changes from 170–230 km in the EEC to 70–80 km in Palaeozoic terranes (e.g. Artemieva 2003; Artemieva et al. 2006; Majorowicz et al. 2003).

A large-scale 3-D gravity modelling (Yegorova & Starostenko 2002; Kaban et al. 2003; Yegorova et al. 2007) suggests that densities in the mantle are different to the both sides from the TESZ. This agrees with the analysis of density distribution in the crust and the upper mantle for the territory of Poland presented by Grabowska et al. (1998). In addition, 2-D gravity modelling along the LT-7, TTZ and POLONAISE'97 wide-angle reflection and refraction profiles (Krysiński et al. 2000) indicates lateral density variations in the mantle beneath these profiles. Thus, there is a lot of previous studies suggesting laterally heterogeneous density distribution in the upper mantle of our study area. These inhomogeneities can be reflected also in the geoid shape.

2 The Total Effect of the Crust on Geoid Undulations

In the first part of this paper (Majdański et al., in press), we estimated the total effect of density inhomogeneities in the crust on the geoid. The 3-D model of the crust of the study area was compiled from a number of 2-D velocity models along deep seismic sounding profiles. The model was divided into several major layers: topography, sediments, crystalline upper crust, lower crust and the upper mantle and extended to the maximum depth of 50 km. The gravity effect of each layer was calculated separately. To assure the best possible fit to the observed geoid (Fig. 2a), the crystalline crust and the upper mantle were subdivided into major tectonic units (see Majdański et al. 2009 for details). The optimal fit of the calculated and observed geoid (Fig. 2c) was achieved for the division presented in Fig. 2(b). The residual presented in Fig. 2(d) shows a number of positive and negative anomalies with amplitude of about 6 m. These anomalies could not be explained by variations in crustal thickness and density only.

Figure 2.

Comparison of the observed geoid (a) with the effect of the 3-D crustal density model to geoid height (c). Subplot (b) shows major tectonic units recognized in the crust. Residual geoid height (d) after eliminating the effects of the topography and the crust shows the large-scale effects that could be interpreted as caused by density inhomogeneities in the upper mantle. PO, Pomerania unit; KU, Kuiavian unit; MP, Malopolska unit.

As the lithosphere in our ‘stable’ continental study area appears to be in a state of local isostatic equilibrium, the geoid anomalies in each point of the study area are proportional to the dipole moment of the density distribution above the compensation depth (Turcotte & Schubert 2002). This moment is defined as

1

where h is a total height of the vertical column over the compensation depth and Δρ(x, y, z) is anomalous density with respect to selected reference density column. From this it follows that the density anomalies in the mantle averaged over the column above the assumed compensation depth can be estimated directly from the observed geoid if the density structure in the crust is known. In this paper, we estimate the lateral variations of average mantle densities in our study area, using inversion of the residual between the observed geoid and calculated effect of the 3-D crustal model.

3 Inversion Procedure

3.1 A priori information: constraints on density distribution in the mantle

To constrain the inversion, we estimated possible variations of density that can be expected in the upper mantle of our study area. Physical properties of continental lithospheric mantle are relatively well known not only from geophysical studies but also from direct measurements, using mantle rocks that have been overthrusted or exhumed to the surface by various tectonic and magmatic processes. These studies have shown that the upper mantle is composed mainly of peridotites (cf. Nixon 1987; Griffin et al. 1998; Downes 1997). The main rock forming minerals for them are olivine, orthopyroxene, clinopyroxene and spinel. The olivine content in the lithospheric mantle can vary from 50 per cent to nearly 100 per cent (Downes 1997).

Another important type of rocks in the continental upper mantle is represented by eclogites. Generally, they were formed in the past as a result of tectonic processes involving subduction of the crust into the mantle. These rocks may be present in the upper mantle in a form of remnants of subducted slabs. Generally, mantle eclogites may have densities up to 0.2 g cm⁻³ higher than those of peridotites. However, estimates of amount of eclogites in the upper mantle are controversial. Schulze (1989) proposed that their amount in the upper mantle is insignificant (several per cents), as more dense rocks would sink into the deep mantle and become recycled. Recent studies of upper-mantle xenolithes from Precambrian areas show, however, that proportion of high-density eclogites in kimberlite pipes may be significantly larger than several per cents (Kopylova et al. 2004).

The lithospheric mantle beneath continents has undergone significant (probably multiple) melting through geological time and is depleted in such components as Fe, Al, Ca and Ti (cf.Griffin et al. 2003; Lee 2006). The conversion from fertile to depleted mantle material is expressed by a decrease in clinopyroxene and orthopyroxene and relative increase in olivine and also by an increase of Mg content. Griffin et al. (2003) suggested that the mean geochemical composition and density of the subcontinental lithospheric mantle beneath various tectonic units generally depend on the tectonothermal age—younger lithospheric mantle is less depleted and hence is denser. Thus, the lithospheric mantle beneath the young (Phanerozoic Western Europe) has the average density of 3.37 g cm⁻³ at a depth of 70 km, whereas the older (Proterozoic) mantle is less dense (3.34 g cm⁻³; Gaul et al. 2000; Griffin et al. 2003).

Compared with the lithospheric mantle, the asthenosphere is fertile and chemically more homogeneous due to convection. The main component of a partially molten material in the asthenosphere is tholeiitic basalt that is produced by melting of upper-mantle peridotite (Herzberg 1995). The density of the asthenosphere depends upon the density and the content of partially molten material that can be estimated from seismic velocities and seismic attenuation obtained by teleseismic tomography (Sobolev et al. 1997; Sato et al. 1998; Petit et al. 2002). In accordance with these studies, the content of melt in the asthenophere is rather small (less than 5 per cent).

The other constraint on the amount of partial melt in the asthenosphere is the stability condition for the olivine-molten basalt mixtures. As the density of the molten basalt is less than that of the solid olivine matrix, the olivine-molten basalt system is mechanically stable only in the case when the melt is concentrated within isolated inclusions (pockets). If the inclusions form the interconnected network, the flow of melt is initiated, and the melt tends to concentrate in large volumes, like it is observed beneath the mid-ocean ridges and volcanoes. Theoretical modelling of elastic and electrical properties of the olivine-molten basalt system (Kozlovskaya & Hjelt 2000) demonstrated that in some cases, only 5 per cent of melt inclusions is enough to form a perfectly interconnected network. Thus, assuming the melt content less than 5 per cent and the density of the molten basalt under upper-mantle pressure and temperature conditions as 2.72 g cm⁻³ (Manghnani et al. 1986), the lithosphere-asthenosphere density contrast can be estimated as being less than 0.04 g cm⁻³.

The density contrast between the asthenosphere and mantle lithosphere beneath the TESZ can be roughly estimated from the absolute values of S-wave velocity (V_S) in the northern part of the TESZ presented by Cotte et al. (2002). They showed that the V_S in the mantle lithosphere in the TESZ is 4.57 km s⁻¹, whereas the V_S in the asthenosphere is 4.36 km s⁻¹. Using the scaling factor relating decrease in V_S to the decrease in density, which is equal to 0.05 at a depth of about 100 km (Deschamps et al. 2001), the asthenophere density beneath the northern part of TESZ is 3.34 g cm⁻³.

The estimates of the density contrast at the lithosphere-asthenosphere boundary can be obtained also from regional (medium-wavelength and long-wavelength) gravity studies in Europe. Marquart & Lelgemann (1992) used the density values of 3.23–3.32 g cm⁻³ for the lithospheric mantle and 3.2 g cm⁻³ for the asthenosphere in the study of medium-wavelength geoid anomalies over Europe along the European Geotraverse. Yegorova & Starostenko (2002) assumed that the density contrast between the European lithospheric mantle and the asthenosphere is 0.04–0.05 g cm⁻³. Cella et al. (1998) obtained the densities in the range of 3.27–3.33 g cm⁻³ at a depth of 60–80 km for the upwelling asthenosphere beneath the Central Mediterranean. Summarizing all the estimates of the density in the asthenosphere cited above, we can conclude that the range of possible values of the asthenosphere density for Europe is 3.2–3.34 g cm⁻³.

This abundant information about possible density variations in the upper mantle from previous studies was used as a priori information in the inversion problem.

3.2 Model parametrization

In the inversion problem, we fixed the density distribution in the crust down to the Moho boundary and inverted for density distribution in a mantle averaged over the depth range from the Moho boundary to the assumed compensation depth.

According to classic flexural isostasy hypothesis, response of the lithosphere to loading can be modelled as that of mechanically strong lithosphere (MSL) overlying a mechanically weak asthenosphere (Watts 1978, 2001; Turcotte & Schubert 2002). Under such assumption, the depth of compensation (that is, the depth where the lithosphere is isostatically adjusted) coincides with the lithosphere-asthenosphere boundary (LAB). At present, this concept is more sophisticated, because several definitions of lithosphere exist (e.g. seismic lithosphere, electrical lithosphere, thermal lithosphere, rheological lithosphere, see Martinez & Wolf 2005, for review). In young continental lithosphere, the LAB is associated with a layer of low seismic velocities due to presence of partially molten material (seismic LAB) and can be estimated directly from results of passive seismic experiments. For the western part of the TESZ, the depth to the LAB was estimated by Cotte et al. (2002) and Plomerova et al. (2002). They interpreted teleseismic data of the TOR passive experiment and showed that the depth to the LAB varies between 50 and 140 km beneath the TOR array.

Unfortunately, the seismically defined LAB was not available for our study area, as the TOR array was located outside of it. In addition, the TOR experiment showed that the low velocity layer (seismic LAB) disappears to the northeast from the TESZ beneath the Precambrian EEC (Cotte et al. 2002; Shomali et al. 2006), where the seismic lithosphere is very thick (Shapiro & Ritzwoller 2002; Funke et al. 2003; Eaton & Jones 2006; Plomerova et al., in press). In such areas, the asthenosphere is not associated with the partially molten material, but it is defined as a convecting layer where thermal gradient is adiabatic (thermal LAB; McKenzie & Bickle 1988; O'Reily & Griffin 1996; Lee 2006). However, the thickness of mechanically strong part of the lithosphere (rheological lithosphere or elastic lithosphere) there is generally lower than the depth to the thermal LAB. For example, Milne et al. (2001) obtained a glacial isostatic adjustment model for the Fennoscandia with a 120 km thick strong layer overlying a viscous upper mantle. This estimate agrees well with result of Kaikkonen et al. (2000), who determined that the thickness of the MSL varies between 70 and 140 km in the central part of the Fennoscandia.

Taking these results of previous estimates of lithosphere thickness into consideration, we can assume that isostatic adjustment of the lithosphere in our study area is achieved at a depth of no more than 140 km, either at a seismic LAB (southwest from TESZ) or at the boundary of MSL (northeast from the TESZ). However, the ‘true’ compensation depth associated either with the LAB or with the MSL is unknown for our study area.

As was shown by Kaban et al. (1999), any compensating mass distribution can be approximated by a thin layer placed above the effective compensation depth. Therefore, in our study, we approximated the real distribution of density in the upper mantle by an inhomogeneous layer of fixed thickness of 40 km placed at the compensation depth of 140 km. This depth is in a good agreement with the maximum compensation depth inferred by Martinec (1994b). In addition, it is lower than the depth to the seismic LAB in the TESZ and southwest from it (Cotte et al. 2002), and it is lower than the thickness of the MSL in Precambrian part of Europe (Kaikkonen et al. 2000; Milne et al. 2001).

Following the idea of multigrid parametrization of inverse problems (Kissling et al. 2001), we parametrized the density in the model by two 3-D grids. A regular fine-scale grid (0.5° × 0.333° × 1 km) was defined to calculate the forward problem (forward problem grid). In addition, the study area was extended to infinity to minimize the boundary effects and take into account the unknown density distribution outside the study area (see Majdanski et al. 2009 for details) in forward problem calculations. The second irregular grid (inversion grid) was defined over the forward problem grid. For this, the model was extended vertically to 140 km in depth. The uppermost part of the model was the layer from the top of the model down to the depth of 60 km, representing the density distribution in the crust. To model the density in the mantle, a 40 km thick and laterally heterogeneous layer was placed above the depth of 140 km (Fig. 3). This configuration gives the highest influence (see Fig. 4) on calculated anomalies; that is why we didn't use other layers. The layer was subdivided into blocks representing individual tectonic units of the lithosphere. To improve the stability of solutions and to obtain realistic values of densities, the densities for the individual blocks were defined with respect to the stable and cold EEC, for which the constant density value of 3.35 g cm⁻³ was assumed. These relative density values with respect to the constant value of 3.35 g cm⁻³ formed a vector of model parameters m for the inversion problem. Solution to the inverse problem is the distribution of average density values in individual blocks.

Figure 3.

Schematic diagram of two levels of parametrization of the model. The fine grid marked with thin lines with block size of 0.5° × 0.333° × 1 km were used for forward calculations, whereas division of the layer in the mantle into units was used in the inversion. The crust structure (green) was fixed. The grey colour indicates position of the equivalent layer in the upper mantle above the assumed compensation depth of 140 km.

Figure 4.

Sensitivity test for three different units: Malopolska Unit (a), small unit in the TESZ (b) and the EEC (c) showing the effect of the relative density contrast of +0.05 g cm⁻³ in selected unit at different depth: 0–60 (without the crust), 60–100, 100–140, 140–180 and 180–220 km. Colour scales range is constant for the each unit.

3.3 Sensitivity of the geoid to density variations in the upper mantle

To investigate the sensitivity of the geoid to density variations in the mantle, we calculated the effect of relative density anomalies of order of 0.05 g cm⁻³ located at different depth beneath tectonic units of different size. As demonstrated in Section 3.1, the residual density of 0.05 g cm⁻³ corresponds to the density inhomogeneity due to presence of asthenosphere (low-density of 3.3 g cm⁻³) or eclogitic rocks with high content of garnet (high density of 3.4 g cm⁻³). The results of sensitivity analysis are presented in Fig. 4.

The analysis demonstrates that density anomalies of such a scale can produce geoid undulations of an order of several metres. Therefore, the anomalies of geoid height seen in the residual geoid (Fig. 2d) can be modelled using the selected model parametrization.

3.4 Inversion problem formulation

The vector m can be obtained by inverting residual geoid anomaly that form the observation vector d_obs = (d₁, d₂, …d_Q), where d_i, i = 1, …, Q are values of the residual in the points of the selected grid on the observation surface. The inversion procedure maximizes the a posteriori conditional probability density of the vector m after the experiment (Tarantola 1987; Menke 1989):

2

where p(m) is the a priori probability density function (PDF) of the vector m,M is a set of feasible solutions, p(d_obs ǀ m) is the conditional probability density of the experimental data and p(d_obs) is the marginal probability density of the experimental data. As p(d_obs) does not depend on m, it can be assumed to be equal to some constant b. Then the maximum of the conditional PDF (2) gives an estimate for the parameter vector m (maximum a posteriori solution). The PDF p(d_obs ǀ m) depends on an operator g that allows us to calculate the data values from known values of model parameters:

3

Following the formalism proposed by Kozlovskaya (2000), we used a possibility functions Pos(m) instead of the a priori PDFs to represent the a priori information about possible range of densities in the upper mantle. An example of possibility function for density variation of ith individual block is

4

where m_upper and m_lower, are upper and lower constraint on model parameters. Individual possibility functions are combined into multidimensional possibility function for parameter vector m using fuzzy logic rules (Kozlovskaya 2000).

In our study, we assume that the a posteriori conditional probability density of the vector m is generalized Gaussian of order 2, with diagonal covariance matrices. In this case instead of maximization of (2), it is possible to use minimization of the following function:

5

4 Inversion Results

A serious problem for our study was absence of detailed seismic model of the upper mantle. Because of it, we performed inversion of the residual geoid for several trial subdivisions of the equivalent layer in the upper mantle into blocks. Examples of them are shown in Fig. 5 together with the corresponding calculated geoid and residual between the latter and observed geoid. In each model the boundaries of separate lithosphere blocks were defined taking into account seismic, tectonic and thermal data, and the fit to the observed geoid was analysed.

Figure 5.

Examples of the calculated geoid (left-hand panel) and residuals (middle panel) for different inversion grids with the final density values (right-hand panel). The boundaries of each blocks were defined by boundaries of major tectonic units (division a), based on the velocity structure in the crust (b), including information about the heat flow distribution (c) and final model including all information (d). EEC, East European Craton; PU, Pomerania unit; KU, Kuiavian unit; MU, Malopolska unit; PP, Palaeozoic Platform; SM, Sudetes Mountains; BM, Bohemian Massif; CPT, Carpathians.

Our initial model (Fig. 5a) has been constructed under assumption that the structure of the mantle beneath our study area is relatively simple, as has been discovered by global seismic studies mentioned in Chapter 1. Such structure can be represented by four basic blocks: PP; Teisserye-Tornquist Zone (TTZ); Carpathians (Ca) and EEC. For such a simple model, we were not able to achieve a good fit to the observed geoid, however. It is seen (Fig. 5a) that the residual between observed and calculated geoid varies between −6 and 6 m, which is of the same order as the residual between the observed geoid and the effect of the crust. This suggested that the upper mantle in our study area is more differentiated laterally than can be concluded from global seismic studies.

To improve the fit to the observed geoid, we modified the initial model and included several additional blocks into the mantle (Fig. 5b). The boundaries of the blocks were defined taking into account mainly the structure of the crust and the uppermost mantle from seismic data. For example, it was necessary to introduce an additional unit (denoted by a question mark) to explain the increase of the amplitude of the residual to 4 m in the SE part of the TESZ. Position of the unit and its spatial extend were defined in accordance with the depth to the Moho boundary, which has a pronounced depression in this place. On the basis of the differences in the P-wave velocity model (Janik et al. 2002; Dadlez et al. 2005), three minor lithosphere blocks were separated in the TESZ zone. Namely, division of the crust into Kuiavian, Pomeranian and Małopolska blocks was extended into the subcrustal lithosphere. Taking into account the thickness of the sedimentary layer (Majdański et al. 2006), the Bohemian Massif was separated as a distinct unit in the PP, as a region characterized by generally thin sediments. Fig. 5(b) shows that this kind of model parametrization is better than the previous one, because the absolute value of the residual between the observed and calculated geoid has decreased significantly (down to 3 m), except of marginal regions of the model.

To further improve the fit to the observed geoid, we added several minor tectonic units to the model. The units were distinguished basing on tectonic, seismic and geothermal data (Fig. 5c). During inversion, the boundaries of separate units were modified to obtain the best fit to the observed geoid. The final density model and its gravitational effect are shown in Fig. 5(d). The calculated geoid has the shape similar to the observed, one and the absolute values of the residual are less than 2 m. Further increase of the number of blocks and modification of block boundaries did not result in improvement of the fit to the observed geoid. The size and shape of the remaining residual anomalies (Fig. 5d) suggests that they can be artefacts due to uncertainties in the model of the crust discussed in the first part of our paper (Majdanski et al. 2009). Therefore, further improvement of the crustal model with the use of the Bouguer anomaly would result in more consistent model of the mantle.

5 Resolution Analysis

Table 1 represents correspondence between model parameters for the final model; it means between components of the vector m and geological units marked in Fig. 5(d). The quality of the final solution is illustrated by Figs 6 and 7. Fig. 6 shows coefficients of correlation between separate model parameters. It is seen that correlation between them is generally weak, as the absolute values of correlation coefficients are less than 0.4. Such a weak correlation can be neglected in resolution analysis (Sambridge 1999a,b), and resolution for each model parameter can be characterized by projection of multidimensional a posteriori PDF to axes corresponding to separate model parameters (Fig. 7). As seen, the a posteriori PDF is non-Gaussian.

Table 1.

Average density values in the upper mantle obtained by inversion of the residual geoid with assumption of 3.35 g cm⁻³ density for the EEC.

Figure 6.

Coefficients of correlation between separate model parameters for the final density model shown in panel d. The values of coefficients are shown by colour scale. The correspondence between the model parameters and units in the model is presented in Table 1.

Figure 7.

Multidimensional a posteriori PDF (5), illustrated by its projections to axes corresponding to separate model parameters.

Fig. 7 shows that parameters 1 (western part of the Pomeranian unit) and 6 (Unit 1 inside the EEC) are not well resolved within the range of density variations assumed in our study. The density in the west of the Pomeranian unit tends to be lower than the minimum density for this block assumed in our study. As this unit is located on the margin of our model, it can indicate that the area with low density is actually larger. Most probably, it continues to the NW, outside the study area. On the contrary, the density of the Unit 1 tends to be significantly higher than the maximum density corresponding to eclogitic rocks. As this unit is located in the middle of our study area, this can indicate that either the thickness of this block is larger than 40 km or the crustal density above this block was underestimated.

The well-resolved parameters are 2 (Kuiavian unit), 3 (Malopolska unit), 4 (Palaeozoic Platform), 5 (Carpathians), 7 (Bohemian Massif), 8 (Sudetes), 9 (east of the Pomeranian unit) and 10 (Unit 2 in the TTZ).

6 Discussion

6.1 Origin of lateral density variations in the upper mantle

Fig. 5 (right panel) shows the density values in the equivalent layer for different number of units used in the inversion and the density values for the final model. Although these values are confined to the assumed compensation depth, they may be treated as a first-order approximation of lateral density variation in the upper mantle of our study area.

The result for the inversion with simplified structure of the upper mantle (4 units) shows very low densities of 3.26 g cm⁻³ in the TESZ area and relatively low densities of 3.29 g cm⁻³ for the PP. However, this density distribution in the mantle did not improve the fit to the observed geoid. This suggests that the sharp boundary between the EEC and TESZ seen in the crust (Majdański et al., in press) does not continue to the mantle.

The models with finer division of the upper mantle (Figs 5b–d) fit better to the observed geoid and demonstrate good qualitative agreement with the heat flow. Moreover, all models show the same trend: namely, high density beneath the southwestern part of the TESZ (Kuiavian unit) and low density in its northwestern part (western part of the Pomeranian unit). This means that this feature is stable and independent of details of model parametrization. Our final model with the best fit to the observed geoid (Fig. 5d) shows that the density in the mantle beneath Kuiavian and Malopolska Units does not differ from that of the EEC. This suggests that the upper mantle of the EEC continues beneath the SE part of the TESZ and is cold and mechanically strong.

The high-density values beneath the Bohemian Massif and Unit 1 in the EEC are also invariant with respect to model parametrization. One explanation for these high values is presence of eclogitic rocks posed into the upper mantle during subduction processes in the past. However, the high density in the upper mantle of the Unit 1 may also be an artefact resulting from underestimated crustal density in this unit. In addition, this parameter is not well resolved (Fig. 7).

Lower density in the upper mantle beneath the PP and Carpathians are two other features independent on model parametrization. They can be explained by presence of partially molten asthenosphere (cf.Babuska & Plomerova 1992; Ladenberger et al. 2006). Our final model suggests, however, that the structure of the upper mantle beneath the PP is heterogeneous, because the low density appears only beneath Sudetes. On the contrary, the density in the northern part of the Palaeozoic Platform (denoted as the PP in Fig. 5d) is high and approaches that of the EEC. Grad et al. (2003) demonstrated that the uppermost mantle beneath this unit has very high P-wave velocities of 8.4 km s⁻¹. Such a high values of velocity can correspond to eclogitic rocks with the density of about 3.4–3.5 g cm⁻³ (Kobussen et al. 2006). However, the density value for the Palaeozoic Platform obtained in our study (Table 1) is lower than the density of eclogites. The reason for this apparent disagreement between high P-wave velocities and relatively low densities may be presence of low-density asthenosphere beneath the high-density upper mantle lid. As we estimated the density in the upper mantle averaged down to the assumed compensation depth, such variations of density with depth compensate each other in our model and cannot be discriminated. The presence of the asthenosphere below the Palaeozoic platform agrees well with the high heat flow in this area.

The slightly low density value in the East Pomeranian unit (Fig. 5) may be an artefact, because this unit is spatially coincident with the area of the largest vertical subsidence in Poland caused by mining activity (Kowalczyk 2006). Due to these movements, the condition of local isostatic compensation is violated there, and the correspondent geoid anomaly is projected into the upper mantle.

6.2 Analysis of local isostatic compensation

As our density model was obtained under the assumption of local isostatic compensation, we checked whether this condition is satisfied in our final density model. The calculations were performed in similar way as in publication by Ebbing et al. (2006). For each postulated depth of compensation, an average density was calculated every 1 km in a column with respect to a mean density at given depth of the whole model. Such an estimate is proportional to the dipole moment of the density distribution above the compensation depth (eq. 1), but it is more convenient for visualization.

The results presented in Fig. 8 shows that condition of local isostasy is almost fulfilled in our study area, except for three units that still can be distinguished in a plot corresponding to 150 km (Bohemian Massif, Pomerania unit and Unit 1 in EEC). However, comparison of plots in Fig. 8 to the relief of major interfaces in the crustal model (Fig. 9) suggests that compensation of various loads (or mass deficiencies) occurs at different depth. The load due to topography is not compensated in the uppermost 25 km of the crust, because areas with high topography (Sudetes, Bohemian Massif and Carpathians) are clearly visible in plots corresponding to 10 and 25 km. However, these features almost disappear in the lower crust at the depth of 35 km. This depth is shallower than the Moho boundary beneath these units (Fig. 9d), but it is lower than the depth to the lower crust (Fig. 9c). This indicates that compensation of topographic loads in these units occurs in the lower crust (but not at the crust-mantle boundary).

Figure 8.

Estimation of the isostatic compensation calculated as an average density residual in a vertical column of the 3-D density model relative to the average density value. Calculations were performed under assumption of the compensation depth at 10, 25, 35, 50, 100 and 150 km, respectively.

Figure 9.

Major boundaries in the 3-D crustal model: (a) shows the topography map according to NOAA (1988); (b) is the thickness of the sedimentary cover; (c) shows the depth to the lower crust in the 3-D seismic velocity model (Majdanski et al. 2009) defined as a layer with velocities of 6.6 km s⁻¹ < V_P < 7.9 km s⁻¹; (d) shows the Moho boundary. Numbers in (a) according to Table 1.

At the depth of 25 km, one can also distinguish non-compensated mass deficiency in the crust due to thick sediments in the TESZ. However, the effect of sediments in the NW part of the TESZ (Pomeranian Unit) disappears almost completely at a depth of 35 km, whereas the effect of sediments in the SW part of the TESZ (Kuiavian and Malopolska Units) is still seen at a depth of 35 km. This indicates that compensation of thick sediments in the TESZ occurs in the lower crust, as the depth to the Moho boundary is 35–38 and 40–45 km in the NW and SE part of the TESZ, respectively (Figs 9c and d).

The major feature in the crust in our study area is the strong contrast in density and thickness of the crust between the EEC and younger units to the west from the TESZ. This feature is not compensated at the Moho boundary because the difference between these areas is clearly visible in the plots corresponding to 50 and 100 km compensation depth. However, this difference disappears almost completely below the assumed compensation depth (150 km), where the effect of inhomogeneities in the upper mantle was taken into account (Fig. 8).

The diversity in the depth of compensation of various loads in tectonic units of different age indicates also diversity in rheological structure of the lithosphere in our study area. Recently, Thatcher & Pollitz (2008) suggested that isostatic adjustment in continental lithosphere at long timescales (∼10⁶–10⁷ yr) occurs with most lithospheric stress supported by an upper crust overlying a much weaker ductile substrate that includes both the lower crust and upper mantle. Our study does not support this uniform mechanism of isostatic adjustment of the lithosphere. Instead, it shows that different loads and mass deficiencies are compensated at different levels of the lithosphere. The loads due to topography in Sudetes, Bohemian Massif and Carpathians as well as thick sediments in the TESZ are compensated at a depth corresponding to the lower crust, which suggests that the lower crust is mechanically weak there. However, compensation of variations in crustal thickness and structure requires the mechanically weak layer in the mantle also, although we cannot define the true compensation depth in our study. Generally, our result agrees well with the classical ‘jelly sandwich’ distribution of strength in the lithosphere (Ranalli 1995; Burov & Watts 2006) that implies that the lithosphere is alternation of mechanically strong and weak layers, both in the crust and upper mantle, and that different loads can be compensated at different depth.

7 Conclusions

The main conclusions from our study can be summarized as follows:

• (1)

  Variations of thickness and density of the crust in the area adjoining the TESZ can explain the major part of the observed geoid undulations but not all of them. The residual between the observed geoid and undulations caused by density heterogeneities in the crust is of an order of several metres. It can be explained by density heterogeneities in the upper mantle.

• (2)

  The condition of local isostatic compensation is almost fulfilled in our study area at the assumed compensation depth. However, compensation of different loads and mass deficiencies occurs at various depths in different tectonic units, both in the crust and upper mantle.

• (3)

  Compensation of topographic loads in Sudetes, Bohemian Massif and Carpathians, as well as thick sediments in the TESZ, occurs in the mechanically weak ductile lower crust.

• (4)

  Strong contrast in density and thickness of the crust between the EEC and younger units is not compensated at the Moho boundary. The total compensation of all loads in the lithosphere occurs in the upper mantle.

• (5)

  The low densities in the upper mantle below the Sudetes, Carpathians and Pomeranian Unit correlate well with the regions of increased heat flow. They can be interpreted as indicating presence of partially molten material (asthenosphere).

• (6)

  The upper mantle of the EEC extends beneath the SE part of the TESZ up to the Carpathians front. It is indicated by densities in the upper mantle of the Kuiavian and Malopolska units that are similar to the density of the EEC. This implies different tectonothermal evolution for SE and NW parts of the TESZ and suggests that the upper mantle there is cold and mechanically strong.

• (7)

  Relatively high densities observed in the northern part of the PP do not correlate with the high heat flow in this unit, but correlate well with the high P-wave velocities in the uppermost mantle. One possible explanation of this phenomenon can be presence of high-density eclogitic rocks in the uppermost mantle underlain by material with low densities (probably, asthenosphere).

• (8)

  Deviations from the state of local isostatic equilibrium and high density in the mantle observed in the Bohemian Massif may suggest that these units are supported by mechanical strength of the underlying lithosphere.

Acknowledgments

We are very grateful to Dr J. Ebbing and the second anonymous reviewer, whose constructive comments helped to improve the earlier version of our paper. This work was partially supported by MNISW grant No 2P04D05428. MM and MŚ's stay at Sodankyla Geophysical Observatory of Oulu University was financed by the Academy of Finland (grants No. 112603 and 122645).