Interpretation of geoid anomalies in the contact zone between the East European Craton and the Palaeozoic Platform—I. Estimation of effects of density inhomogeneities in the crust on geoid undulations

Summary

In our paper, we compile a detailed 3-D velocity model of the crust for the territory of Poland and investigate the effect of a 3-D density distribution in the crust on geoid undulations. The major tectonic feature in the study area is the Trans European Suture Zone (TESZ) separating the Precambrian East European Craton from younger tectonic units of Palaeozoic Europe. High-resolution seismic models of the crust and upper mantle from recent controlled-source experiments were interpolated into a 3-D P-wave velocity model for the whole Poland. The model was converted to a density model assuming different equations connecting seismic velocities to densities for crustal rocks, sediments and the uppermost mantle. The crustal density model was used to calculate synthetic geoid undulations for the territory of Poland and to estimate influence of separate layers (topography, sediments, crust and uppermost mantle) on the geoid. To calculate the gravity effect of a large-scale 3-D density model, we developed a code based on a point mass algorithm as well as prisms algorithm that uses some easy to implement optimizations. The synthetic geoid undulations were compared to an existing gravimetric quasi-geoid solution. The modelling showed that the difference in crustal structure between the East European Platform and Palaeozoic Platform explains significant part of the observed geoid undulations in the study area. The remained residual could be related to density inhomogeneities in the upper mantle.

Gravity anomalies and Earth structure, Crustal structure, Europe

1 Introduction

In Europe, lithospheric intraplate deformations are defined by two major geodynamic processes: one is convergence of Africa and Eurasia and second one is glacial isostatic adjustment (GIA) of Fennoscandia. As shown by Marotta and Sabadini (2004), the first process plays the dominant role in the part of Europe extending from the Alps up to about latitude 52°N, while in latitudes higher than 58°N the effect of GIA prevails. In this separation of stress field, the southwestern margin of the East European Craton (EEC) (Fig. 1) is crucial, as it shields the northeastern regions from the compressional effects of Africa–Eurasia convergence. In the area adjoining this margin (it is called the Trans European Suture Zone, TESZ), the strain rate is clearly affected by both mechanisms (Marotta & Sabadini 2004; Marotta et al. 2004), therefore, both these processes should be taken into account in modelling of geodynamic processes there. Until recently, models of GIA were based mainly on 1-D rheological models of the Earth, with radial distribution of mantle viscosity and density. Recent studies showed, however, that the applicability of 1-D earth models to GIA modelling is limited (Whitehouse et al. 2006), and realistic 3-D distribution of density and viscosity inside the Earth is necessary.

Figure 1.

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

In the present paper, we investigate the 3-D density structure of the lithosphere in the area adjoining the TESZ. There are 18 published high quality 2-D P-wave velocity models based on results of new experiments and re-interpretation of old data of wide-angle reflection and refraction experiments LT (Grad et al. 2005), TTZ (Janik et al. 2005), POLONAISE'97 (Grad et al. 2002), CELEBRATION 2000 (Guterch et al. 2003) and SUDETES 2003 (Grad et al. 2003b). These models make a good basis for compiling a 3-D P-wave velocity model of the crust in the territory of Poland. However, one needs to know also density, to make proper tectonic and geological interpretations and modelling of geodynamic processes. The 3-D density distribution in the deep Earth can be obtained by two different ways. The first one is recalculation of known seismic P-wave velocity models into the density models using some relationship between density and seismic velocity (e.g. Birch 1961; Krasovsky 1981). The second way is to use inversion of the gravity field constrained by a priori information from seismic models (Kozlovskaya et al. 2004a). In both cases, usage of seismic models is a necessary condition.

Many publications devoted to regional-scale 3-D density modelling are based on interpretation of the Bouguer anomaly (Artemjev et al. 1994; Glaznev et al. 1996; Lefort & Agarwal 2002; Yegorova & Starostenko 2002; Kaban et al. 2003; Kozlovskaya et al. 2004a; Ebbing et al. 2006). The gravity field after Bouguer reduction represents the first vertical derivative of the gravitational potential caused by an irregular mass distribution within the Earth. This type of gravity field is very sensitive to the density distribution in the crust and is very convenient for interpretation of shallow structures, like sedimentary basins. Another type of observable gravity effects caused by a 3-D density distribution in the crust and upper mantle are variations of height of geoid. The geoid is defined as a surface with the constant value of gravity potential. While the Bouguer anomalies are inversely proportional to the square of the distance between a measurement point and an anomalous mass, the geoid, as an equipotential surface, is inversely proportional to the distance from the anomalous mass. Because of this, one can use geoid undulations to examine variations of density at larger depth, caused by uplift of the asthenosphere, deep roots beneath shield areas, mantle plumes, etc. Examples of such studies can be found in papers by Marquart and Lelgemann (1992) and Ebbing et al. (2006). Therefore, combined use of Bouguer anomalies and geoid undulations is a power tool to study the density structure of the lithosphere (cf.Ebbing et al. 2006). Such a modelling was not possible to do for our study area, however. The main problem was absence of digital Bouguer anomaly map compiled with uniform reduction density and resolution for the whole our study area. In addition, the seismic velocity models used in our study were not detailed enough for precise modelling of the Bouguer anomaly.

