Three-dimensional seismic modelling of the crustal structure between East European Craton and the Carpathians in SE Poland based on CELEBRATION 2000 data

Summary

In this paper, we present results of 3-D first-arrival tomography applied to data recorded in SE Poland during the CELEBRATION 2000 seismic experiment. The target area covers ca. 500 × 500 km² and represents a complex geological setting from old Precambrian platform (East European Craton, EEC), through the crustal blocks (terranes) that form the Trans-European Suture Zone (TESZ) to the young Alpine orogen—the Carpathians. Contrasting velocity distributions in various geological units makes the tomographic inversion challenging. For this reason, much attention was paid to the inversion methodology. We tested ‘multioffset’ and ‘multiscale’ approaches. By increasing the offset range in the ‘multioffset’ inversion we have independently constrained different depth ranges of our model. We have found that the ‘multiscale’ method, that is, the gradual stepping from bigger to smaller model cells produced the preferred solution. Resolution of the resultant crustal model was determined by performing checkerboard and restoration resolution tests. Lateral resolution of the order of 25 km is inferred down to 10 km depth and 50 km down to 25–30 km depth. We have also applied a quasi-Monte Carlo analysis to determine the absolute errors of model parameters, which to our knowledge, is the first such attempt in case of a 3-D wide-angle data set. The mean uncertainties of our tomographic velocities are 0.1 km s⁻¹, but some anomalies are recovered with larger errors. We also advocate that the results from the central part of the Małopolska Block (MB) are biased by the influence of crustal anisotropy as confirmed by an independent study. The overall picture of the transition from the EEC to the TESZ resembles an Atlantic-type passive margin, with several high velocity bodies occurring along the EEC edge and relatively low velocities (V_P < 6 km s⁻¹) extending down to ca. 20 km. Of the crustal units investigated, the Upper Silesian Block (USB) fully deserves to be called a terrane. It is distinct in terms of seismic velocities and its eastern edge (Kraków–Lubliniec Fault) is clearly visible and continues beneath the Carpathian nappes. The differences between the overall crustal structure of the MB and the Łysogóry Unit (LU) are minor and thus the ‘terrane’ concept does not necessarily apply to them.

Inverse theory, Controlled source seismology, Seismic tomography, Crustal structure

1 Introduction

Modern seismic acquisition capabilities allow the deployment of dense networks of long-range refraction/wide-angle reflection profiles in order to obtain detailed 3-D images of the crust and uppermost mantle. The world's largest deep seismic sounding project called CELEBRATION 2000 (Central European Lithospheric Experiment Based on Refraction) (Guterch et al. 2001, 2003) aimed at recognition of the deep geological structure from the Alps through the Pannonian Basin, Carpathian Mountains, Trans-European Suture Zone (TESZ) and the East European Craton (EEC). Interpretation of the CELEBRATION 2000 data started from 2-D modelling along the main profiles (see the results in Hrubcová 2005; Janik et al. 2005; Malinowski et al. 2005; Grad et al. 2006; Środa et al. 2006), but the final goal of the experiment was to derive a 3-D crustal model. Until now, the 3-D modelling based on the CELEBRATION 2000 data has been performed for the Eastern Alps (Behm et al. 2007) and for the Pannonian Basin area (Bodoky et al. 2003; Kovacs et al. 2004).

In this paper, we present results of first-arrival tomography applied to the data from the area of SE Poland (Fig. 1) where the densest CELEBRATION 2000 recordings were made (Fig. 2a). The resultant model together with the models based on the POLONAISE' 97 and SUDETES 2003 data sets (Środa et al. 2002; Majdański et al. 2007) provide the most comprehensive 3-D crustal image in Central Europe.

Figure 1.

Main geological units in the study area. HCF—Holy Cross Mountains Fault, GRF—Grójec Fault, KLF—Kraków - Lubliniec Fault Zone. Dots mark receiver positions of the CELEBRATION 2000 profiles. Red circles denote shotpoints. The inset shows study area on the background of the main tectonic units of Central Europe. Modified from Malinowski et al. (2005).

Figure 2.

(a) Area of SE Poland subjected to the 3-D tomographic modelling. Points mark receiver positions of the CELEBRATION 2000 project. Main profiles are marked by green. Red stars—shotpoints; (b) Ray paths for all shot-receiver pairs used in the tomography (14 719).

