Evaluation of the Possibilities Validation of Interval Velocity Models Using Non-Seismic Data and Its Impact on Geological Interpretation of PreSDM Results

Abstract

The paper presents the problem of generation and validation of Velocity Interval Depth (VID) models with the application of non-seismic geophysical and geological data. The study area is a part of the Carpathian Foredeep located close to its contact with the Carpathian Overthrust. In this area of complicated geological structure, hydrocarbon deposits have been successfully explored for decades with seismic methods and drilling. The research applied the Simultaneous Joint Inversion (SJI) of independent geophysical data, which is a modern methodology of geophysical data processing, that is still under development. Such an attempt was necessary due to the lack of a sufficiently dense grid of wells in the study area, in which seismic velocities would be correctly recorded. Such data would be then applied for the generation of relevant VID models, which in turn, could be used to perform the Prestack Depth Migration (PreSDM) procedures. The application of procedures taking advantage of independent geophysical and geological data enabled researchers to control the generation process of the spatial VID model in the areas without wells. The analyses aimed to verify the correctness of VID model evaluation and its influence on the quality of seismic imaging in the area of the Carpathian Overthrust. Precisely, the influence was tested of such non-standard generation procedure of seismic velocity fields, not only on the PreSDM results but also on the geological interpretation of both the Rączyna and the Jodłówka gas deposits. The latter aspect of the presented results seems to be crucial to the effectiveness of petroleum exploration in the transition zone between the Carpathian Orogen and the Carpathian Foredeep.

1. Introduction

The paper presents the problem of construction and validation of seismic wave Velocity Interval Depth (VID) models with the application of geophysical and geological data. One of the validation methods of VID models is the evaluation of spatial positioning errors of selected seismic horizons after the application of the PreSDM procedure in relation to the interpretation of well log datasets. This is the standard Quality Control (QC) procedure executed during the processing of reflection seismic data. It is correct in the areas where the wells have already been completed and the relevant well logs have been recorded. Additional applicability criteria of this validation method are depth interval and quality of recorded geophysical data. Unfortunately, the available geophysical datasets commonly do not cover the full depth interval. The lack of suitable well logs contributes to attempts to validate the VID models with the independent geophysical data [1,2]. Taking into account the professional experience of the authors and the available publications [1,2,3], it is assumed that in the areas where the well-log markers are unavailable, it is necessary to use the processing results of data acquired with the geophysical methods other than reflection seismic, i.e., gravity and magnetotelluric. The geophysical method commonly used as complementary for seismic surveys is gravity. Its application is justified by close mathematical relationships of seismic wave velocity, bulk density of rock formations and their elastic parameters [4,5]. In addition to gravity, the processing and interpretation of seismic data can be supported by the results of magnetotelluric surveys, which ensure the relatively great depth of geoelectric recognition of rock formations [6,7] as well as a depth range independent of the distribution of seismic waves velocities. Unfortunately, the close physical relationship between the seismic wave velocity and the electric resistivity does not exist [4,5,8]. Thus, if the methods based on the resistivity distribution within the rock formations must be applied, the empirical statistical relationships are in use [4,5,9]. In practice, the determination of correlations between resistivity and seismic wave velocity is possible for genetically defined, relatively homogenous rock formations but, unfortunately, it often is rather ambiguous. Modern statistical methods are quite helpful [4] together with machine learning methodologies [10,11,12,13]. In the project reported below, the authors applied the modern method of integrated inversion of independent geophysical data, which is still under development.

The definition of this method is quite broad and is often equated with the SJI, i.e., one of the options of the Cooperative Joint Inversion (CJI) methodology [9,14]. Another variant of this method is the CJI with the application of structural (or parametric) bonds, which replace the not-always correlations between various physical parameters (velocity, density and resistivity) [3,9]. In the Polish literature, this CJI variant is known as the integrated inversion (interpretation) with constraints (Integrated Inversion Interpretation) [6,7,15]. In this variant, more emphasis is placed on the structural image than on ensuring the correlation of physical parameters at a satisfactory determination level R². Due to the progress in the development of these methodologies and access to relevant software tools, this CJI variant has become the additional validation method of velocity field determinations. The SJI method provides automatic and simultaneous computation of inversion for two or more sets of geophysical data, which are completely independent of each other [1,2,3,9]. It must be emphasized that the procedure described in this paper has already been the subject of interest by leading working groups in the world involved in the development of this methodology [1,2,3,9,14]. Also, the research teams from Poland focused on the following subjects [15,16,17,18,19,20]. The results presented below are the continuation of years-long research focused on the application of independent data used to construct seismic wave velocity fields [1,2,21,22,23]. In addition to a number of publications dealing with this problem, attention must be paid to the solution proposed by the working group from the Polish Oil and Gas Co. (POGC) for seismic, gravimetric and magnetotelluric datasets and applied to the comprehensive interpretation of regional seismic data acquired within the project of Transcarpathian Seismic Profile (POLCRUST) [7,24,25,26]. For that solution, several analyses and tests were performed with the application of various seismic wave velocity fields [7,25,26]. Particularly, the modeling of full wavefield was carried out in order to evaluate the possible advantages provided by the application of such inversion to seismic data processing [26]. The analyses presented below aimed to test the influence of standard generation procedures of velocity fields not only on the PreSDM results but more importantly on the effects of geological interpretation of the Rączyna gas deposit area [27]. The latter aspect of the presented research seems to be crucial to the effectiveness of petroleum exploration in the transitional zone between the Carpathian orogen and the Miocene Carpathian Foredeep. The research included both the geophysical and geological datasets originating from the study areas of the following project: “The new 3D seismic-geological model testing the perspectives of discovery of hydrocarbon unconventional deposits within the Miocene formations of the Carpathian Foredeep” performed within the INGA 2 research project. Both were co-financed by POGC and the National Centre for Research and Development.