In our paper, we investigate anomalies of geoid height caused by density variations in the crust and upper mantle beneath Poland. Complicated structure of the crust in Poland has been revealed recently by new controlled-source seismic studies (e.g. Grad et al. 2003a, 2006; Majdański et al. 2006). However, knowledge about the structure of the mantle and the lithosphere-asthenosphere boundary is limited to several studies (Ladenberger et al. 2006; Pushkarev et al. 2007). From the other hand, results of TOR temporary passive seismic array in 1996–1997 demonstrated variations of seismic velocities and lithosphere thickness beneath Palaeozoic Platform (PP), TESZ and Fennoscandian Shield (Cotte & Pedersen 2002; Plomerová 2002; Shomali et al. 2006). Thus, one can expect that the asthenosphere uplift and density variations between different units of the lithosphere would be reflected also in the geoid shape. In addition to seismic studies, gravity modelling based on interpretation of the Bouger anomaly demonstrated that density inhomogeneities can exist in the mantle beneath Poland (due to asthenosphere uplift or to compositional variations). For example, 2-D gravity modelling by Krysiński et al. (2000) along POLONAISE and LT profiles suggested the presence of density anomalies in the mantle beneath these profiles, although this study did not explain origin of these anomalies. Yegorova and Starostenko (2002) and Yegorova et al. (2007) found positive density anomaly in the upper mantle beneath the Polish Trough. The analysis of the inhomogeneities in the crust and upper mantle in the territory of Poland was also presented by Grabowska et al. (1998).

Thus, the main motivation for our geoid study was to access the 3-D density structure of the crust and upper mantle in the area of TESZ. Recently, the similar problem was addressed in more large-scale studies (e.g. Artemieva et al. 2006; Tesauro et al. 2007). The other motivation was a great number of recently published information about the crustal structure obtained by the deep seismic sounding experiments (e.g. Grad et al. 2006; Majdański et al. 2006), together with a high interest in the geophysical community to regional gravity studies.

In our study, we demonstrate how geoid heights can be used instead of Bouguer anomalies in lithosphere studies. We describe the results of our study in two separate papers. The present one is devoted to modelling of the effect of inhomogeneous crust beneath Poland on the geoid undulations and the second one investigates the inhomogeneities in the upper mantle.

2 Calculation of Geoid Undulations from Large-Scale 3-D Density Distributions

2.1 Parametrization

One of the most common ways of parametrization of regional 3-D density models are rectangular regular grids (Glaznev et al. 1996; Sobolev et al. 1997; Kozlovskaya et al. 2004a). The gravity potential of models described with such grids can be calculated in the first approximation as a potential of a mass point. For precise calculation, we could use the simplified prism algorithm (Hjelt 1974). This algorithm is more time consuming, but it gives better results than point mass calculation for sources close to observation points, because it include also a shape of the blocks.

These two schemes can be easily implemented to calculate total gravity potential from given mass distribution. However, it cannot be compared to the observed values of gravity potential, because the former includes influence of masses outside the known mass distribution. To compare observed and synthetic values, it is necessary to take into account also the effect of mass distributions outside the study area. The simplest and the fastest way is calculating the average density in given layer of the model grid and subtracting it from the values in each block (Kozlovskaya et al. 2004a). The whole process is repeated in each layer of the model separately. This allows us to calculate directly the residuals of the potential caused by the residual mass distribution with respect to regional layered model.

The final step of calculation of the geoid is an implementation of the Bruns formula:

1

where V—calculated potential, γ—gravity acceleration according to reference ellipsoid, which gives the geoid values in metres. In this work, we use the GRS80 reference ellipsoid (Geodetic Reference System 1980). We also assumed that all calculation points, in which we compare the calculated and observed values, were located on the reference ellipsoid. To find the influence of the whole model in the single point on the ellipsoid, we need to sum influences of all blocks. Obtained values of the geoid can be finally compared with the observed gravimetric geoid data.

2.2 The Cartesian grid

