Data fusion in geophysics: Seismic tomography and crustal structure in Poland as an example
Published:January 01, 2006
- PDF LinkChapter PDF
Matthew G. Averill, G. Randy Keller, Kate C. Miller, Piotr Sroda, Tiffni Bond, Aaron Velasco, 2006. "Data fusion in geophysics: Seismic tomography and crustal structure in Poland as an example", Geoinformatics: Data to Knowledge, A. Krishna Sinha
Download citation file:
The Trans-European Suture zone in eastern Europe is a complex region that marks the suture of Phanerozoic western European terranes with the Precambrian East European craton. Because of sedimentary cover and the Baltic Sea, the nature of this suture is known primarily from geophysical studies. Seismic velocity and gravity models from earlier experiments indicate changes of crustal thickness from 28–35 km to 42–47 km across the Teisseyre-Tornquist zone from the Paleozoic platform of western Europe to the East European craton. Tectonic models suggest the presence of a Precambrian–early Paleozoic passive margin beneath the Teisseyre-Tornquist zone. The Holy Cross Mountains in southeastern Poland represent an anomalous crustal block whose origin and interaction with the East European craton is unknown. We have employed a new approach to quantitatively integrate (i.e., data fusion) industry-acquired drilling and seismic reflection data, as well as refraction data from the CELEBRATION 2000, POLONAISE ‘97, and Teisseyre-Tornquist zone experiments in a tomographic inversion. We find that it is feasible to integrate data of different resolutions in a tomographic inversion. Secondly, we constrained upper-crustal velocities within our study area to create a more completely resolved velocity structure. Thirdly, we modeled a uniform gravity data set to validate all available data and construct models of crustal structure across the Holy Cross Mountains region.
The term data fusion has been used to describe many types of data integration efforts. In this study, we present a new approach for quantitative fusion of three fundamentally different types of seismic data in a joint tomographic inversion to determine the velocity variations in Earth in three dimensions. We applied our data fusion technique to a region in Poland where we had access to a variety of data sets. Our study area was part of the Trans-European Suture zone. This geologically complex and interesting feature extends 2000 km from the North Sea to the Black Sea (Fig. 1) and has evolved as the result of several phases of continental-scale deformation (Pharaoh, 1999). In our case, the inversion produced a volume in which seismic velocity is determined for blocks whose dimensions are on the order of a few kilometers, but smaller-scale investigations should work without difficulty. Our ultimate integrated analysis also employed gravity data and a variety of geologic observations. Thus, we will briefly illustrate the final geologic implications of the tomographic result that we obtained.
The data incorporated into this study include (1) results from a series of large-scale seismic experiments (Fig. 2) that targeted the deep structure of the Trans-European Suture zone: POLONAISE ‘97 (POlish Lithosphere ONsets—An International Seismic Experiment) (Guterch et al., 1999), CELEBRATION 2000 (Central European Lithospheric Experiment Based on Refraction) (Guterch et al., 2001), and the Teisseyre-Tornquist zone experiment (Grad et al., 1999), (2) industry reflection seismic profiles with stacking-velocity tables acquired from Apache Oil Company, (3) industry well log data including check-shot surveys, downhole velocity surveys, and stratigraphic columns acquired from the Apache Corporation, (4) gravity from the Polish onshore gravity database provided by Stanislaw Wybraniec at the Polish Geological Institute, and (5) geologic data from previous regional studies and the Geologic Atlas of Poland. The goals of this study are to integrate the available seismic (travel-time) data into a three-dimensional tomographic inversion to create a model for the velocity structure within this region and to provide a better understanding of the origin and tectonic evolution of this and adjacent crustal blocks in relation to the eastern Trans-European Suture zone and East European craton. In turn, the seismic data are correlated with gravity and geologic data to constrain structural models across the study area.
The Trans-European Suture zone is composed of a series of sutures created during the amalgamation of Phanerozoic European terranes with Proterozoic eastern European terranes (Guterch et al., 1999). This zone extends 2000 km from the North Sea to the Black Sea (Fig. 1) and has evolved as the result of Caledonian, Variscan, and Alpine composite orogenic cycles (Pharaoh, 1999). The Caledonian and Variscan orogenies resulted in accretion of Phanerozoic European terranes to the rifted southwestern margin of Baltica (East European craton) from 490 Ma to 290 Ma (Keller and Hatcher, 1999). Phanerozoic Europe consists of an unknown number of microcontinent-sized crustal blocks possibly derived either from tectonic transport around the margin of Baltica (having Baltica affinity) or from the Gondwanan supercontinent margin (Pharaoh, 1999). Terrane boundaries are generally poorly exposed oceanic sutures, primarily marked by geophysical lineaments and defined using faunal affinities (Pharaoh, 1999). The Teisseyre-Tornquist zone marks the eastern boundary and southeastern extent of the Trans-European Suture zone adjacent to Baltica (East European craton), as well as the coalescence of several tectonic features (Fig. 1). Trending northwest-southeast from the Baltic Sea to the Black Sea, the Teisseyre-Tornquist zone is also thought of as the contact zone between the eastern and western European platforms (Lamarche et al., 1999). This zone and its adjacent terranes are marked by distinct upper-crustal features as well as deeper structure constrained by geophysical techniques such as refraction seismic, gravity, magnetic, and heat flow (e.g., Guterch et al., 1999).
Large-Scale Seismic Experiments
The Teisseyre-Tornquist zone experiment targeted a key portion of the Trans-European Suture zone adjacent to the axial part of the Polish trough (Grad et al., 1999) (Fig. 2). Along this profile, a thick (up to 20 km) layer of sediments and metasediments was found, and the thickness of Earth's crust increases from northwest to southeast from 35 to 41 km. In turn, profiles from the POLONAISE ‘97 project were designed to parallel and cross this previous experiment to further constrain the crustal structure southeastward along the Trans-European Suture zone.
In May of 1997, the POLONAISE ‘97 (Guterch et al., 1999) seismic refraction/wide angle reflection experiment took place predominantly in Poland and included extensions of the P4 profile into Germany and Lithuania. This experiment utilized 800 recording instruments to record along five profiles (Fig. 2). The layout of the four main interlocking profiles was designed to produce significant three-dimensional coverage, and tomographic analysis has been undertaken on these data in both two (Guterch et al., 1999) and three dimensions (Sroda et al., 2002). These results indicate the presence of a deep, up to 20 km, asymmetric basin consisting of sedimentary to metasedimentary or volcanic rocks with velocities of <6.1–6.2 km/s overlaying a higher-velocity (6.7–7.5 km/s) lower crust with the base varying significantly in depth across the Trans-European Suture zone.
The goals of CELEBRATION 2000 (Guterch et al., 2001) were to (1) investigate the deep structure of the southwestern margin of the stable continental interior of Europe, (2) delineate the major terranes that compose western Europe in the region, and (3) construct a three-dimensional model of the lithospheric structure and evaluate and develop geodynamic models for the evolution of the area. This experiment produced a very large data set by recording along a series of interlocking profiles whose total length was ∼8900 km. These profiles were designed to form a three-dimensional array to augment the 5400 km of traditional profiles (Fig. 2). The analysis of data from this experiment and fusion of these data with those from other sources was the primary motivation of this study.
Seismic tomography is a relatively new approach to determine the velocity structure of the earth (e.g., Nolet, 1987). In this approach, velocity anomalies are detected by analyzing the travel times of many seismic rays through a region of the earth. Usually some reference earth model (velocity structure) is assumed, and variations from this model are determined and typically expressed as percent change in velocity and displayed in color. This approach is similar to that used in medicine to form images from X-rays or ultrasound, such as the CAT scan (computerized axial tomography), and is substantially different from the traditional approaches used in seismic reflection and refraction studies (e.g., Yilmaz, 2001). The situation in the earth is more complex because seismic rays bend as they travel and because the distribution of seismic recorders and seismic sources is not uniform. As in most other geophysical techniques, tomography can be applied on either a local or global scale. Thus, seismic tomography is an active area of research that is producing interesting new images of the velocity variations in the earth at a variety of scales.
The essence of seismic tomography is to use a large number of observations for waves that have traveled a variety of paths in order to mathematically determine a velocity model. The most commonly employed approach is to approximate the earth by a three-dimensional mesh of blocks (e.g., Hole, 1992). The travel-time integral above is used to sum the travel times through each block intersected by a raypath. Many observations then lead to a large system of linear equations where the unknowns are the velocities of the blocks. The quality of the solution depends on the number of observations and their accuracy, the distribution of raypaths, and the degree to which the block model can actually approximate the true velocity structure.
In this study, we were in the unique position of having access to a considerable amount of additional data in the form of borehole check-shot surveys and stacking velocities derived from the processing of seismic reflection data. Thus, we devised a method to include these data in the tomography in order to obtain better constraints on the velocity structure of the upper 20 km of crust in the study area.
Seismic Tomography Algorithm
The tomographic approach that was used for this study was a three-dimensional seismic travel-time tomography or back-projection method for determining three-dimensional velocity structure from first-arrival travel-time data (Hole, 1992). It uses a linear tomographic inversion algorithm that minimizes computational time, provides a relatively high-resolution model, and can handle structures with large velocity contrasts (Hole, 1992). It is roughly a two-step process: (1) First-arrival travel times are calculated for each node in the grid by solving the eikonal equation (an approximation to the wave equation) with a finite-difference method (Vidale, 1988), and (2) travel-time residuals are inverted for changes to the velocity model. As a starting model, we used a three-dimensional expansion of a simple one-dimensional initial velocity distribution for the region that was based on known geology and previous seismic studies. This model is then expanded to a uniform grid with equal spacing in all three dimensions.
Once travel times have been calculated through the initial model, the second or inverse step begins. The model is parameterized in terms of constant slowness (reciprocal of velocity) cells. The goal of this step is to invert the differences between observed and calculated travel times (residuals) for changes in slowness that minimize the travel-time residuals. Next, the slowness of each cell in the initial model is determined by the average of the slowness of the cell's eight bounding nodes (Zelt and Barton, 1998). The computations produce two sums for each cell that are needed in order to update the slowness. The first is the raypath-averaged slowness perturbation for each ray passing through the cell, which equals the travel-time residual (observed minus calculated) of the ray divided by the ray's total length through the cell. The second sum is equal to the number of rays passing through each cell. Once these calculations have been performed for each ray that traverses the model space, the first sum is divided by the second to obtain the slowness update for each cell through which at least one ray passes. The update is then added to the model and the whole process starts again once a moving average filter is applied to stabilize or smooth the updated model. The moving average filter is input into the program as a starting parameter and works most efficiently when its size is decreased after each run. For a more complete discussion of the computational method, see Hole (1992). A typical run consists of four to seven iterations. As a result, models subject to smaller smoothers define progressively smaller wavelength features (Zelt and Barton, 1998). Since it is not reasonable to update the model until the travel-time residual equals zero because of picking errors, a model is considered to be acceptable when the RMS (root-mean-square) of the travel-time residuals is roughly equal to the picking error.
Description of Travel-Time Data Sets
In order to define the boundaries of our velocity model, the distribution of shots, receivers, check-shot surveys, and stacking-velocity determinations were plotted relative to the tectonic features of interest (Fig. 3). Based on the compilation of all these data, a model space was constructed to maximize the data coverage with respect to the Trans-European Suture zone and Polish trough (rectangle in Fig. 3).
The refraction data set consisted of observations of seismic arrival times for receivers and sources that lie within the bounds of the model space. This large data set included 500 arrival times from the Teisseyre-Tornquist zone experiment (from six shots and 155 receivers), 1520 arrival times from the POLONAISE ‘97 experiment (from 16 shots and 253 receivers), and 5806 arrival times from the CELEBRATION 2000 experiment (from 33 shots and 616 receivers). Since the Teisseyre-Tornquist zone and POLONAISE ‘97 data within our model space have been analyzed in several previous studies (Grad et al., 1999; Guterch et al., 1999; Jensen et al., 1999; Sroda et al., 1999, 2002), the first-arrival times (picks) were well established. As a result, these data were obtained directly from the Institute of Geophysics of the Polish Academy of Sciences in a format appropriate for input to the inversion code. By contrast, the CELEBRATION 2000 data set has only been analyzed along individual profiles, and the results are preliminary. Therefore, for consistency in this study, we picked first arrivals both for individual profiles and for fan shots within the established model space, and formatted them for use in the inversion. An example of the seismograms for a profile is shown in Figure 4. Dashed lines mark the approximate location of the first arrivals that were input into the inversion code on this figure.
The “nonrefraction” data set is a compilation of well data and velocity data from reflection profiles acquired from the Apache Corporation. The well data were in the form of check-shot surveys in which a seismometer is lowered into the well and records a source that is located near the well. This technique provides a very accurate measurement of the velocity distribution in the vicinity of the well. However, this type of survey is expensive and thus relatively rare. The processing of seismic reflection data involves the summing of seismograms from shot-receiver pairs that share the same approximate subsurface reflection point. This process is called common depth point (CDP) stacking. Before the stacking can be done, the seismograms must be aligned so that the reflections will add in a constructive fashion at the expense of random noise. This alignment requires the estimation of seismic velocity as a function of depth called stacking velocities (e.g., Yilmaz, 2001).
Check-shot surveys from individual wells consisted of either depth-time pairs or depth-velocity pairs for the position within the well. For purposes of the inversion, these data were converted to a one-way travel time from a shotpoint at the surface to a receiver at depth. The geometry for each well location was prepared by assigning a shotpoint number to the well. In addition, a receiver number was assigned to each depth for which there was a travel time. The travel-time picks available were for very short depth intervals (30–50 m), given the large spatial scale of our study. Thus, these data were decimated to consist of picks at only 2 km intervals in order to avoid numerical problems within the 2 km model cells.
The travel-time data that we determined from the processing of seismic reflection data consisted of travel times derived from the tabulation of stacking velocities available at intervals along the seismic profiles. These data were provided by the staff of Geophysika Torun (2003, personal commun.). Although it would have been ideal to completely reprocess the reflection data and thus have direct knowledge of the reliability of these velocity determinations, the original data were not available. Since stacking-velocity data take the form of time-velocity pairs, they had to be converted to depth-time pairs at each CDP location for input to the inversion. Once this was accomplished, each CDP location was assigned a shotpoint number. Each depth-time pair at each location was then assigned a receiver number in a fashion similar to the approach taken with the check-shot data. As with the well data, travel-time picks derived from stacking velocities were decimated to 2 km intervals. In addition, CDP data were cut off at 2.5 s one-way time due to the decreasing resolution of stacking velocities with increasing depth.
The total “nonrefraction” data set consisted of 213 arrival times from check-shot surveys conducted in 72 wells, and 525 arrival times derived from stacking-velocity determinations made at 234 CDP locations (Fig. 3).
Implementation of the Tomography Algorithm
The geographical information for all of the shot points and receivers of each dataset (Teisseyre-Tornquistzone, POLONAISE ‘97, CELEBRATION 2000, wells, and CDPs) was transformed into the X and Y coordinate system of the model space via shell scripts. These scripts are used to produce a series of output files that contain receiver offsets and travel-time picks derived from interpretative first-arrival picking, receiver locations, and shot-point locations (X, Y, Z). In our case, the receiver offsets for the well and CDP data are only in the Z dimension, with the surface location treated as the source. Therefore, velocities are only calculated along vertical raypaths below the source. These files are then used as input into the inversion.
The main parameters set during the inversion process include the model dimensions, cell size, and the iteration schedule. The model dimensions were set to 440 km along the X axis (perpendicular to the axis of the Trans-European Suture zone), 240 km along the Y axis (the axis of the Polish trough), and 70 km along the Z axis (depth). The cell size within the model was set to 2 km for each dimension. The iteration schedule was set to iterate six times at each of five moving average smoothing windows decreasing in dimension and incorporating all data for receivers within 350 km of the corresponding seismic source. The smoothing windows were set by the number of cells in the X, Y, and Z dimensions at 64 × 64 × 16, 32 × 32 × 8, 16 × 16 × 4, 8 × 8 × 4, and 8 × 8 × 2.
Results of Tomographic Inversions
Several steps were used to develop three-dimensional velocity models for the study area. First we derived a number of one-dimensional velocity models to be used as starting models for the inversion, which took into consideration previous models for the area as well as the new CELEBRATION data. These one-dimensional models were tested by inverting the refraction data only. This effort provided information on the dependence of the inversion result on the initial model as well as control models for the refraction data that were used to observe the effects of adding the industry CDP and well data. Based on these tests, a model with a smooth increase in velocity with depth was chosen. This model has a generally smooth velocity gradient that steepens at ∼22 km depth, followed by a gradual decrease in the velocity gradient to 62 km. A negative velocity gradient at the base of the model is used to prevent rays from escaping at its base.
An RMS error was calculated and reported for each shot-point based on travel-time residuals—differences between calculated and observed travel times—for first-arrival picks. The overall RMS error for a given iteration is the sum of the RMS for each shotpoint scaled by the number of picks for the shot. Initial RMS values near 0.75 converge to final values of 0.14 s for the refraction data, 0.14 s for refraction data with well data included, and 0.18 s for all the data. In order to monitor the inversion, two-dimensional vertical slices through the three-dimensional velocity model and horizontal slices at various depth intervals were examined for individual runs. The vertical slices were chosen based on the data density and geologic considerations and correspond to profiles at X = 300 km, which is perpendicular to the axis of the Trans-European Suture zone and crosses the Holy Cross Mountains and Lublin Trough area, and at Y = 100 km, which follows the axis of the Polish trough (Fig. 3).
Results Using Only the Refraction Data
The first step in the process was to undertake a conventional inversion of first-arrival picks for the combined Teisseyre-Tornquist zone, POLONAISE ‘97, and CELEBRATION 2000 refraction experiments. After some experimentation, a starting model with a smooth increase of velocity was employed because it produced the most stable inversion. A vertical slice (X = 300 km) of the result of this inversion is shown in Figure 5. This slice crosses the Teisseyre-Tornquist zone and, in the upper 20 km of the crust, shows distinct changes along the profile. The spatial distribution of these velocity variations is shown in a horizontal slice through the model at a depth of 15 km (Fig. 6). These views of the velocity models show velocities between 5.0 and 6.0 km/s forming a layer that quickly thickens into a deep (∼20 km) basin from the southwest (at Y = 0) toward the northeast into the Trans-European Suture zone. In Figure 5, this basin thins northeastward and develops a secondary shallower basin (∼12 km) before thinning to less than 10 km at the northeastern end of the profile (at Y = 240 km). In turn, a thin (0–5 km) cover of velocities less than 5 km/s develops midway along the profile and then thins to the northeast. The interval with velocities between 6.0 and 6.8 km/s is clearly thickest (10–20 km) to the southwest and dramatically thins to the northeast along the profile.
Results with Industry Data Included
Industry well and reflection data sets were independently added into the inversion with the refraction data. Well data, consisting of check-shot surveys, were added to the refraction data set, and the inversion procedure was repeated. In turn, CDP data, based on industry reflection lines, were added to the refraction and well data, and the inversion procedure was repeated again. Results for inversions including CDP data showed instabilities, and model RMS values did not decrease below 0.18 s, compared to values of 0.14 s for both refraction data and the combined refraction and well data. Although this result is disappointing, we attribute some of the errors to inherent difficulties in calculating stacking velocities from reflection seismic data. We suspect that careful velocity analysis, rather than the production mode approach used when all that is sought is a good stacking result, would have given us a good result, but we cannot quantify this assertion. Another consideration is that the CDP data were very closely spaced and thus could contain travel-time anomalies from local geologic features that were too small to be resolved by the much more widely spaced refraction and check-shot data (Fig. 3). We calculated the differences between the different inversion results and found that some groupings of CDPs with high RMS values showed strong correlation with areas where drilling and seismic reflection data demonstrate the presence of large bodies of salt, which have high seismic velocities. This indicates that some of the problems encountered by including the CDP data are not due to erroneous data but are due to our attempt to fuse data with spatial resolutions that are too different to be compatible.
The inclusion of the check-shot surveys, however, improved our inversion results. These well data, limited to the uppermost crust (<6 km), provide valuable and reliable velocity information for sedimentary rocks within the Polish trough. Plots showing the differences between the velocity models for the refraction data alone and for the refraction and well data combined show the effects of adding the well data (Fig. 7). In this figure, positive values indicate a decrease in velocity when the well data are added. This corresponds to better resolution of the low-velocity sediments in the upper crust (<5 km). In turn, velocities at deeper levels are affected by changes in the model in the uppermost crust, and velocity gradients are steepened beneath the basins. The difference plots helped delineate low-velocity regions as possible basin structures and demonstrated the positive effects of constraining the upper crust in addition to refraction data.
Vertical slices (X = 300 km; Y = 110 km) through our final velocity model resulting from the fusion of the refraction and well data are shown in Figure 8. These diagrams show that while this study includes a very a large amount of data, the extensive size of the study area provides only limited three-dimensional tomographic ray coverage. As a way of displaying only those areas of the model with good resolution, areas lacking “ray coverage,” i.e., that are not traversed by raypaths, were set to a null value, “masked,” in this model. As discussed in Averill (2003), this model has many interesting geologic implications that are beyond the scope of this paper. However, we will show how it can be qualitatively fused with a totally different type of geophysical data in the section below.
GRAVITY DATA AND MODELING
Measurements of Earth's gravitational field are a very cost-effective source of information about the structure of the lithosphere. These very precise measurements are processed so that known variations in the gravity due to factors such as elevation and latitude are removed. The remaining anomalies are the manifestation of geologic structures and are mapped and modeled to determine density variations. For many years, the scientific community has cooperated to compile gravity measurements into databases on a global and regional basis (e.g., http://paces.geo.utep.edu). The database for Poland is particularly good, and the gravity map of Poland (Fig. 9) contains several major anomalies of interest (Wybraniec, 1999) that can be compared with our velocity models. These include an anomalous block of high gravity values in the southeastern portion of our study area and the northwest-southeast-trending series of highs and lows that extend across central Poland. In fact, examination of the gravity anomalies helped to determine the bounds for our tomographic analysis. In unison with tomographic modeling, our goal was to jointly constrain the crustal structure of this region with gravity data. The gravity data used for this study were provided by the Polish Geological Institute (S. Wybraniec, 2003, personal commun.) in the form of a 5 km spaced gridded data set that was reformatted for input into our gravity modeling and mapping programs.
The gravity anomalies clearly delineate several tectonic features in the region (Fig. 9). The Holy Cross Mountains in the southeastern region of Poland as well as the Lublin Trough clearly coincide with gravity highs, shown by the hotter colors. A northwest-southeast–trending gravity low corresponds to the deep (>10 km) sedimentary deposits of the Polish trough. The Sudetes region in southwestern Poland is the core of a crystalline uplift in the Variscan terrane in this region and is clearly marked by a gravity high. Features such as the Odra fault zone and the Holy Cross dislocation (a fault that is a major tectonic boundary in the Holy Cross Mountains [Valverde-Vaquero et al., 2000]) are marked by lineations in the gravity anomalies and truncation of gravity highs. The Grojec fault (a subsurface feature that trends northeast across the study area at approximately x = 230 km; Fig. 9) marks the northwestern boundary of the Holy Cross Mountains block (the large gravity high that occupies the southeastern portion of the study area) and is clearly a major feature associated with the sharp truncation of the gravity high in that region.
By overlaying two-dimensional velocity depth slices onto our gravity map, we were able to relate velocity structures to the observed gravity anomalies. For example, at a depth of 15 km (Fig. 9), a northwest-southeast-trending low-velocity (<6.0 km/s) zone correlates with the trend of the Polish trough that is marked by low gravity anomaly values and represents an accumulation of sedimentary strata. A high-velocity (6.4–7.0 km/s) anomaly located over the Lublin Trough (northeastern half of the Holy Cross Mountains block) clearly accounts for the high gravity anomaly and suggests the presence of dense mafic material. The steep velocity gradient bounding the Holy Cross Mountains block on the northwest correlates with the Grojec fault, showing that it is a major, deep-seated boundary.
Gravity modeling requires plotting a profile of gravity values against distance, and a profile was constructed along the line shown in Figure 3. During modeling, the goal is to match this gravity profile of observed values with a calculated gravity profile produced by forward modeling of crustal features. Inverse methods can also be used to model gravity anomalies, but we chose forward modeling order to be able to integrate a wide variety of constraints. As this is a nonunique process, it is necessary to start with a well-constrained model, given known geologic relationships along the profile and independent geophysical constraints.
To constrain our models, we included locations of geologic and tectonic features, well and seismic data provided by industry, interpreted seismic refraction lines, and slices from tomographic modeling discussed above. The first step was to correlate available data with their location along the strike of the profile line. In turn, these data were used to create an initial density model that honors the constraints at any designated location. This initial model is then adjusted, honoring the constraints to best fit the observed profile. The result is an integrated structural model along the strike of the profile that jointly satisfies the seismic, gravity, and geologic data available.
In order to prepare the model, crustal structure is represented by individual polygons within the model. These polygons are defined by nodes in model space and are assigned a density value based on constraints such as velocity and lithology from wells, seismic surveys, and known velocity-density relationships. Our computer code, which is based on the method of Cady (1980), uses the dimensions and density of these polygons to produce a calculated profile for the model structure. Due to the integrated process used to derive these models, they should be considered as crustal-scale cross sections, as opposed to crustal models derived from gravity anomalies alone.
The gravity model across the Holy Cross Mountains (Fig. 10) clearly shows a transition from a Paleozoic platform terrane (Upper Silesian block) to the East European craton. At this transition, the thin (31–33 km) crust of the Paleozoic terrane is in juxtaposition with the three-layer, 43–45 km thick crust of the East European craton. A large mafic body marks the edge of the rifted cratonal margin and accounts for the large gravity high over the ∼15 km deep Paleozoic basin within the Holy Cross Mountains block. A thickened portion of the mafic body underlies and accounts for the Lublin Trough gravity high. A thickened carbonate Paleozoic section (10–15 km) within the transition zone underlies the Holy Cross Mountains region. These features are suggestive of a volcanic continental margin formed during the breakup of Rodinia. These bodies were modeled with lateral variation in and out of the model, corresponding to the limited lateral extent of the associated gravity anomalies. The wedge-shaped feature modeled within the upper crust of the East European craton margin is interpreted to correspond to deformation of the margin. Deformation of the Paleozoic section along the East European craton margin corresponds to uplift of the Holy Cross Mountains. The lack of large-scale orogenic features such as major basement deformation and uplift suggests an incomplete or indirect collision during the Variscan orogeny.
The lateral extent of some of the features is shown in Figure 9, where a horizontal slice of the velocity model is shown along with the gravity anomaly map. A region of high velocities (∼7.0 km/s) coincides with the thickened mafic body underlying the Lublin Trough and marks the northeastern edge of the Teisseyre-Tornquist zone. A corresponding velocity high to the southwest corresponds to the southwestern edge of the transition zone and the edge of the Paleozoic European platform. A trough of velocities less than 6.0 km/s corresponds to the thickened sedimentary strata. The 6.0 km/s contour roughly coincides with the top of the crystalline basement across the transition, and the Grojec fault shows up as a discontinuity in the velocity structure at 15 km, confirming that it is a major crustal feature. The result of this integrated modeling thus sheds additional light on the deep structure of this region and represents a form of data fusion in which a two-dimensional gravity model is used as the medium for the integration of a wide variety of data.
The results of this study include a successful demonstration of the benefits and drawbacks of quantitative fusion of small-scale industry-acquired well and reflection data with large-scale refraction data as well as an integrated interpretation of the structure and evolution of a portion of the Trans-European Suture zone and the surrounding region in southeastern Poland. Using the results of this study and previous work we can conclude the following:
The integration of travel times from small-scale data sets, such as well and reflection data, into a large-scale tomographic inversion can provide additional constraints on the structure of the uppermost crust and in turn more confidence about the resolution of deeper structure. Our approach has the advantage of including all travel-time data in one inversion.
The presence of a mafic lower crust, indicated by a thickened layer of high-velocity (6.8–7.4 km/s) material and a distinct high gravity anomaly is consistent with observations from modern rifted passive continental margins (e.g., Lizarralde and Holbrook, 1997).
In the Holy Cross Mountains block, the presence of a high-velocity zone (6.6–7.2 km/s), modeled as a high-density (3.1 g/cm3) mafic body, accounts for the observed gravity anomaly (greater than +30 mgal) located beneath an ∼15 km deep basin with velocities <6.0 km/s. These velocities are consistent with Paleozoic sedimentary rocks of East European craton margin affinity. In turn, they are part of the Holy Cross Mountains block.
Across the Trans-European Suture zone, there is a transition from a thick (>40 km), three-layer crust to a thinner (31–33 km) crust that is consistent with previous work.
Lack of large-scale orogenic features across the transition and the observed deformation of Paleozoic marginal sediments in the Holy Cross Mountains suggest indirect or incomplete Variscan collision/suturing along the East European craton margin.
Our integrated gravity modeling is a quantitative but less formal type of data fusion in that the various constraints are used in a qualitative fashion. We are presently investigating ways in which other types of seismic data can be used in a mathematically formal fusion via tomographic inversion.
We would like to particularly acknowledge the assistance of John Hole, who helped us understand his tomography program and provided a helpful review. This work is part of an ongoing collaboration between the University of Texas, El Paso, and the Institutes of Geophysics of the Polish Academy of Sciences and the University of Warsaw. Apache Petroleum provided most of the industry data used in this study. Our participation in the seismic data collection was funded by grants from the National Science Foundation (NSF Celebration: OISE-0001356) and the State of Texas Higher Education Coordinating Board. Our research on data fusion is funded by the GEON project (NSF grant EAR-0225670) and by NASA through PACES (cooperative agreement NCC5–209). IRIS/PASSCAL is supported by the National Science Foundation (NSF) and provided the majority (∼700 of 1200) of the instruments for this experiment, and most of these instruments were purchased through grants to the University of Texas, El Paso (State of Texas Higher Education Coordinating Board, NSF/MRI, and U.S. Department of Defense).
Figures & Tables
Geoinformatics: Data to Knowledge
- Central Europe
- data processing
- digital data
- geophysical methods
- geophysical profiles
- geophysical surveys
- gravity anomalies
- gravity methods
- inverse problem
- reflection methods
- refraction methods
- Russian Platform
- seismic methods
- seismic profiles
- suture zones
- Swiety Krzyz Mountains
- three-dimensional models
- thrust faults
- Tornquist-Teisseyre Zone
- velocity structure