The article has been divided into six chapters. In Section 1, “Introduction”, a general description of the purpose of this material was provided. Section 2, “Location of Research Area”, includes information about the location and geological characteristics in relation to archival seismic surveys conducted in the studied area. Section 3, “Methods and Materials”, is a crucial part of this work as it describes the applied methodology for validating independent velocity models used in the seismic data processing. Additionally, the seismic and gravimetric data used for calculations are described in this chapter. To facilitate a better understanding of the subject matter, Section 3 has been divided into five subsections. Section 4, “Research Results and Discussion”, presents the obtained results along with a discussion on the possibilities offered by the application of independent velocity models in the seismic data processing. The last section, “Summary and Conclusions”, includes a summary of the described procedure presented by the authors. Both the benefits and limitations of this methodology were indicated.

2. Location of Research Area

The study area is located in quite a well-recognized region of the Rączyna and the Jodłówka gas deposits [27], south of Jarosław city, in the Fore-Carpathian District. The PreSDM procedure was tested along a fragment of the POLCRUST Trans-Carpathian Seismic Profile [7,24,25]. Datasets acquired from this profile were reprocessed within the two research projects co-financed by the POGC and the National Centre for Research and Development. The INGA 2 project aimed to improve seismic imaging in the areas of the Carpathian Overthrust and the Carpathian Foredeep (M-5 traverse). The precise location of the study area is depicted in Figure 1 in connection with the structural map of the Miocene formations’ basement, which was interpreted based on both 2D and 3D seismic sections [27]. Additionally, the map contains the localization of the most important wells, which provided data used in analyses as well as the ranges of reservoir rocks in both the Rączyna and the Jodłówka areas. The seismic data processing was focused on analyzing the Middle Miocene autochthonous sediments (Figure 2A,B—No. “2”), which correspond with the described Rączyna and Jodłówka gas deposits [27].

The geological structure of the study area was illustrated through the interpretation of two archival seismic profiles crossing the Miocene autochthonous sediments, which are natural reservoirs for the Rączyna and Jodłówka gas deposits. The autochthonous Miocene formations directly cover the highly varied, erosional top surface of the Precambrian basement immediately beneath the Carpathian-Stebnik Overthrust. The development of sedimentary formations infilling the Foredeep was strongly influenced by the paleorelief of the Precambrian basement top surface (Figure 2A,B—No. “3”) dominated by deeply incised river paleovalleys. Similar to analyzed Miocene formations, the geometry of the sub-Miocene basement was finally shaped by thrusting movements of both the Flysch Carpathian nappes and the allochthonous Stebnik Unit (Figure 2A,B—No. “1”) [27,28]. Most part of the study area was covered by the sediments involved in the Carpathian-Stebnik Overthrust (Figure 2A,B—No. “1”). The Precambrian formations extended throughout the study area provided the basement of deposition of Miocene formations, which infill the Carpathian Foredeep [28]. These three structural stages are temporarily and spatially interrelated due to the dynamics of the development of the Carpathian Foredeep eastern segment [27,29]. In the study area, the Precambrian is represented by low-grade-metamorphosed (phyllitized) clay tones and mudstones sporadically intercalated by quartzites and quartzitic sandstones. All these sediments show high dip angles, which indicates strong tectonic activity whereas the erosional, morphologically varied Precambrian top surface dips generally to the south [27,28,29]. This top surface is directly covered by autochthonous Miocene formations (Figure 2A,B—No. “2”), which leveled the relief. The Miocene formations represent the Lower Badenian (Baranów Beds), the Upper Badenian (including the evaporites) and the Sarmatian (the youngest unit of the autochthonous Miocene succession) [27,28]. The formations comprise clayey-silty shales with fine-grained sandstones. Less common are fine-grained sandstones interbedded by shales.