For small areas (distances smaller than 100 km), we could use the Cartesian coordinates system assuming the flat Earth's surface. In this case, one can use cubic or rectangular prism model blocks and easy way to implement the prism algorithm. As mentioned above, the point mass scheme is much faster than the prism algorithm, and the differences in obtained values are lower than computation precision for distances between block centre and calculation point much larger than the block size. To optimize calculation time and final precision, we used switch of algorithms at given distance. For distances between block centre and calculation point smaller or equal to 3.5 of the vertical block size, we used the prism algorithm, while for larger distances the point mass algorithm was applied.

The second optimization is possible, if we use the fact that a single block influence in the point mass algorithm depends only on the distance and density of the block, but it does not depend on location of the block (Fig. 2a). In the case of a regular XY grid, there will be several blocks with the same distance to the calculation point. Thus, it is necessary to calculate only influence of one column of blocks (blocks on different depths) and then multiply it by proper density values from each column. In the worst case, it is necessary to calculate influences of a half of nodes in the area (triangle of blocks) and blocks along a diagonal (see Fig. 2a) in one calculation point in the corner, to estimate all possible distance combinations. This optimization is very efficient, if it is combined with fast sorting (‘quicksort’) and searching algorithms.

Figure 2.

Regular grid (a) has a strong symmetry. There are the same distances to eight blocks from one point. The gravity potential only from the distance and density, so it is necessary to calculate only one influence for all eight of blocks marked with arrows. To find all possible distance combinations, it is necessary to calculate distances from the first block in the corner (black dot) to blocks marked with empty circles. In the spherical coordinates (b) model is regular when using angles between blocks. It is possible to find symmetry also in this parametrization, but it is much weaker. (c) True distance (AC) between the centre of the model block and the calculation point at the ellipsoid was calculated using simple geometry. Earth's radii for points A and B include ellipsoidal curvature difference from GRS80; ang—angle calculated using Heversine formula.

Implementation of both schemes gives a very fast and precise calculation algorithm. Testing it with the model with 21 × 24 × 50 blocks showed that it is 30 times faster than the standard prism algorithm. Such a scheme can be used in interpretations, in which a lot of forward model calculation is required (inversions).

2.3 The spherical grid

The flat earth model cannot be used in large areas. In this kind of interpretation, it is necessary to use a real shape of the Earth. However, it is more convenient to use the spherical coordinates system or the geographical coordinates in our case. In this parametrization, the model is not regular anymore in the sense of distances. For simplicity, we chose regular grid in geographical coordinates. To implement point mass algorithm, we need to calculate distance between centre of a model's block and calculation point at the reference ellipsoid (GRS80) and volume of a block. The volume was calculated as

2

where ϕ is latitude of the block centre, R_i is distance to the block centre from the Earth's centre and dx and dy are block dimensions given in radians. The calculation of distances between block and calculation point was much difficult. The Cosinus law normally used is such a case cannot be used because of limited computer precision in handling very small angles, observed for tenths of kilometres on the Earth' surface. To solve this problem we used the Heversine formula,

3

which gives an angle from the Earth's centre between two points on the Earth's surface. Knowing the angle between these points and the ellipsoid radii in both points, it was possible to calculate the distance using a simple geometry (Fig. 2c). The ellipsoid radii were calculated as

4

where f is 0.003352810681225 according to GRS80. Because of irregular distances between blocks, the calculation optimization used in Cartesian coordinates system cannot be used. However, also in this case some symmetry can be noticed. In Fig. 2(b), we see the same distances between blocks located at the same latitude. It is possible to pre-calculate the influence of all blocks in the one point and use these values for all calculation points at the same latitude. However, the volume of each model block is different and has to be calculated every time. That is why this optimization is much less efficient.

3 Geoid for the Territory of Poland

The gravimetric geoid for the territory of Poland was elaborated on the basis of satellite measurements (Łyszkowicz & Denker 1994). The set of distances of the observed geoid from the reference ellipsoid GRS80 was defined in the nodes of the grid 1′ × 1′, and it included the area between 48°00.5 and 55°59.5′N and 13°00.5 and 24°59.5′E.

The layout of the observed geoid is presented in Fig. 3. The values of the observed differences between the geoid and the ellipsoid GRS80 in the study area vary between 23 and 48 m and in Poland itself between 28 and 43 m. The area of Palaeozoic Europe is characterized by larger anomalies (N > 40 m) than the area of East European Platform (N > 30 m), whereas the gradient zone between them corresponds to the TESZ area. As the size of the area is approximately 8° and 12° in N–S and W–E directions, respectively, the longest wavelength of a geoid signal there is l = 30. This suggests that only harmonics with numbers of l > 30 should show signals totally inside the study area, while harmonics with smaller numbers would produce a linear trend or a constant shift from the average.

Figure 3.

Observed geoid undulations from the reference ellipsoid GRS80 (according to Łyszkowicz & Denker 1994). The gradual change between the EEC and the Palaeozoic Platform is observed in the area of TESZ. Black lines marks position of profiles (see Fig. 4).