The primary objective of our work were to provide a well-constrained 3-D velocity model of the crust and uppermost mantle. This model should constitute a basis for more complex and integrated interpretations in the future, including joint refraction/reflection tomography (Malinowski 2006, Malinowski et al. 2008) and potential field modelling. In this paper, we focus primarily on the methodological aspects of the modelling as well as on the detailed discussion of model resolution and uncertainties. For estimating uncertainties a quasi-Monte Carlo method was successfully applied. Finally, we also provide the preliminary geological interpretation of our results.

2 Tectonic Background and Previous Studies

Deep seismic sounding investigations carried out in Poland since the late 1960s revealed an anomalous structure associated with the transition zone from the EEC to the Palaeozoic Europe (Guterch et al. 1986) later referred to as the TESZ (Berthelsen 1993).

The study area is located within the southern end of the Polish segment of the TESZ. The TESZ sensu lato constitutes a collage of crustal blocks which are usually interpreted to form terranes of different rank and provenance (either Gondwanan or Baltican). The terranes docked with Baltica at different times during the Palaeozoic (e.g. Pożaryski et al. 1992; Lewandowski 1993; Franke 1995; Winchester et al. 2002). In the study area three such units are usually distinguished: Łysogóry Unit (LU), Małopolska Block (MB) and Upper Silesian Block (USB) (Fig. 1).

USB continues further SW towards the Brno Block and the two are known collectively as Brunovistulicum (Dudek 1980; Buła et al. 1997)—one of the first tectonic units of Central Europe to be defined as a ‘terrane’ (Brochwicz-Lewiński et al. 1986) and still the best documented terrane in the TESZ. The USB itself constitutes a crustal-scale unit with a well defined eastern boundary—namely the Kraków–Lubliniec Fault Zone (KLF) as imaged recently by the CELEBRATION 2000 data (Malinowski et al. 2005). The southern part of the USB is covered by the Carpathian nappes.

The nature of the two remaining units is much more speculative. Both of them are accessible in outcrop only in the area of the Holy Cross Mountains and their southern part is buried below thick Miocene sediments of the Carpathian Foredeep. The Holy Cross Mountains Fault (HCF), which is believed to separate MB and LU, is not discernible in seismic data (Środa et al. 2006; Malinowski et al. 2005). A recent study of Środa (2006) shows that both units are affected by 8–10 per cent P-wave anisotropy confined to the upper crust.

In palaeogeographic reconstructions MB is taken as a fragment of Baltica displaced by ca. 1000 km in post-Emsian times (Lewandowski 1993) or as a part of Gondwana that was approaching Baltica in the mid-Cambrian (Belka et al. 2002). LU is envisaged either as a suspect Gondwana-derived terrane located close to the MB and the EEC in the Late Cambrian (Belka et al. 2002) or as a proximal Baltica-derived terrane which has not changed its position relative to Baltica since Late Precambrian times (Dadlez et al. 1994; Dadlez 2001). Recent analysis of sedimentological data (Jaworowski & Sikorska 2006) suggests however that both the LU and the MB were not detached from Baltica and thus they do not form allochtonous terranes. According to this interpretation they form together part of the Baltica (EEC) passive margin.

Results of 2-D modelling along CELEBRATION 2000 profiles that run through the central part of the MB and LU (Malinowski et al. 2005; Środa et al. 2006) have shown that those units are almost undistinguishable in terms of their seismic properties. The only exception is the area of Radom-Kraśnik High, an easternmost part of the LU (Dadlez 2001) which seems to have a well established crystalline basement (Janik et al. 2005; Malinowski et al. 2005). The remaining parts of the MB and LU are characterized by an unusual low velocity (V_P < 6 km s⁻¹) upper crust that is up to 18 km thick and is similar to the low-velocity TESZ infill imaged in northern Poland along other deep seismic profiles (Guterch et al. 1994; Janik et al. 2002; Grad et al. 2003). This low-velocity crust was attributed to the East Avalonia accretionary wedge in the northern part (Grad et al. 2002) or to the Neoproterozoic passive-margin metasediments in the southern part (Malinowski et al. 2005). What is also common for both the northern and southern segments of the TESZ is that its middle/lower crust resemble severely attenuated EEC crust (Grad et al. 2002). However, fundamental differences exist concerning the nature of the EEC edge itself. In northern Poland it is imaged as subvertical contact (e.g. Grad et al. 2003). In contrast, in our study area it forms a gentle step-like contact with characteristic occurrence of high velocity bodies (HVB) as interpreted from the 2-D profiles (Malinowski et al. 2005; Grad et al. 2006; Środa et al. 2006). Those HVB are also well reflected in the gravity and magnetic fields (Grabowska et al. 1998; Grabowska & Bojdys 2001). Hence, the overall geophysical signature of this area seems to be the analogue of an Atlantic-type passive continental margin observed at present (e.g. Holbrook et al. 1994).