The Carpathian Overthrust (Figure 2A,B—No. “1”) involved sediments of various ages but, due to the similarity of their tectonic styles, all are presented together. The Stebnik Unit, which underlies the sediments involved in the Carpathian Overthrust structure, comprises shales and mudstones intercalated by sandstones and minor conglomerates. In the bottom parts of this succession, thicker horizons of anhydrites, gypsums and halites are observed. In the Rączyna area, the Skole Unit, which belongs to the Flysch Carpathians, includes Lower Cretaceous Spas Beds, Upper Cretaceous Variegated Beds and Siliceous Marls, Upper Cretaceous/Paleogene Inoceramus Beds, Eocene Variegated Shales and Hieroglyphic Beds, Oligocene Menilite Beds and Oligocene/Lower Miocene Upper Krosno Beds, all showing variable thicknesses [27]. The Quaternary sediments cover the autochthonous Miocene deposits but, locally, these strata rest upon the strata involved in the Carpathian-Stebnik Overthrust. These are sandy weathering clays, gravels and sands together with paludal and glaciolacustrine deposits of total thickness in the range from 0 to 60 m.

3. Methods and Materials

3.1. Velocity Model Validation Method

For many years studies on various aspects of multivariant integrated inversion of independent datasets, based upon the same or physically related petrophysical parameters, have been carried out [1,2,3]. However, an important issue lies in integrating data that represent diverse properties of rock formations, even though they lack direct physical connections, making it impossible to establish precise analytical links for such datasets [4,5]. The solution to this problem was proposed by Meju [3], who developed a methodology for effective linking of non-correlated physical properties of the rock formations in an analytical manner, based upon the so-called structural resemblance of parameters distribution [30]. Their concept was built upon the Electric Resistivity Tomography (ERT) methodology combined with the Refraction Seismic Tomography (SRT) [3,31]. The analysis included 2D models of seismic wave velocity and resistivity. The idea of a 2D inversion solution proposed by Meju et al. [3] took advantage of the fact that even in the absence of analytical links between various physical properties of the rock formation; the possibility was expected of their structural resemblance within the fragments or the entire of the analyzed geological formation. This resemblance resulted from the statistical correlation of the changes in their distribution trends in the zones relatively uniform structurally and genetically [30,31]. The usefulness of this correlation methodology of geophysical models was confirmed by Jegen et al. [21] and Moorkamp et al. [9,30], who developed the code for 3D inversion of magnetotelluric, gravity and seismic data [9,30]. The physical properties of geological formations vary with their spatial position and changes may occur in any direction due to the phenomenon of anisotropy. Hence, at each point of the analyzed space, changes may be characterized by their intensity and direction [3]. The similarity of these changes determines the structural convergence of the distribution of the two parameters in the rock formation. These parameters are mathematically expressed as a gradient field vector of the so-called Cross Gradient (CG) [1,2,3,9,30].

In order to test the SJI as a tool for validation of the interval velocity model applied to run the PreSDM, the results of the long-offset refraction survey and seismic survey were used, provided within the INGA 2 project. The workflow for this procedure is presented in Figure 3. The implementation of the test was enabled by the application of a schema generally known and described by other authors, allowing the utilization of archival seismic and gravity data. It should be emphasized that, until now, these data have not been used in the seismic data processing within this part of the Carpathians. The available schemas were supplemented with a procedure for creating independent velocity models, which serve to address ambiguities resulting from the exclusive use of velocity fields that are intermediate results of the seismic data processing itself. As is known, the acquisition process and initial processing procedures are subject to certain statistical errors and necessitate simplifications in algorithms for the applied procedures.

Refraction datasets were determined for 0–6 km offsets using the reflection seismic records. Initially, at the stage of construction of shallow velocity models the short-offset part of the data was used for calculation of static corrections. In the SJI procedure described above, these datasets were responsible for the stabilization of the model. Moreover, they also provided the proper conditions for calculations of the Cross Gradient (CG) [1,2,3,9,30]. The application of that function caused that for the part of the model, where the sensitivity of the gravity method was the highest [30,32], the velocity was simultaneously determined by two mutually agreed inversions within the SJI. This contributed to the detail of the VID model, where the strongest impact of incorrectly defined interval velocities on the PreSDM results was found. The obtained data allowed modifying the model in depth interval from several hundreds of meters down to 3 km below the surface, as confirmed by the result of ray tracing. However, the final penetration depth of seismic refraction depended on both the nature of the data and the geological conditions. The second type of dataset used for validation of the VID included residual anomalies of the gravity field. Taking into account that calculations were carried on in a 2D domain, the values of residual anomalies were extracted from the Bouguer anomaly map along the seismic traverses [32].

3.2. Joint Inversion as a Tool for Velocity Models Validation

The SJI procedure applied to the validation of the VID was carried on according to the scheme shown in Figure 4, which demonstrates the general operational principle of the discussed algorithm. In this algorithm, two independent geophysical datasets, i.e., the residual values of the gravity field and the times of first breaks for 0–6 km offsets were synchronized to simultaneous inversion. Precisely, the mentioned above datasets were processed simultaneously with a single, coherent algorithm. The application of the SJI procedure was mainly aimed at generating a credible distribution model of petrophysical parameters of rock formations based on available data not analyzed previously in a simultaneous mode. In this SJI solution, the link function of two independent datasets was the spatial distribution of the Cross Gradient (CG) for two different model vectors (Equation (1)) [1,2,3,9].

$$\mathrm{C}\mathrm{G}\left(m_1,m_2\right)=\nabla m_1\times \nabla m_2$$

(1)