Therefore, the apparent linear trend in the observed geoid (Fig. 3) corresponds to harmonics with wavenumbers of l < 30. However, such long-wavelength gravity anomalies are not necessarily caused by deep structures beneath the lithosphere or outside our study area. They are often associated with local factors, like strong variations in crustal structure (Wang 1998) or areas of GIA (Simons & Hager 1997). Thus filtering of the long-wavelength component (or linear trend) from the geoid signal can remove the contributions from long-wavelength shallow sources. In our study area, such a local source can be a large contrast in thickness and structure of the crust in the TESZ revealed by seismic studies (Grad et al. 2002, 2003b, 2005; Guterch et al. 2003; Janik et al. 2005). In the next chapters, we estimate the effect of this feature on geoid shape.

4 3-D Velocity Model of the Crust for Poland

4.1 Seismic results used to construct a 3-D velocity model

There are 18 published high quality 2-D velocity models of the crust obtained by reflection and wide-angle refraction experiments LT (Grad et al. 2005), TTZ (Janik et al. 2005), POLONAISE'97 (Grad et al. 2003a), CELEBRATION 2000 (Guterch et al. 2003; Grad et al. 2006) and SUDETES 2003 (Majdański et al. 2006). Location of profiles used in our study is shown in Fig. 4. We collected all of them and read the P-wave velocity values every 10 km along profiles and every 1 km in depth starting from the sea level down to 50 km. As shown in Fig. 5 (top) for S02 profile model, the velocity grid, defined in such a way, is dense enough to sample important structures in the crust, thus using a denser grid was not necessary. Poorly resolved areas in some of the models (like velocities of 7.6 km s⁻¹ in the lower crust in S02 profile) were not taken to further interpretation. As a result, we gathered the data set of 60 794 values of V_p and interpolated this data set into a 3-D velocity model of the crust in Poland. The model was parametrized as a regular grid model in geographical coordinates with block sizes of 0.5° in WE direction, 0.333° in NS direction and 1 km in vertical direction. The WE block size in kilometres changes depending on its location from 37.1 km at latitude 48°N to 31.9 at 55°N. The denser grid was not necessary, because this size (about 30 km) corresponds to the velocity inhomogeneities that can be recognized in the crust using wide-angle reflection and refraction experiments data, where the distance between shots was about 30 km. The whole model has 24 × 21 × 50 blocks and samples the area between 13°E and 25°E and 48°N and 55°N. The ‘0’ depth value was attributed to the sea level.

Figure 4.

Location of deep seismic sounding wide-angle reflection and refraction profiles in the study area. The points every 10 km along profiles where 1-D velocity models (with values on every 1 km depth) were sampled, are indicated by black dots.

Figure 5.

The model for the profile S02 is presented (top) with 10 × 1 km grid. The grid nodes show the location of velocity values taken into 3-D model interpolation. Comparison of 2-D ray tracing model along P4 profile (middle) and a cross-section through 3-D model (bottom) obtained from all profiles by interpolation. The position of the cross-section corresponds to exact location of the P4 profile.

For interpolation, we used the BLOXER software (Pirttijärvi 2004), in which a two steps scheme is implemented. In the first step, it assigns values for those blocks that contain at least one point, or a mean value, if there is more than one point in a block. In the second step, blocks that contain no data points are filled with values using nearest neighbourhood algorithm. Using such an algorithm ensures continuity of the model velocities in the areas between profiles. To check the interpolation precision, we compared oblique slices along the profiles with the model for the profiles. An example for the P4 profile is shown in Fig. 5 (bottom).

It is hard to estimate the precision of this interpolation for every point of our 3-D model, however. The precision of 2-D ray tracing models based on wide-angle experiment data is estimated as ±0.1 km s⁻¹ for the velocities and ±1–2 km for boundaries (Janik et al. 2002; Grad et al. 2003a, 2006). Because our model is smooth (without boundaries), we estimate that the precision of our velocities cannot be better than ±0.1 km s⁻¹ along profiles. Obviously, this value increases with the distance from profiles, although it cannot be estimated precisely.

4.2 Major boundaries in the crust