3 Modelling

3.1 Seismic data

During the CELEBRATION 2000 experiment data were recorded not only along in-line profiles but also along cross-lines which provided good 3-D coverage. However, the real coverage is less than would be expected for all possible shot-receiver pairs due to the three separate deployments. Despite this, the ray path density in the most geologically complicated area is excellent (Fig. 2b). Data were recorded mainly by the single-channel RefTek-125 ‘Texan’ units, equipped with 4.5-Hz geophones, deployed on average spacing of 3 km. More than 60 shotpoints (of which 56 were used in this study) with charges of 100–800 kg explosives provided seismic energy.

The quality of cross-line record sections is comparable or sometimes even better than inline ones (compare Fig. 3 and Figs 4–5). In addition to the crustal refracted phases, both from the sedimentary cover (P_sedim) and consolidated crust (P_g), several mantle phases are observed including well-developed reflected waves. In this paper, we focus however on first-arrivals only. These were traced up to nearly 350 km distance from a source.

Figure 3.

Examples of inline record sections from different CELEBRATION 2000 profiles. Reduction velocity 8 km s⁻¹. Bandpass filter 2–15 Hz was used.

Figure 4.

Examples of cross-line record sections. Positions of the shotpoints and recording profiles are marked on the map. Reduction velocity 8 km s⁻¹. Bandpass filter 2–15 Hz was used.

Figure 5.

Examples of cross-line record sections. Positions of the shotpoints and recording profiles are marked on the map. Reduction velocity 8 km s⁻¹. Bandpass filter 2–15 Hz was used.

The input to the inversion was 14 719 hand-picked first-arrival traveltimes, representing P_sedim, P_g and P_n refractions. The picking uncertainty was determined by looking at the energy ratio in a 0.5 s long window centred at the picked arrival (Zelt & Forsyth 1994). The mean picking error amounts to 0.07 s with a standard deviation of 0.03 s, thus we assumed a constant picking uncertainty of 0.1 s. The scatter in the plot of all traveltimes (Fig. 6) illustrates that the large differences in crustal structure might be expected in the investigated area. Dispersion of traveltime data is highest for offsets between 100 and 150 km and amounts to ca. ±2 s, which might be partly explained by the difficulty and non-uniqueness in determining the first-arrival picks in that offset range. In order to prepare 1-D initial model we fitted a fifth order polynomial to our traveltimes and subsequently we calculated traveltime curves using program LAUFZEIT (Kaminski & Müller 1979). By means of a trial-and-error procedure we derived the starting model presented in the inset of Fig. 6. High gradient zone represents the sedimentary layer with velocities between 2 and 5.5 km s⁻¹. Moho was set as a transition zone between crust and mantle velocities in the depth range of 30–40 km.

Figure 6.

Plot of 14 719 first-arrival traveltimes used in the inversion (reduced at 8 km s⁻¹). Green line—traveltime curve. Blue line—traveltime curve calculated for the model shown in the inset.

3.2 Parameters and inversion strategy

For inversion of the first-arrival traveltimes we used the well-known FAST package (Zelt & Barton 1998). In this commonly used software, the forward solution is obtained efficiently using a finite-difference scheme for solving the eikonal equation (Vidale 1990; Hole & Zelt 1995). The inverse problem is regularized by use of the vertical and horizontal ‘roughening’ matrices (second spatial derivatives). The solution of the inverse part is sought by the iterative LSQR solver (Paige & Saunders 1982). The trade-off parameter λ between data residuals and model norm (roughness) is automatically determined at each iteration. For all calculations we used a cubic grid of 1 km spacing between nodes, as this allows efficient determination of traveltimes at the precision needed for the inversion. There are several ‘tuning’ parameters in the FAST code of which the most important is the s_z parameter. It determines the ratio of vertical versus horizontal model smoothness and varies from 0 to 1 (0 means no vertical smoothing and 1 means no horizontal smoothing). After several tests we chose s_z = 0.2 which is consistent with the values considered by Zelt et al. (2003) as representative for wide-angle data. Parameters of the modelling are summarized in Table 1. An additional 8 km were included at the top of the model to properly handle the elevation and the smoothing of model edges, resulting in a total model thickness of 88 km.