A simplified explanation of the simultaneous inversion methodology is demonstrated in the form of Equation (2) of a Join Objective Function (JOF). This formula represents how the minimization in the optimization algorithm is being released below, where the νΦcoupling (m) parameter is responsible for the integration element [3,30,31]. This parameter defines mathematically the relationship between the parameters of the model which can be related to CG:

Φ_joint (m) = Φ_d (m) + λΦreg (m) + νΦcoupling (m)

(2)

where:

Φ_joint (m)—join objective function

Φ_data (m)—fitting coefficient of calculated and measured data

Φ_reg (m)—regularization condition

(m)—model vector

λ—lagrange multiplier used in regularization process

νΦ coupling (m)—parameter, which defines mathematically the relationship between parameters of the models used in calculations.

Application of the SJI with the CG resulted in the generation of two models, which determined the spatial variability of velocity and density. An important problem in the application of the SJI procedure with the previously described datasets for greater depths of rock formations was the incommensurability of sensitivities of gravimetric and seismic refraction methodologies. The depth range of gravity data is much greater but after exceeding the depth range of sensitivity for the refraction method, the joint inversion becomes a cooperative inversion.

In each calculation loop, the average fitting error of calculated and measured values for both parameters was estimated together with the average value of the cross-gradient for both models. Thus, the data quality of obtained results could be evaluated and the range of input model refining could be determined. The aim of the inversion procedure was to minimize the fitting error down to the lowest possible value. In the case of geophysical data, the goal was to reconstruct the measured physical field from a series of simple models. The minimization of cross-gradient values meant obtaining structurally convergent geophysical models [3,9,30]. Therefore, the uncertainty was reduced, caused by the reconstruction of rock formations using the developed geophysical models. The introduction of data recorded independently on the leading methodology enabled the researchers to verify one method with the results of the other [3,9,10,15]. According to the assumed sensitivity reduction of gravity method with depth, the applied weighting algorithm allowed the researchers to filter out the effects resulting from the projection of anomalies distribution at shallower depths. Hence, the effect originating from residual, fast-changing values of gravity anomalies was distributed near the surface whereas the slow-changing values and regional trends were distributed at depth. During the simultaneous procedure of data integration, the SJI was run in order to construct the Multiparameter Geophysical Model (MGM), which was a distribution function of density and interval velocity parameters measured in the wells. As a result, separate distributions of parameters were obtained from the MGM although this was a rather serious simplification of the joint inversion procedure, which required additional parameterization. In practice, it was possible to apply additional weights, which defined the impact strength of one method on another [30,32]. In fact, a model based on one algorithm affected the inversion solution in the other algorithm. Weighting procedure could be used also at the level of datasets where their fitting was analyzed, i.e., a hierarchy was established in which a parameter is to be fitted with the lowest error [30,32]. Such a weighting procedure is used, among others, in joint inversion of data acquired from refraction seismic and gravimetry. Generally, the seismic method was weighted as the leading one because gravity is based upon the distribution of the potential field, which cannot directly determine the vertical variability of density in geological formations [33]. The simultaneous inversion allowed automatic modification of the density field based on seismic tomography in the shallow zone, which was crucial to PreSDM whereas the gravity inversion enabled the researcher to verify the model and to supplement geophysical information in places where seismic survey could not be run in deeper parts of geological formations. When the weights are determined, parameters such as resolution, sensitivity and prospection depth of each method must be considered [30,32,34]. The well-known weighting function W(z) is used for this purpose described below:

$$W\left(z\right)=\frac{1}{{(z+z_0)}^{\frac{\beta }{2}}}$$

(3)

where: z—depth, z₀ and β—parameters depended on model parameters [32].

3.3. Gravity Data Characterizations

The results provided independent data carrying information on the structural variability related to changes in the density contrasts of rock formations [6,9,17,35,36,37]. The characterization of the physical phenomena, which are essential for this geophysical methodology results in its high sensitivity in the near-surface zone but also in a decrease in sensitivity with increasing depth [32,34]. This is a classic phenomenon of geophysical data analysis based on the potential fields. It is associated with a decrease in the resolution of this methods with depth. However, the published results of analyses demonstrated that this limitation did not adversely affect the validation of VID models [1,2,9,30]. The analysis of previous publications and unpublished reports enabled the researchers to conclude that the specificity of the gravity field distribution in both the western and the eastern regions of the Polish Carpathians was diverse [6,32,34,35,36]. The most complicated was the Bouguer anomaly pattern for the eastern region. The gravity survey carried out in the Carpathian area revealed a number of local anomalies of irregular shapes [6,33,35,36] (Figure 5A). These anomalies indicated the presence of geological structures composed of rocks with different densities [6,33,35,36,37]. The analysis of the gravity field pattern, combined with the varying density of particular lithological units [16,35,36], showed that relatively negative anomalies were caused by various formations.

Disturbances of Bouguer anomaly distribution were caused by a strong regional gravity field generated by a deep-seated basement. As a result, anomalies of lower amplitudes and local ranges, which represent structural elements or lithological units of lower order, became less visible. In order to separate and highlight the gravity effects generated by these structures (or units), regional gravity anomalies were calculated (Figure 5A).