The structure of the crust in the territory of Poland is characterized by strong variations. The crust of the EEC (based on CEL05 results; Grad et al. 2006) in the NE part is about 50 km thick with uncomplicated three-layer structure with flat boundaries between them. A very thin sedimentary layer overlies a crystalline crust with velocities 6.1 km s⁻¹ at the depth of less them 1 km. In the crust, there exist two boundaries with increase of velocities from 6.5 to 6.6 km s⁻¹ at 23 km and 6.6 to 6.9 km s⁻¹ at 34 km. The velocity change at the Moho boundary from 7.1 to 8.2 km s⁻¹ at the depth of 48 km is more significant. The EEC type of the crust is rapidly changing when we move to the TESZ (based on P4 interpretation; Grad et al. 2003a). There are thick (over 10 km) sediments with complicated velocity structure, underlain by a layer of metamorphosed sediments (down to 20 km) with velocities of 5.8 km s⁻¹. Two layers in the crystalline crust with velocities of 6.25 km s⁻¹ in the upper crust and 7.0 km s⁻¹ in the lower crust are divided by boundary at a depth of about 29 km. The Moho with extremely large velocity jump from 7.0 to 8.4 km s⁻¹ (along P4 profile) is located at a depth of 39 km. The crust of the TESZ differs between different terranes (Malopolska, Kuiavia and Pomerania) but ever significant change we observe between them and the Variscian crust in the SW Poland. This crust consists of two layers (based on LT profiles; Grad et al. 2005) with average velocities of 6.2 and 6.6 km s⁻¹, subdivided by a boundary at 24 km and covered with 1–5 km thick sediments. The Moho boundary at the depth of 30 km is characterized by a large velocity change from 6.7 to 8.2 km s⁻¹. Further to the SW, within the Bohemian Massif, there is no sedimentary layer and the upper crust velocities vary from 6.0 to 6.15 km s⁻¹ (based on S02 profile; Majdański et al. 2006). The middle crust with velocities of 6.55 km s⁻¹ extends from 20 to 31 km. The velocities in the lower crust are about 6.9 km s⁻¹. The Moho is located at the depth of 35 km and is characterized by change in velocity from 6.9 to 7.9 km s⁻¹. In the southern part of Poland, we observe an orogenic crust of Carpathians and Sudetes. These units have similar about 32 km thick three-layer crusts with velocities of upper crust of 4.4–5.6 km s⁻¹ for Carpathians and 5.8–6.2 km s⁻¹ for Sudetes, a middle crust with velocities from 6.2 to 6.4 km s⁻¹ and a lower crust with velocities from 6.5 to 6.9 km s⁻¹ (based on CEL05 and S03 results; Grad et al. 2006; Majdański et al. 2006). The boundaries between crustal layers are at 19 and 26 km and the Moho boundary is located at 32 km, with the upper-mantle velocities of 7.95 km s⁻¹ in both cases. It is also important for our study that the territory of Poland is rather flat with an average elevation of 170 m, except of the mountain areas in the south, where elevation is locally going up to over 2000 m.

5 The Effect of Major Boundaries in the Crust on Geoid Undulations

When estimating the influence of the crust, it is obvious that we cannot use just one density-velocity relation for the whole model. There are two reasons for this: first, the density-velocity relationship is different for the sediments and for the bedrock (Schön 1996); second, heat flow studies showed that thermal regime for Palaeozoic and Proterozoic geological units of Poland is significantly different (Majorowicz 1978; Majorowicz et al. 2003).

Therefore, we divided the crust into three layers on the base of the P-wave velocities: the most external one, including masses over the sea level (topography), the sedimentary layer and the one corresponding to the crystalline crust. The entire influence of the crust on the geoid is the result of effects of these three layers:

5

where N_topo is the effect of topography, N_sed is the effect of the sedimentary layer and N_cc is the crystalline crust effect. To calculate synthetic geoid undulations caused by a 3-D mass distribution within the crust, we developed an optimized computer algorithm presented in chapter 1.

5.1 The topography effect N_topo

Topography is one of the most crucial factors affecting the shape of the geoid. When defining its influence, the data about topography were taken from the NOAA base (NOAA 1988). The area of the research is characterized by a varied topography. About 3/4 of the area is flat, with the average altitude of 170 m. Mountains located in its southern part, with the highest altitude of over 2000 m. Due to heterogeneous topography, a greater resolution of the model was needed, that is why it was parametrized by the grid of 1 × 1 km along the axes x and y, while the elevation was transformed into 10 m blocks with the constant density of 2.5 g cm⁻³ equal to the average density of the sediments (Grabowska & Bojdys 2004). The topography map together with its influence of the geoid shape is shown in Figs 6(a) and (b). It can be noticed that the topography effect is larger for the mountain area, where it causes the uplift of the geoid up to 18 m. As for the remaining area, the effect is smaller and it is about 5–10 m.

Figure 6.

Topography map (a) (according to NOAA 1988) with calculated topography influence on a geoid (b). The topography model used for calculation had horizontal grid with 1 km spacing and vertical grid of 10 m. The constant density value of 2.5 g cm⁻³ for topography was used. The map illustrating the thickness of the sedimentary cover (c) and calculated effect on the geoid (d). For this calculation, the sediments densities were found using eq. (3), and the rest of the model was filled with constant density 2.7 g cm⁻³. The depth of the Moho boundary based on the 7.9 km s⁻¹ isoline in our 3-D velocity model (e) and its effect on the geoid (f) calculated with constant densities above (2.7 g cm⁻³) and below (3.35 g cm⁻³) of the Moho.