Table 1.

Parameters of the tomographic modelling.

As the study area is large and characterized by complex geology, we paid special attention to selection of a proper inversion strategy. A 1-D starting model might not be representative for the whole area. Thus, to reduce the risk that during the inversion, models would go outside the ‘linearity’ of the starting model, we used two special approaches referred to as ‘multioffset’ and ‘multiscale’.

The ‘multioffset’ strategy was adapted after Środa et al. (2002). Increasing the maximum offset allowed in the inversion, gradually increases the depth within the model affected by the inversion as longer offset refractions penetrate to greater depths. We used evenly spaced offset bins of: 50, 100, 150 and 200 km and for the final run we included all the offsets up to 350 km. For each bin four non-linear iterations were run. We also used two inverse model parametrizations: 10 × 10 × 4 km and then 5 × 5 × 2 km. The overall ‘multioffset’ modelling flow is illustrated in Fig. 7.

Figure 7.

Flow chart of the ‘multioffset’ inversion.

Seismic tomography is a typical mixed inverse problem; that is, simultaneously over and underdetermined. Additionally some of the model cells are crossed by several rays whereas others are not penetrated by rays. The ideal model parametrization should be flexible enough to account for this fact and should be adapted to the expected ray density. However, it is difficult to obtain flexible parametrizations using regular model grids. Alternatively, Voronoi polygons or Delauney triangles might be used (Böhm et al. 2000; Trinks et al. 2005). Another possibility is to use multiscale tomography (MST), in which the solution is sought on different scales of regular grids (Zhou 2003, 2004).

In this study we are not using MST sensu stricto, but a modified form. We are performing ‘single scale’ tomography (i.e. tomography with constant model parametrization), gradually decreasing model cell size. For each model parametrization (50 × 50 × 8 km, 25 × 25 × 8 km and 12 × 12 × 8 km) we ran one non-linear iteration testing 6 λ values. Only for the final parametrization (6 × 6 × 2 km) we ran four non-linear iterations. We call this modelling approach ‘multiscale’ one (see the flow in Fig. 8).

Figure 8.

Flow chart of the ‘multiscale’ inversion.

3.3 Results of the inversion

If one looks at the χ² values only (Tables 2–4), the preferred solution would be the ‘multioffset’ solution (χ² value of 4.7 and 2 for coarser and finer parametrizations, respectively), not the ‘multiscale’ solution (χ² = 5.9). These differences are insignificant if we consider the number of iterations needed to obtain a similar level of data matching (20 or 40 in the case of ‘multioffset’ or only seven in the case of ‘multiscale’ inversion). Slices through the ‘multioffset’ solution (Fig. 9) reveal some ‘spike’ anomalies that are hard to be geologically justified. The ‘multiscale’ model is smoother (see the slices in Fig. 10), but there is generally good agreement between both inversion results. However, if we look at the simple 1-D averages of both solutions and compare them with starting model, the ‘multiscale’ seems to be closer to the initial model (Fig. 11). Thus, on the expense of slightly higher residuals, we choose the ‘multiscale’ solution as the preferred one.

Table 2.

χ² and rms values for ‘multioffset’ inversion (parametrization 10 × 10 × 4 km).

Table 3.

χ² and rms values for ‘multioffset’ inversion (parametrization 5 × 5 × 2 km).

Table 4.

χ² and rms values for ‘multiscale’ inversion.

Figure 9.

Horizontal slices through the ‘multioffset’ solution in the depth range 0–55 km. Velocity contours 0.25 km s⁻¹. Recording profiles are indicated by pale yellow lines. Depths in all subsequent plots are referenced to the sea level.

Figure 10.

Horizontal slices through the ‘multiscale’ solution in the depth range 0–55 km. Velocity contours 0.25 km s⁻¹. Arrows indicate anomalies that were probably recovered with large errors.

Figure 11.

Comparison of the starting model and the 1-D averages from (a) ‘multioffset’ and (b) ‘multiscale’ solutions.

The ‘spike’ character of the ‘multioffset’ solution results also in more patchy ray distribution (Fig. 12). However we cannot directly compare the ray density for both solutions due to their different parametrization (5 × 5 × 2 km cells for ‘multioffset’ and 6 × 6 × 2 km cells for ‘multsicale’ inversion), the ‘multiscale’ solution seems to provide more even ray coverage down to 25–30 km depth (Fig. 13). Both approaches are characterized by highest ray density in the 5–10 km depth range. Models are practically unconstrained below 55 km depth, where only a few rays are found.