In order to illustrate the influence of the basement, distributions of the regional component of the recorded gravity field were determined. The transformation parameters were selected in such a way that extraction was possible of disturbing bodies in the particular depth intervals. Therefore, the filters were tested suitable for the generation of regional maps. The results are presented in Figure 5. On their basis, the residual values were determined and used to carry out the inversion tests. At this stage of the project, the aim was to filter out the effects of the deep basement and the isostatic movements related to the lithosphere/asthenosphere boundary [35]. Taking into account the fact that the analysis included also the basement structures, it became necessary to ensure that the gravity effect generated by differences in the density of basement formations is included in the data input into the SJI. The results of the tests enabled the researchers to select Option No. 1 of residual Bouguer anomalies (Figure 5B) for inversion calculations. Data analyses published earlier demonstrated that the values of residual anomalies calculated for various depth intervals increased with the increasing depths of measurements [6,33,35,36,37], as illustrated in Figure 5. Variability analysis of residual anomalies with the depth of measurements allowed us to conclude that the source of the Carpathian gravity depression [36] (Figure 5A,B anomaly I) is the mass deficiency occurring within a significant depth interval, from flysch formations to deep crystalline basement. This deficiency determines both the shape and the amplitude of regional gravity anomaly. The anomalous decrease in the value of the gravity field in the Carpathian foreland reflects the diverse depth of the sub-Tertiary basement covered by less dense Miocene molasses although Bojdys et al. [35] proposed alternatively the isostatic effects as a factor modifying the gravity field in that area. Despite the interpretations, the influence of the basement effect on the modeling of the gravity field had to be minimized. Although the study area is limited (Figure 5), additional analyses were carried on in order to include the gravity effects generated even by far distant zones (Figure 5A,B anomaly I) and, thus, to ensure proper boundary conditions for gravity modeling. This resulted from the fact that the gravity method is very sensitive to superposition effects [6,35]. For this reason, the anomaly described in the literature as ”the Carpathian gravity depression” [35] (Figure 5A,B anomaly I) strongly influenced the values of the Bouguer anomaly in the entire analyzed area. Particularly, it influenced the results of modeling in the marginal parts of the Carpathians where the boundary line between the Carpathian structures and their foreland [38] was most clearly marked by the presented distribution of residual values of Bouguer anomaly along the analyzed seismic traverse (Figure 5A,B anomaly II). To the NE of this line, the strikes of structures built of Miocene formations in the Carpathian Foredeep radically changed their directions in relation to the structural directions typical of the Polish Eastern Carpathians. The block structures developed within the sub-Miocene basement were separated by zones of discontinuous deformations [33,35,36]. Residual anomalies were characterized by symmetric shapes of transversal cross-sections and high-intensity changes within the main Carpathian gravity depression. However, these changes clearly decreased to the NE from the marginal zone of the Carpathian Overthrust [29,38]. The central part of the profile was dominated by anomalous zones corresponding to structures of a smaller depth range whereas in the marginal zones, the structures of a greater depth range were more common. However, regardless of the type of transformation, the directions of their strikes and their positions on the profile were highly consistent.

3.4. Velocity Model Building

The principal generation method of the optimal VID model is the application of MVA. Based upon the residual hodograph curvatures occurring on trace gathers after the PreSDM, the Residual MoveOuts (RMO) are calculated, from which the VID initial model is updated using the iteration method [39,40]. Hence, the resultant VID model is strongly dependent on both the seismic data quality and the processing methodology. In order to solve this problem, within the INGA 2 research projects, the generation of independent velocity fields was applied. Based upon many tests, it was assumed that the initial velocity field for PreSDM can be the Parametric Velocity Model (PVM) presented in Figure 6. The input datasets for model construction were the interval velocities calculated from average velocities recorded in selected wells from the Eastern Carpathians and the Carpathian Foredeep. Considering the regional nature of the velocity field and the varying quality of average velocity records, the static parametric modeling was carried on for layers up to 1000 m thick (Figure 6A). This enabled the interpreters to apply block averaging of velocity values and to eliminate the errors resulting from the varying quality of velocity measurements in wells.