5.2 The sedimentary layer effect N_sed

Sediments are an important part of the Earth's crust, especially in the TESZ zone where their depth reaches even 20 km. Thus, their influence is crucial for proper modelling the gravity field. In our study, we defined the thickness of the sedimentary layer using the values of P-wave velocity in the uppermost crystalline crust from the 3-D velocity model (Grad et al. 2006; Majdański et al. 2006). Namely, the sediments were defined as rocks with P-wave velocities lower than 5.8–6.0 km s⁻¹, depending on the area. The map illustrating the thickness of the sedimentary layer (Fig. 6c) shows a deep sedimentary basin running through the middle of the model that corresponds to the TESZ.

The sediments density depends on many factors. The mineralogical composition, porosity and quantities of the material that fills the pores (saturation condition) have the greatest effect. There is a strong correlation between the rocks density and porosity. The increase in porosity leads to decrease in density and is controlled by the density of fluids filling the pores. The relative increment in density diminishes with depth, due to compaction of the material and closing the pores under the influence of lithostatic pressure. At a depth of approximately 4 km, the effect of porosity is small and the density of sedimentary rocks approaches a certain asymptotic value (Dortman 1976; Kopf 1980; Jelic 1984; Schön 1996).

To specify the influence of sediments on the gravimetric geoid, we used a general dependence of sediments density on depth that assumes an exponential increase of density with depth:

6

We estimated coefficients A, B and C in eq. (6) using the density in sedimentary layer along the profile P4 obtained by Grabowska and Bojdys (2004), who used the data from boreholes in their study. The density values were taken at every 50 km along the profile and at every 0.5 km depth. As the P4 profile crossed the most important crustal units of Poland, with different type of deep structure, this density-depth relationship can be applied to the whole model of the sedimentary layer.

Fig. 7 shows the distribution of the values of density with depth taken from the model for the profile P4, along with the line that corresponds to the elaborated density-depth relation. The estimated dependence between the density of the sedimentary layer and depth is

7

Figure 7.

The density-depth relation for sediments calculated for P4 model (according to Grabowska & Bojdys 2004). Black points mark the values read from P4 model every 50 km along and every 0.5 km depth. The line shows the relation (7).

The calculated coefficient of correlation between the values taken from the model by Grabowska & Bojdys (2004) and those calculated using the formula (7) was 0.81.

To calculate the gravity effect of the sediments (Fig. 6d), we transformed the velocities to the densities using eq. (7) and filled the rest of the model with constant density of 2.7 g cm⁻³. As can be seen (Fig. 6d), areas with maximum thickness of the sediments correlate well with the areas where the surface of the geoid is undulated with respect to the reference ellipsoid. One can immediately notice the TESZ, where this difference reaches up to −2 m. Another distinctive area is the Bohemian Massif. Lack of the sedimentary layer in this region results in the uplift of the surface of the geoid with respect to the reference ellipsoid, by about 5 m.

5.3 The effect of the crystalline crust N_cc

Similar velocity criterion was used to define thickness of the crystalline crust. We attributed to the crystalline crust the blocks with velocities between 6.0 and 7.9 km s⁻¹. As the crystalline crust in Poland consists of several geological units with different composition, age and tectonic history, it cannot be interpreted as a single layer. This problem required a more detailed consideration.

5.3.1 The influence of the shape of the Moho boundary

The gravitational effect of the crust depends mainly on topography of boundaries having a sharp density contrast. That is why the influence of the Moho boundary is particularly important. The depth of the Moho was calculated on the basis of the velocity isoline of 7.9 km s⁻¹ (Fig. 6e). It is seen that the research area can be generally subdivided into two major parts: the crust of the EEC, with its typical thickness of 45–50 km and a thin crust of the Palaeozoic formations in the Western Europe, where crust is typically 30–35 km thick.

The effect the Moho boundary on the shape of the geoid is presented in Fig. 6(f). To exclude effects from other sources, we assumed that the density in the part of the model above the Moho (including sediments) has the constant density of 2.7 g cm⁻³ and the part below the Moho has the constant density of 3.35 g cm⁻³. This strong contrast was applied to show the maximum possible influence of the shape of the Moho boundary and it is not present in the final model. As can be seen, the shape of the Moho boundary is an important factor determining the synthetic geoid and it is the main reason for the linear trend in the shape of the geoid anomaly. However, this factor is not sufficient to explain all the details of the observed geoid shape. It can be concluded that heterogeneous composition of the crust is quite crucial, although its gravitational effect is weaker compared to the effect of shape of the density discontinuity.