Figure 12.

Horizontal slices through the ray hit-count distribution for the ‘multioffset’ solution in the depth range 0–55 km.

Figure 13.

Horizontal slices through the ray hit-count distribution for the ‘multiscale’ solution in the depth range 0–55 km.

In the vertical slices (Fig. 14) we marked the approximate Moho location following the 7.6 km s⁻¹ isoline. It is connected with the fact, that in a continuous tomographic velocity field an interface should be placed between mean velocities above and below an interface. For the crust–mantle boundary, assuming mean velocities of 7.2 km s⁻¹ in the lower crust and 8 km s⁻¹ in the sub-Moho mantle, we obtain 7.6 km s⁻¹ as a ‘tomographic Moho’. This ‘Moho’ is largely consistent with the results of the detailed 2-D forward ray trace modelling along the main CELEBRATION 2000 profiles and 3-D reflection/refraction tomography (Malinowski 2006; Malinowski et al. 2008).

Figure 14.

W–E vertical slices through the ‘multiscale’ solution. Dotted white line marks approximate Moho position (see the text for explanation). Velocity contours 0.25 km s⁻¹

3.4 Indications of crustal anisotropy?

The high level of rms residuals after several non-linear iterations (see Tables 2–4) indicates a very low convergence of the inverse problem. For instance, when we decreased the final cell size in the ‘multiscale’ solution to 3 × 3 × 2 km and we ran an additional four iterations, the χ² dropped from 5.9 to a value of 4. Such a low convergence might be a result of the data itself or other ‘structural’ factors. As noted by Hobro (1999), it may also indicate the influence of anisotropy. We suggest that this is the most plausible reason for our relatively large residuals. Środa (2006) analysed the azimuthal distribution of P_g arrivals from the MB and showed that the area of the MB is characterized by 8–10 per cent P-wave velocity anisotropy in the crust down to 15 km depth with the fast propagation axis at 115°, coincident with the structural axis of folds that outcrop in the Holy Cross Mountain area. The azimuthal plot of the final rms residuals for the ‘multiscale’ solution clearly indicates a similar fast direction (Fig. 15). The effect of the anisotropy can actually be traced in the 3-D model by comparing the inline results from various azimuths. In 2-D models from profiles roughly perpendicular to the fast axis (CEL01 and CEL02), we do not observe any increase in velocities in the central part of the MB (Malinowski et al. 2005; Środa et al. 2006). However, a high-velocity upper crust (ca. 6.2 km s⁻¹) is observed along the CEL14 profile, which runs along the fast axis. Thus, we believe that the velocity anomaly of 6.2 km s⁻¹ occurring near the crossing point of all these profiles at 10 km depth (Fig. 10), likely originates from this anisotropy.

Figure 15.

Azimuthal distribution of the rms traveltime residuals for the ‘multiscale’ model.

3.5 Resolution and uncertainty analysis

To assess the resolution of the velocity model, we used both the checkerboard (CRT) and restoring resolution (RRT) tests. In the CRT we added sinusoidal velocity anomalies to the starting 1-D model with the maximum amplitude of ±0.25 km s⁻¹ and dimensions of 25 × 25 or 50 × 50 km. The synthetic data were generated for both test patterns using the actual shot and receiver geometry. Gaussian noise with σ = 0.1 s (equal to the picking uncertainty) was added to the synthetic data set which was then subsequently inverted. Results of the CRT indicate that the resolution of at least 25 km is maintained down to ca. 10 km depth (Malinowski 2006), whereas lateral resolution of 50 km is inferred down to 25 km in the central part of our model (Fig. 16a).

Figure 16.

Horizontal slices through the results of resolution tests: (a) CRT; (b) RRT. Velocity contours 0.1 km s⁻¹.

The potential weakness of the CRT is connected with the fact that the ray paths for the anomalous pattern might be far from those predicted in the preferred model. To overcome this limitation, the RRT method might be used (Zhao et al. 1992). Using the RRT method we inverted a synthetic data set which was generated for the preferred final model using the actual geometry. As noted by Zelt & Barton (1998)‘any feature of the preferred final model that is not present in the recovered model cannot be resolved’. Results of such an approach (Fig. 16b) indicate that the differences between the preferred and recovered solutions are of the order of 0.1 km s⁻¹. However, some small-scale anomalies within the SE part of the MB might not be well-resolved. For example, the difference between the two models at X = 325, Y = 125 km at 10 km depth amounts to 0.8 km s⁻¹.