Modeling applied the anisotropic kriging methodology with the exponential variogram [41] and the preferred interpolation azimuth 120 degree, i.e., consistent with the regional strike of geological structures in the study area [Analysis based upon geological maps provided by the PIG-PIB [42] and on the basic geological map after Jankowski [38]. Moreover, a regional trend of velocity variation with depth was determined, which was then used as a secondary variable. Such a procedure was necessary due to the relatively small number of wells together with the steeply dipping bedding, which would cause an underestimation of velocities in the deep zone and their overestimation in the shallow zone while ignoring the depth trend of the velocity. For model construction in the depth domain (Figure 6B), the interval velocities were used, obtained from the measurements of average velocities together with shallow velocity models derived from the Seismic Refraction Tomography (SRT) survey. Therefore, datasets provided by the shallow velocity model could be included in the petrophysical modeling [3,17,18,31]. Such an attempt enabled a better representation of the velocity variations in the shallow zone.

3.5. Methodology and Parametrization of Kirchoff PreSDM

Prestack depth migration was performed using the Kirchhoff summation algorithm which performs a migration by applying a Green’s function to each CDP location using a traveltime map. Kirchhoff Migration is commonly implemented following the technique outlined by Schneider [43] for the integral solution of the scalar wave equation. The migrated image is constructed by summing weighted amplitudes along diffraction curves. These curves are determined by traveltimes from surface to subsurface scatters. Kirchhoff depth migration estimates traveltimes from a velocity model (supplied as vertically and laterally variant interval velocity field in-depth domain), using ray-tracing or solving an expression such as the Eikonal equation.

The Kirchhoff PreSDM procedure was performed in two steps. In the first step, the travel times were calculated from an imaged underground point to shot and receiver points using the given model of interval velocities. In the second stage, the weighted summation of trace amplitudes was computed along the supposed trajectory of the diffraction surface [39,40]. The calculation method of arrival times plays a key role in the results of migration because it determines the correctness of the shape of the diffraction surface. The calculation algorithm for traveltimes was based on the solution of the eikonal equation, which is a high-frequency approximation of the wave equation [44]. It comprises the two methods: the seismic ray tracing [44,45] and the previously developed finite difference [46].

For the purposes of our task, the input data consisted of fully preprocessed CMP gathers, located on floating datum, after offset regularization and sorted into the specified group bin size. Migration parameters were determined based on tests, especially migration aperture, maximum dips, maximum frequency, traveltime computation method, offset groups and output depth step.

Finally, the following main parameters were used: offset group: 120 m, aperture size: 7000 m, maximum dip up to 75 deg, maximum offset: 10,080 m, ray tracing method: Eikonal solver, frequency range: 8–90 Hz, output sampling rate: 5 m.

Additional processing was performed on the migrated gathers. First, the angle mute up to 45 deg was conducted, and then residual moveout correction RMO based on interpreted horizons was computed and applied. Maximum moveout for nominal offset was put as +/−800 ms. Trim statics calculation was the last step performed on CRP gathers before stack. Poststack processing consisted of signal-to-noise enhancement realized using coherency dip filters in the F-X domain.

The imperfection of Kirchhoff PreSDM methodology in the areas of complicated tectonic structure results from the fact that under such conditions, the values of Green’s function cannot be unambiguously determined due to the presence of multiplied wavefronts generated by a large number of discordant reflectors [40,44]. Additionally, in the geological structures strongly refracting the acoustic waves, the true reflectors are accompanied by false reflection horizons resulting from the migration procedure itself and unrelated to the errors contained in a given seismic velocity model. If such false horizons are accepted as true ones they can be erroneously included and analyzed in the estimation procedures of seismic velocities [40,47,48]. If seismic datasets originated from the Carpathians and the Carpathian Foreland are considered the main obstacles in their usage are heterogeneities in data recording procedures. It is revealed by a large number of successive omitted shotpoints, which additionally hinder the correct calculation of velocity-depth models then applied to the computation of static corrections. Their erroneous values may cause unreal changes in the course of hodographs, which, in turn, may deteriorate the imaging of analyzed rock formations with the PreSDM procedure using the floating datum migration variant. This datum is an equivalent of the terrain level, which represents the varied topography. In the analyzed zone of the Carpathian Overthrust, significant variability of seismic velocities is observed in the Low Velocities Layer (LVL), represented by high values. Thus, the correct attempt at the problem of changes in the near-surface zone as well as the correct determination of first-order static corrections [48] is crucial to the successful execution of both the Migration Velocity Analysis (MVA) and the mentioned above PreSDM procedure [39,40].

4. Research Results and Discussion

According to the methodology presented by Colombo [1,2], if the seismic sections after PreSDM meet predetermined quality standards, the extracted velocity model is regarded as final. Otherwise, the model should be modified with the iterative process of MVA or the procedure should be returned to the SJI stage but with changed parameters. However, in the described approach, this methodology was applied to correct and verify the velocity model developed with the parametric modeling. The SJI procedure performed on both the refraction data and the residual anomalies of the gravity field enabled the researchers to develop the VID model. Therefore, the VID model (PVM), which is based on parametric modeling, was validated. Unfortunately, the parametric model has some drawbacks resulting from the limited amount and irregular distribution of well logs as well as from the extrapolation of velocity values measured in wells to areas where no wells exist, thus, no local values are available. The results of quality analyses of obtained velocity fields are shown in Figure 7.

Figure 7A presents the raw VID model constructed based on parametric modeling. After converting this model to density using Gardner Equation [49,50], the gravity effect was calculated and shown as a red line in Figure 8. It can be observed from this graph that the model provides a similar gravitational response compared to the actual data. Relatively good agreement is visible on the sides of the area of interest (Carpathian region and marginal zone of the pre-Carpathian basin). However, discrepancies are evident within the Area of Interest (AoI). These observations necessitate the need to correct the P-wave velocity model to reduce its ambiguity. Despite the observed deviations, it still served as a good initial model for further analysis. Three zones marked with numbers “1” to “3” are indicated on the velocity cross-sections. These zones represent the main areas of changes introduced into the initial P-wave velocity model. Based on these zones, it can be observed that using tomography alone leads to a smoothing effect on the velocity distribution, which is not the desired outcome. The application of SJI with the same parameters as the seismic tomography provided a more stabilized solution and allowed for velocity corrections, especially in Zone “1”.

For additional verification, a Single Gravity Inversion (SGI) [30,32] was performed, which means that only gravity data were used for inversion. As seen in Figure 8, a very good fit was achieved between the calculated and observed data. However, this approach has its limitations, as the inversion algorithm adjusts the model to converge the misfit error at the expense of obtaining a realistic final model [9,30]. Therefore, this result is presented as an example of quality control. Consequently, the misfit error in the gravity data during SJI is larger (approximately 1 mGal). This is closely related to the structural constraints imposed by the CG during the joint inversion of refraction and gravity data [3,9,30]. This algorithm ensured the stabilization of the inversion solution and proper regularization for the employed models, as well as the resolution and sensitivity of the geophysical data used [30,32,34].

To verify the solution and check whether the combined inversion process preserved the velocity distribution as defined in the initial model PVM, a simple comparison of interval velocities with the values extracted from models obtained through SRT and SJI was carried out. Additionally, the extracted curve from the initial model PVM was plotted on the panels. These comparisons are presented in Figure 9, showing the velocity distribution for two wells in the Rączyna deposit area [27]. The plots demonstrate relatively good agreement between the curves while maintaining the trends of vertical velocity variation. Particularly good agreement is observed at greater depths.

Based on the models presented in Figure 7, PreSDM migrations were performed without MVA corrections, as the purpose of the described test was to verify the utility of the PVM models as initial models for the first iteration of MVA. Additionally, the suitability of these models as the main solutions for VID fields intended for PreSDM was verified. The conducted test demonstrated that the developed models (Figure 7) perform well in depth migration. Figure 10 presents the results of PreSDM for a section of seismic profile M-5 (Figure 1). An analysis of the quality of the obtained results was conducted based on the comparison of the seismic sections after PreSDM migration. At this stage of the work, the focus was on evaluating the seismic structural image, which is an important aspect of data processing that leads to the proper identification of the geological medium. Several studies carried out in this area allowed for the verification of the effectiveness of the applied velocity field validation methodology [7,24,25]. It should be emphasized that the data were recorded along a regional profile aimed at deep structural reconnaissance. The quality of the obtained seismic image is comparable to the results of processing along seismic lines where the acquisition parameters were designed to recognize the Jodłówka and Rączyna reservoir zones [27]. It is worth noting that this result was obtained in the first iteration of PreSDM.

The obtained results indicate that the procedure of creating independent models for PreSDM, supported by the combined inversion of independent geophysical data, is a good tool for assisting seismic data processing. In Figure 10, Panel A shows the seismic sections created in the CDP domain from the seismic gathers after PreSDM, calculated using the basic velocity model PVM (Figure 7A). It exhibits the most consistent phase (signal) representation of reflections, mainly in the flysch part of the thrust zone, compared to the other sums (Figure 10B,C). However, the analysis of the seismic data quality for exploration purposes, specifically in the Jodłówka and Rączyna [27] reservoir zones (Figure 2), indicates that the structural image in this section is the weakest. This is particularly evident in the northern part of the section, where versions B and C (Figure 10) show continuity of horizons, while version A (Figure 10) shows them being scattered and difficult to interpret. These structures are associated with the Jodłówka and Rączyna reservoirs. The seismic sections presented in Panel B (Figure 10) provide a correct and comparable structural image compared to the production processing versions (Figure 2). Panel C, in comparison to Panel B (Figure 10), presents a more consistent image of deeper parts of the Carpathian Flysch (south of the thrust), with fewer attenuated areas in the seismic record. In this result, the amplitudes observed along the Miocene reflections are more stable, whereas in Panel B, significant amplitude variations are observed. Additionally, north of the R-8 well, version C exhibits a seismic image characterized by the presence of continuous reflection sequences, which are disturbed in version B (Figure 10 and Figure 11). In terms of accurately representing the structural variability of Miocene deposits associated with reservoir presence, reproducing the amplitude variations and geometry of reflections is crucial.

In the details, attention should be mainly focused on the quality and continuity of the reflectors. For this purpose, enlarged sections of seismic profiles (Figure 10) are presented in Figure 10. Arrow symbols indicate zones where seismic imaging improvement is observed. Zone “I” (Figure 10 and Figure 11) indicates a location in the northern part of the area where significant improvement in the continuity of seismic reflectors. This is likely due to the proper stabilization of the velocity field validation process in this part of the profile. These changes are marked as anomaly “I” in Figure 10 and Figure 11. Other crucial areas (“II” and “IIa”) present changes within the R-8 well area in the “Rączyna-Jodłówka” reservoir zone [27], where the distinct lining of seismic reflections defining this reservoir is observed (Figure 2). A similar situation occurs for the zone “III” in the “Rączyna S” reservoir area [27]. It corresponds to the previously described zone where an improvement in the continuity of the seismic record is observed. Considering the methodological nature of the conducted work, such a state of affairs allows for a precise interpretation of reservoir zones. The last observed zone “IV” is associated with the zone of reflections at the contact between flysch deposits and the crystalline basement. Using the validated VID model with SJI, a structural seismic image was reconstructed.

It is worth emphasizing that a significant improvement in the seismic image has been achieved, as clearly visible in comparison to the archival results presented in Figure 2. The analysis of the above figures demonstrates the crucial role of the VID velocity field in depth processing. It can also be seen that by using the described methodology, the obtained result approaches the quality of 3D seismic processing, which provides significantly better results than 2D seismic. The results presented in the above analysis were compared to previous processing results along the Transcarpathian seismic profile (POLCRUST), confirming the correctness of the adopted data processing sequence [7,24,25,26].

In order to obtain a smoother seismic image, an additional post-processing step was performed. Figure 12 shows a comparison of examples of CDP gathered from the Carpathian thrust area. Significant improvement in signal coherence can be observed on these gathers, as well as the effects of procedures such as Angle Mute, RMO update, Trim Statics, and FX filtration [40].

Figure 13 presents the result of PreSDM migration based on the seismic components shown in Figure 12, which are the result of additional processing incorporating procedures to improve signal coherence and residual RMO corrections, enhanced by the application of trim statics. This approach allows for obtaining a seismic image that is suitable for geological and reservoir interpretation at an early stage of depth imaging. Key locations on the seismic image (Figure 13) are marked with circular ellipses, confirming the correctness of the performed analyses and validating the use of independent procedures for creating VID models. Regions of significant improvement in imaging shallow and steep outcrops of flysch formations are indicated by orange ellipses. Imaging such areas is often challenging due to acquisition gaps, which can cause difficulties in velocity determination in a classical processing approach. The red ellipse represents a shallow zone beneath the thrust, where reflections indicating its boundary are accurately depicted. The blue and yellow ellipses mark a similar zone with similar structural characteristics related to middle Miocene levels. In standard processing results, this zone is characterized by muted seismic recordings. Similarly, in the zone marked by the green ellipse, an improvement in the quality of seismic recording is observed. The grey ellipse illustrates the improvement in the quality of reflections from autochthonous Miocene sediment layers wedging on the crystalline basement. The results of advanced processing of the 2D data presented above are comparable to the results of 3D seismic processing (Figure 2B).

5. Summary and Conclusions

Due to numerous limitations associated with the process of building the VID model, it should not be expected that the velocities determined using seismic methods are free from errors. Therefore, they should be approached with limited confidence. Consequently, it is difficult to claim that the obtained values are true physical parameters, especially when there is no possibility of verification using high-quality well data. However, by using modern and advanced inversion methods (such as SJI) supported by MVA and considering velocity changes in the LVL, it is possible to construct a correct and realistic distribution of rock velocities. It should be noted that both validation procedures used for the VID models, namely SJI and MVA, like any inversion procedure, are inherently ambiguous and heavily depend on the initial model. The obtained results prove that modeling procedures (i.e., SJI) may cause fewer errors in comparison to standard algorithms. This is particularly relevant in the case of complex geological media and limited well information. In such cases, parametric modeling takes into account the three-dimensional variability of the geological medium during the construction of the initial velocity model. As a result, more accurate estimates of horizontal and vertical velocity changes in the medium, are obtained. Furthermore, the presented methodology eliminates the problem of misalignment between 2D profiles resulting from separate, traditionally applied analyses. In such a situation, similar to the time-to-depth conversion process, using the PreSDM procedure allowed us to obtain a solution that was validated by independent geophysical data. Nevertheless, the final decision regarding its correctness always depends on the person performing the geological interpretation. The presented methodology does not negate typical seismic reflection processing steps as they are still methodologically correct. However, this analysis should be regarded as a procedure for preparing a proper starting velocity field for the first iteration of MVA. Validated VID models using the SJI inversion procedure undoubtedly constitute a valid component introduced into the scheme of depth seismic imaging.

In conclusion, it should be emphasized that despite the presented advantages and disadvantages of the described methodology, an example of a PreSDM result presented in Figure 13, demonstrates that a single iteration of the VID velocity field validation procedure provides accurate outcomes. As referenced in the article, the comparison with archival seismic data shows an improvement in reflection imaging. Furthermore, enhancement in seismic reflection continuity can be expected through additional iterations using the MVA procedure. It should be emphasized that only a simple processing of seismic cross-sections was performed. As shown in this seismic section, it is feasible to perform an improved structural interpretation of the Carpathian thrust zone. Even in this theoretically preliminary version of PreSDM, the reservoir zones of Rączyna are fully interpretable. Particularly significant for structural interpretation are the well-reconstructed reflectors below the thrust boundary. The fact that achieving this result is not possible with standard procedures, the initial MVA iteration presents the validity of the presented methodology and its favorable results. Furthermore, it can be inferred that if additional geophysical data are available, they should be integrated into the seismic data analysis procedure, even as reference data.