The influence of the boundary shape on the synthetic geoid depends on the accuracy of seismic models. In the case of ray tracing modelling, uncertainty of ±1.0 km is assumed for the depth of major reflecting boundaries. To estimate the influence of such uncertainty of the Moho on the geoid, we calculated the effects of this boundary for two artificial models with the constant densities of 2.7 and 3.35 g cm⁻³ above and below the Moho, respectively. In the first model, the whole Moho was uplifted by 2 km, and in the second one it was lowered by 2 km. The maximum difference between geoid calculated for models with the biased Moho and the ‘true‘ Moho was about 0.8 m.

5.3.2 The influence of the heterogeneous mass distribution in the crust

The aim of our crustal modelling was to choose the most optimal initial model of the crust that could be used later to analyse a density variation in the mantle beneath the Precambrian Platform and much younger structures of the PP and the Carpathians. To create proper density model of the crystalline crust, we had to use different density-velocity relations in different parts of the area, because individual tectonic units and crustal blocks may have different density-velocity relations (Kozlovskaya et al. 2004b). Based on the thickness of the sedimentary layer and the depth to the Moho boundary, we distinguished four major areas with significantly different velocity structures, namely, the EEC, PP, TESZ and Carpathians. To specify the boundaries of these units, we made the following analysis. First, we divided the 3-D velocity model of the crust into areas corresponding to major tectonic units from Fig. 1. Then in all points attributed to the same unit, we calculated 1-D P-wave velocity models and estimated an average 1-D velocity model for each unit. For the model nodes near the edges of the main units, we compared the correlation between 1-D velocity model in each node and averaged 1-D models for four main areas. That helped us to define the boundaries of the area more precisely. In addition, on the basis of the differences in the velocity structure beneath the TESZ, three minor blocks were distinguished: Pomeranian, Kuiavian and Malopolska. The Bohemian Massif was separated as a distinct unit in the PP, as a region characterized by very thin sediments. Schematic division of the area into major tectonic units is presented in the Fig. 8.

Figure 8.

Schematic division of the area into tectonic units based on 1-D velocity models. This division was used to recalculate the velocities to densities separately for each area. MP—Malopolska, KU—Kuiavia, PO—Pomerania.

To recalculate densities into seismic velocities, we used the following relations. The first one was the non-linear relationship by Krasovsky (1981):

8

This relationship proved to be efficient for the crust of cold regions, like shields and platforms, where the influence of temperature and pressure on density is minor compared to the effect of rock composition (Krasovsky 1981).

However, this relationship is not suitable for hot areas in the lower crust, where rocks are transformed into mafic granulites and mafic garnet granulites under high pressure and temperature. For such conditions the relationship by Sobolev and Babeyko (1994) can be used:

9

Selection of proper density-velocity relationship for the lower crust is very important, because it defines the density contrast at the Moho boundary. That is why we started our crustal modelling from testing the influence of density-velocity relationships on calculated geoid. First, we calculated the integrated effect of all the layers of the crystalline crust above the lower crust with velocities over 7 km s⁻¹. The densities in these layers were calculated according to the Krasovsky relation. In the remaining part of the model, the constant density of 3.35 g cm⁻³ was assumed. The estimated effect is shown in Fig. 9(a). After that we estimated the effect of different relations describing density-velocity dependence in the lower crust. Namely, the density in the upper and middle crust was calculated as in the previous case, but the density in the lower crust was calculated assuming density-velocity relationships described by eqs (8) and (9). The results are shown in Figs 9(b) and (c), respectively. Additionally, the effect of assuming the constant density of 2.97 g cm⁻³ in the lower crust was calculated (Fig. 9d). In all cases, a homogeneous density of the 3.35 g cm⁻³ was assumed in the upper mantle. The obtained result indicates that influence of inhomogeneities in the upper and middle crust is less significant than the mass distribution in the lower crust or the density contrast at the Moho boundary. Generally, local heterogeneities in the upper crust compensate each other inside it.

Figure 9.

The synthetic geoid calculated for a constant density of 3.3 g cm⁻³ for the velocities >7.0 km s⁻¹ (a), constant density of 3.35 g cm⁻³ in the mantle with different density-velocity relations for the lower crust: Krasovsky (b) and Sobolev and Babeyko (c), constant density of 2.97 g cm⁻³ in the lower crust (d).