Resolution tests provide only a rough estimate of the scale of the anomalies that can be recovered from the data. It is also desirable to determine the error bounds of each model parameter. The uncertainty of the ith model parameter can be described using the definition of an a posteriori covariance matrix (Tarantola 1987):

1

where

2

C_d, C_m are, respectively, data and model a priori covariance matrices, G is the Frechét derivative matrix, ε is the damping factor and is the Hessian inverse. In our case with ca. 300 000 model parameters, the calculation of Hessian inverse is impossible. Methods that do not involve Hessian inversion are still debated (Nolet et al. 1999; Yao et al. 1999). Thus we adopted the method proposed by Zhang & Toksöz (1998) and Korenaga et al. (2000). It is based on an assumption, that the a posteriori covariance matrix might be approximated as:

3

where denotes the mean of the K realizations (of equal probability). We assumed that K = 100 should be sufficient for a good estimate of . Following Korenaga et al. (2000) we prepared a set of 100 randomized 1-D starting velocity models (Fig. 17a). They were based on the 1-D average of our final model. The free parameters were: the depth of ‘consolidated’ basement and the V_P at 0, 5, 40, 50 and 60 km depths, weighted by the appropriate standard deviations of these parameters (see Table 5).

Figure 17.

(a) Plot of 100 randomized initial 1-D models. (b) Plot of the 1-D averages of 100 inverted Monte Carlo ensembles. Thick lines denote the 1-D average of the preferred model and ±3σ bounds. Insets show histogram plots of χ² values for the initial (a) and inverted set of models (b).

Table 5.

Parameters for choosing 1-D random models.

For the inversion we use also 100 randomized traveltime data sets with Gaussian noise (σ = 0.1 s) added. However, as noted by Zhang & Toksöz (1998), the traveltime picks are usually mutually correlated, thus their error distribution might not be properly described by a Gaussian.

In the next stage we ran 100 inversions using the ‘multiscale’ scheme. The mean initial χ² (see the inset in Fig. 17a) was 250, whereas the minimum value was 44, which is close to the χ² for our actual starting model (χ² = 43). For the inverted models the statistics are as follows: mean χ² = 10, minimum χ² = 6.5 and the maximum χ² = 18. Final χ² distribution is Gaussian-like (see the inset in Fig. 17b).

Slices through the mean model from all 100 realizations are presented in Fig. 18(a). In Fig. 18(b), we show the distribution of standard deviation for each model parameter based on the above quasi-Monte Carlo analysis. These results are quite surprising. The largest errors (apart from model boundaries) are close to 0.5 km s⁻¹ and occur at a depth of 5 km. They diminish with depth. This can be explained as due to the influence of ray density which reaches its maximum at the same depth (Fig. 13). Thus, the errors would arise from ‘accumulating’ the noise added to the data. We cannot however rule out that it is the mixed effect of the data noise and suspected anisotropy. Where the ray density is lower, the dispersion of the 1-D starting models is probably responsible for the observed σ distribution. The mean σ for the Monte Carlo ensembles tends gradually with depth to the mean σ used for creating the 1-D models (Table 5). Apart from the afore mentioned ambiguities, the analysis performed allows question the significance of some of the recovered velocity anomalies. For example, if we compare the anomalies marked by red arrows in Fig. 10 with the same locations in the mean model (Fig. 18a), it is apparent that they have much lower amplitude or even disappear. Their location is correlated with higher σ (Fig. 18b) in comparison to surrounding areas and thus we conclude that these anomalies are not likely essential features of the final model.

Figure 18.

(a) Horizontal slices through the mean model from 100 Monte Carlo ensembles. Velocity contour 0.25 km s⁻¹; (b) Horizontal slices through the standard deviation of 100 Monte Carlo ensembles. Velocity contours 0.1 km s⁻¹.

We observe relatively small dispersion of the inverted Monte Carlo models averaged to simple 1-D (Fig. 17b), despite large initial χ² values (inset Fig. 17a). It could be the result of the assumed σ-weighting which was quite narrow in the middle crust. Larger deviations from the mean, that occur in the shallow part of the model, might indicate the need for a more reliable velocity model for the sedimentary layers.

The mean error of our tomographic velocities was on the order of 0.1 km s⁻¹ (1σ). Due to the ambiguities mentioned, we are probably far from determining the ‘real’ errors of each model parameter using such a simple quasi-Monte Carlo analysis. We suspect that in the case of a 3-D model, 1-D starting models might not properly sample the model space and hence some more complicated models (e.g. 2-D) should be used.

4 Discussion and Preliminary Geological Interpretation

At this stage of interpretation, without integration of near surface geological observables, we can draw only some general conclusions concerning the geological importance of the obtained 3-D model of the whole crust. We attempt however to provide a preliminary interpretation of some of the highlights of our P-wave velocity model. The unresolved influence of the anisotropy on our results may mean that some of the anomalies located within the MB are not properly recovered. Based on the resolution tests, we may assume that the most reliable image of the crust and uppermost mantle structure is confined to the area between X = 100–350 km and Y = 100–300 (350) km. In order to illustrate the correspondence between the velocity pattern and tectonic units we present selected depth slices through the ‘multiscale’ model with a tectonic map overlayed (Fig. 19).

Figure 19.

Horizontal slices at 5, 10, 15 and 35 km depth with tectonic map overlayed. Velocity contours 0.25 km s⁻¹. Positions of vertical slices from Fig. 14 are indicated by yellow lines.

At 5 km depth (Fig. 19a) we observe two main areas of the crystalline basement with elevated V_P > 6 km s⁻¹. These are the Łuków Horst of the EEC and some granitoid massifs within the USB. Notably, the Holy Cross Mountains Fault (HCF) correlates well with the velocity gradient and transition from higher velocities (5.5 km s⁻¹) in the southern part of the Holy Cross Mountains (so-called Kielce unit) to lower velocities (4.75–5 km s⁻¹) in the Łysogóry part. Within the Carpathians we observe higher velocities in the western part (5.5 km s⁻¹) and lower (4.5 km s⁻¹) in the eastern part. The former might represent the basement whereas the latter are typical for the Carpathian flysch rocks. This would be consistent with the inferred deepening of the flysch strata towards the east (Ryłko & Tomaś 1995). Looking at the velocity gradient we deduce that the KLF Zone continues beneath the Carpathian nappes. Some small-scale anomalies within the TESZ correlate with the location of the Grójec Fault (GRF).

The boundary of the USB with the MB and the crystalline crust of the EEC is even more visible at 10 km depth (Fig. 19b). In the same slice there is a good reflection of the Lublin Trough. The axis of this unit is the EEC edge and it is bounded by the 6 km s⁻¹ isoline. West of the Lublin Trough we observe a velocity anomaly of 6.5 km s⁻¹ associated with the surface location of the Radom-Kraśnik High (Dadlez 2001) which is usually considered as a part of the LU. The picture of the EEC edge becomes more complicated to the south, but these anomalies are suspect in terms of their reliability. The positive velocity anomaly at the crossing of CEL01 and CEL02 profiles is not present in the 2-D forward ray tracing model (Środa et al. 2006; Malinowski et al. 2005). Crustal anisotropy might be the likely cause of its occurrence as discussed in Section 3.4. This area is also associated with large modelling errors as indicated by the Monte Carlo analysis (Fig. 18b).

At 15 km depth (Fig. 19c), the USB is visible as a compact area with velocities above 6 km s⁻¹. In the central part of the MB, we observe velocities below 6 km s⁻¹ reaching as deep as 20 km. Such low velocities were reported along the CEL02 profile (Malinowski et al. 2005) and interpreted as representing Neoproterozoic metasediments. Higher velocities observed in the eastern part of the LU suggest its EEC affinity and provides further support for distinguishing the Radom-Kraśnik High. A zone of HVB with V_P of ca. 6.5 km s⁻¹ is associated with the edge of the EEC. Those HVBs are extending down to ca. 20–25 km depth and their magnitude increases to 6.8–7 km s⁻¹. They were reported previously along the 2-D CELEBRATION 2000 profiles (Janik et al. 2005; Malinowski et al. 2005; Środa et al. 2006). These anomalies are also well tied to the observed Bouguer gravity anomalies (Bielik et al. 2006) forming the Lublin Gravity High (Grabowska et al. 1998). Based on the gravity modelling, they were interpreted as mafic intrusions (Grabowska & Perchuć 1985; Grabowska & Bojdys 2001). The velocity pattern again is more complicated in the SE part of the TESZ. We notice a kind of ‘triple junction’ with the velocity anomaly of 6.5 km s⁻¹ at the prolongation of the HCF which joins the zone of HVBs. Monte Carlo analysis suggests however that this anomaly is recovered with large errors (Fig. 18b).