To find the proper density-velocity relation for each tectonic unit, we used the assumption that Moho in the density model (defined by the density of 3.35 g cm⁻³) should be at the same depth as the Moho boundary in the seismic model (V_p > 7.9 km s⁻¹). For this purpose, we calculated average 1-D density models for four major tectonic units using different density-velocity relationships (8) and (9) and also constant density and compared them with correspondent 1-D P-wave velocity models. For example, the average thickness of the crust for PP was estimated to be 30–35 km (Grad 2003a). Assumption of the constant density in the lower crust leads to the appearance of the densities typical for the mantle at 40 km depth. When assuming Krasovsky's relation, these densities appear at 35 km. That is why Sobolev and Babeyko's relation seems to be the most appropriate for the lower crust of this region.

After such an analysis, we came to the conclusion that the relationship (8) suits better for the upper and middle crust and also for the lower crust of the EEC, while the relationship (9) suits better for the lower crust of the remaining part of the area. Calculating the difference between the effects obtained with different relations allows estimation of the uncertainty of synthetic geoid due to different density-velocity relation. The maximum value of the difference computed in this way does not exceed 3 m, and it is the estimate for uncertainty resulting from assuming an improper density-velocity relation.

5.3.3 The complete crust influence on the geoid

Using the empirical relations (8) and (9), we transformed the 3-D P-wave velocity model into the density model and calculated the synthetic geoid by summarizing the effects from the individual layers, namely, topography, the sedimentary layer and the crystalline crust. Fig. 10 presents comparison of the observed geoid (a) with the calculated effect from the 3-D crustal model (b), as well as the difference between them (c).

Figure 10.

The comparison of observed geoid (a), calculated effect of the whole crust (including the effect of topography, sediments and crystalline crust) (b) and the difference between them (c). Observed differences are still larger than a sum of the uncertainties of the effects of all crust components, therefore, we suggest the presence of the density variations in the upper mantle.

As can be seen, the shape of the calculated geoid reproduces rather well the main features of the observed geoid, including the long-period ones. The computed geoid undulations vary between 49 m for the Variscian Europe and 22 m in the Precambrian Platform region. The change of the crust structure from a relatively thinner for the PP and the Carpathians to much thicker in the EEC region results in a gradual change of the geoid height that occurs in the TESZ. The sedimentary layer and local heterogeneities in the upper crust are compensated isostatically inside it, as no correspondent local small-scale anomalies are seen in the geoid.

6 Discussion And Conclusions

In our study, we demonstrated that variable density structure of the crust has significant impact on the geoid shape in our study area. Density contrasts exist not only at the Moho boundary, but on all depth levels of the crust. These internal density inhomogeneities compensate each other to some extent, but not completely. Generally, contrasting crustal structure to the both sides of TESZ is the main factor controlling large-scale variations of the geoid undulations.

Our study demonstrates that geoid heights can be used not only for evaluation of deep mantle density inhomogeneities in global studies, but also in regional-scale crustal studies. In combination with the modelling of the Bouguer anomaly, such an analysis can possibly help to overcome potential inconsistency in gravity anomaly data sets. This approach will be used in the future, especially with the data from GRACE and GOCE satellite projects (Tapley et al. 2004), that would help to improve the Bouguer anomaly database for discussed area.

Our density model depends greatly on the initial velocity model. Its uncertainty generally results from the seismic interpretation along profiles, 2-D to 3-D interpolation, assumption of a density-velocity relation or subdivision of the model into tectonic units. The uncertainty due to topography effect estimated in our study is less than 1 m, while the uncertainty due to non-precise major boundaries in the crust is of about 1 m for the sediments boundary and about 2 m for the Moho boundary. However, all these uncertainties are still lower than the residual in Fig. 10(c), where one can see two significantly anomalous regions. The area with the largest positive residual (up to +6 m) is observed in the PP, while the largest negative residual (down to −8 m) is seen within the EEC. Having assumed that the 3-D velocity model of the crust is close to the true one, we may conclude that we need the assumption of a heterogeneous mantle to explain the observed gravity effect. Thus, the positive anomaly in the PP may be due to the presence of the layer characterized by a lower density (probably, asthenosphere) in the upper mantle. The negative anomaly in the area of the EEC could indicate a cold, thick and dense lithosphere. However, the main crustal suture in our study area (TESZ) is not seen in Fig. 10(c), which means that this structure has no continuation in the upper mantle. Thus, the separate density modelling is necessary to perform to explain the anomalies in the upper mantle. This will be described in the second part of our paper.

Acknowledgments

We would like to thank Prof. Adam Łyszkowicz for providing us with the data set of the geoid undulations in Poland. We are very grateful to Dr M. Kaban and Dr J. Ebbing, whose constructive comment helped to improve the earlier version of our paper. This work was partially supported by MNISW grant no. 2P04D05428. Majdański and Świeczak stay at Sodankylä Geophysical Observatory of Oulu University and were financed by the Academy of Finland (grants no. 112603 and 122645). For all geographical coordinates calculation and maps preparation, we were using GMT package (Wessel & Smith 1998).