At 30 km depth we may notice a difference in lower crustal velocities between the EEC (below 7 km s⁻¹) and the USB (above 7 km s⁻¹). In the western part of the MB, a zone with V_P = 7.5–7.6 km s⁻¹ is present, which should be associated rather with crust–mantle boundary. The thinning of the crust is more visible at 35 km depth (Fig. 19d). However, in the results of the joint refraction/reflection tomography using the P_mP phase (Malinowski 2006), the shallowest Moho occurs directly beneath the USB (32 km). Such a discrepancy might be explained by the fact that the P_n waves used in this study might not necessarily refract at the same interface as where the P_mP reflections originate (Zelt et al. 2003).

5 Conclusions

In this study, we presented the application of 3-D first-arrival tomography to the data recorded during the CELEBRATION 2000 experiment. The complex geological setting, as well as the large model volume, make the inversion challenging. Thus, we used two inversion approaches with the final conclusion that the preferred solution is obtained by gradually stepping from larger model cells to the smaller ones (‘multiscale’ inversion). As far as we know, the quasi-Monte Carlo analysis presented in this study, is the first such application in case of a large 3-D model based on deep seismic sounding data. This analysis allows us to present the approximate distribution of the velocity errors and hence to validate some of the recovered velocity anomalies. The highlights of our modelling results can be summarized as follows.

• (i)

  We may distinguish two crustal units with well-defined crystalline basement: USB and the EEC.

• (ii)

  Relatively low crustal velocities (V_P < 6 km s⁻¹) are observed down to 20 km depth in most of the MB and the LU. The only exception is the area of the Radom-Kraśnik High where higher velocities are reported.

• (iii)

  Since there are only slight differences in the velocity structure of the MB and LU, a ‘terrane’ subdivision based on the seismic properties alone is not viable.

• (iv)

  The Holy Cross Mountains Fault could be traced in the shallower crust (at 5 km depth), but its continuation at depth is not confirmed.

• (v)

  The USB continues beneath Carpathian Mountains with the KLF Zone as its eastern boundary.

• (vi)

  Several HVB were located along the edge of the EEC, which together with the observed low velocity infill of the TESZ led us to the conclusion that the EEC edge resembles an Atlantic type of passive margin (e.g. Holbrook et al. 1994). This seems to be in agreement with recent sedimentological data (Jaworowski & Sikorska 2006).

• (vii)

  The thinnest crust (with the ‘tomographic’ Moho at 30 km depth) is occurring between the USB and MB in a NW–SE trending zone.

The complex pattern of anomalies in the MB is likely resulting from anisotropy. It would thus be desirable to perform anisotropic tomographic modelling in the future. However, it would be faced with several problems, like choosing an efficient eikonal solver in an anisotropic medium. Also, the parametrization which properly handle the effects of anisotropy would make our inverse problem even more underdetermined.

Results of the joint refraction/reflection tomography using the P_mP reflections are promising (Malinowski 2006; Malinowski et al. 2008) and should substantiate our model with the detailed Moho shape. This would be the starting point for analysis combining other geophysical fields, like 3-D gravity modelling or thermal and rheological modelling. The ultimate goal is to integrate existing geological data with modelling results and produce a self-consistent geodynamic model for the complex area of SE Poland. In this sense, this study may be considered as the first essential step leading to such a model.

Acknowledgments

Funding for the CELEBRATION 2000 experiment has come from several institutions: in Poland, this project was supported by the State Committee for Scientific Research (KBN, grant PCZ 06/21/2000), the Ministry of the Environment, The National Fund for Environmental Protection and Water Management, the Polish Oil and Gas Company through the Association for Deep Geological Investigations in Poland (ADGIP). IRIS/PASSCAL is supported by the U. S. National Science Foundation (NSF). The instruments were provided through grants to the University of Texas at El Paso (State of Texas Higher Education Coordinating Board, NSF/MRI, and DoD). The public domain GMT package (Wessel & Smith 1995) was used to produce some of the figures. MM is supported by the Foundation for Polish Science through Young Scientist's Scholarship programme. Comments made by C. Ebinger and an anonymous reviewer helped us in improving our manuscript. Special thanks goes to D. White for language corrections.

Author notes

*

Now at: Geological Survey of Canada, 615 Booth St., Ottawa, ON K1A 0E9, Canada.

†

CELEBRATION 2000 Working Group comprises: P. Środa, W. Czuba, T. Janik, E. Gaczyński (1), G. R. Keller (School of Geology and Geophysics, University of Oklahoma, Norman, USA.)