Selected aspects of modern seismic imaging and near-surface velocity model building in the area of Carpathian fold and thrust belt

Despite the increasing technological level of the reflection seismic method, the imaging of fold and thrust belts remains a demanding task, and usually leaves some questions regarding the dips, the shape of the subthrust structures or the most correct approach to velocity model building. There is no straightforward method that can provide structural representation of the near-surface geological boundaries and their velocities. The interpretation of refracted waves frequently remains the only available technique that may be used for this purpose, although one must be aware of its limitations which appear in the complex geological settings. In the presented study, the analysis of velocity values obtained in the shallow part of Carpathian orogenic wedge by means of various geophysical methods was carried out. It revealed the lack of consistency between the results of 3D refraction tomography and both the sonic log and uphole velocities. For that reason, instead of the industry-standard utilization of tomography, a novel, geologically-consistent method of velocity model building is proposed. In the near-surface part, the uphole velocities are assigned to the formations, documented by the surface geologic map. Interpreted time-domain horizons, supplemented by main thrusts, are used to make the velocity field fully-compatible with the litho-stratigraphic units of the Carpathians. The author demonstrates a retrospective overview of seismic data imaging in the area of the Polish Carpathian orogenic wedge and discusses the most recent global innovations in seismic methodology which are the key to successful hydrocarbon exploration in fold and thrust regions.


INTRODUCTION
The main focus of this study are mutually-dependent aspects of seismic imaging and velocity model building in fold and thrust belt areas. The presented work mainly refers to the near-surface part of the Carpathian orogenic wedge, where the velocity model is difficult to predict. In seismic data processing scheme, an industry-standard approach is to build it using first-break refraction tomography. Nevertheless, in complex geological settings, this approach is debatable (Habeebuddin & Pandey 2006). Other related methods, such as: seismic modelling or illumination studies, require a less-detailed but geology-consistent velocity field. In these cases, as presented below, a more desirable solution is the utilization of interpreted thrusts, surface geologic map and external sources of information (e.g. the uphole measurements).
Seismic exploration projects performed in fold-and-thrust belts (FTBs) encounter numerous difficulties. A typical and frequent problem is the failure of conventional data processing and poor imaging of the thrust structures caused by the limitation of time-domain procedures to account for lateral velocity variations. Under simplified assumptions, the generation of false structures and mispositioning of steeply dipping seismic events is unavoidable. These issues become particularly important when the medium underlying the thrust is the main exploration target (Rigatti et al. 2001). Complex wave paths through the thrust-sheets, together with their high P-wave velocities, result in artificial pull-ups of the sub-thrust horizons and distortions of the seismic image.
Another problematic aspect is the presence of noise in the recorded data. Source-generated noise, such as: direct, refracted, guided and surface waves, may significantly interfere with the shallow reflections and dominate the shallowest portion of the seismic data (Schmelzbach 2007). Due to similar dominant frequencies to the primary reflections, the guided and direct SV-wave trains may be difficult to suppress during the data processing stage and may consequently be misinterpreted as real, reflected events.
The undesirable effect of the surface waves can be analysed by studying the results of the conventional time processing performed on short-spread seismic data (Yilmaz et al. 2010). Such images are overwhelmed by chaotic, random response and an evident lack of a coherent signal. Limited-offset acquisition, conducted under complex subsurface conditions, is the main cause of the poor imaging. Seismic processing becomes even more complicated when one must deal with anthropogenic noise (Curia et al. 2017, Liu et al. 2018. Such noise, produced by wells, electrical transmission lines, pipelines, and traffic, is difficult to attenuate due to the similar frequency range to P-wave reflections. A low signal to noise ratio is a common feature of seismic data sets, acquired in orogenic wedge zones (Gittins et al. 2004, Curia et al. 2017. Due to all those reasons, the data processor must discriminate between the noise and useful data throughout the whole project, which is both time-consuming and subjective. Highly variable elevation is the next element that hinders seismic imaging in FTBs (Habeebuddin & Pandey 2006, Liu et al. 2018. The modelling studies indicate that in the presence of rugged topography, even if the overburden is considered homogeneous, 3D orthogonal acquisition geometries may correctly image only the deep targets and fail with shallow ones (Gerea et al. 2011). Varying topography may cause the scattering of the bodywave refracted energy and generate reverberations. Such noise is spatially aliased in all domains and difficult to suppress, particularly in the case of lowfold sparse 3D surveys. It carries practically no information about the near-surface medium and has amplitude higher than desired signal.
Laterally-changeable conditions complicate the corrections for near-surface effects and the determination of the true velocity model (Gittins et al. 2004, Habeebuddin & Pandey 2006). An analysis of the first arrival refractions, routinely used in many seismic processing projects to estimate the geometry of the low-velocity layers (LVL), is no longer considered appropriate. This technique fails due to the assumption of vertical ray paths, which is not the case in the discussed regions. Ray path bending occurs in the presence of folds and steeply dipping thrust planes.
Seismic data from FTBs may also suffer from poor illumination (Soleimani et al. 2012). The phenomenon may be caused by dipping, near-surface ductile rocks, such as: salt, shale, or gypsum. The total energy transmitted through these rock types is low, therefore the reconstruction of the true image below them, is a challenging task.
There is no single remedy for the mentioned problems. Various data processing solutions, verified positively in recent projects, can be categorised as follows: pre-stack depth migration (PreSDM), stacking techniques, advanced statics solutions, forward modelling, inversion, the use of gravity data, and diffraction imaging.
Numerous publications document the advantage of PreSDM over time-domain migration (Schmid et al. 1996, Rigatti et al. 2001, Issac & Lawton 2008. Isotropic PreSDM provides better fault definition and leads to structurally more accurate images of the thrust-sheets and subthrust intervals. The anisotropic form of the depth migration has the potential to further improve the lateral positioning of the geological target (Vestrum & Lawton 1999). Nevertheless, this process is more time-consuming and requires the iterative selection of Thomsen's parameters: epsilon and delta. The choice of the above values is crucial, but in many practical situations it could only be made by the subjective assessment of the gathers and migrated seismic sections, with limited well control.
Instead of standard Kirchhoff-type migration, Common Reflection Angle Migration (CRAM) may be utilized to enhance the imaging of strongly deformed orogenic wedges (Shiraishi et al. 2019). This method generates common image gathers directly in the local angle domain and has the advantage of working better in the areas of complex wave propagation and poor illumination (Koren & Ravve 2011).
The PreSDM method involves a great amount of work in the process of velocity model building. The initial model should be consulted with an experienced interpreter and incorporate as many geological constraints as possible (Vestrum 2007). The role of the structural geologist is not limited to performing well-to-seismic correlation, picking time horizons, and preparing first-guess velocity volume. The interpreter should also participate in the process of the selection of migration parameters, which involves trade-off decisions, whether to honour the flatness of the gathers or the improved stacked image (Vestrum 2007, Issac & Lawton 2008. Creating the uppermost part of the PreSDM velocity model is equally important and demanding procedure because seismic data typically contain very little information about near-surface velocity. As mentioned earlier, first-break tomography is considered improper for complex geological settings, therefore other techniques, such as: time-variant static corrections based on wave-equation datuming, can be used instead (Habeebuddin & Pandey 2006).
As an alternative to common midpoint (CMP) method, that has been the key element of seismic data processing for the past several decades, other stacking processes, such as: Multifocusing (MF) (Curia et al. 2017) or Common Reflection Surface (CRS) (Mann et al. 2007, Soleimani et al. 2012 can be performed. In highly deformed FTPs, the statistical effect of CMP summation is insufficient because too few traces are stacked, hence the signal to noise ratio remains low. The MF process is applied to conventional time-domain gathers, and as a multiparameter procedure, operates on a portion of the Fresnel zone around each original CMP (Berkovitch 1999, Berkovich et al. 2008. This process can be considered a partial stack along the calculated surface that generates the enhanced gathers, for their further use in seismic migration. CRS stack technology can similarly be regarded as an extension of the CMP stack and velocity analysis (Mann et al. 2007). Whereas conventional velocity analysis uses a stacking trajectory described by a single parameter (stacking velocity) within a gather, CRS method utilizes adjacent CMP gathers and is defined by three parameters. The incorporation of a greater number of traces than in traditional stacking leads to the enhancement of the useful waveforms.
The development of processing procedures may not be sufficient to improve the imaging of the complex targets, if they don't go hand in hand with the appropriately-planned acquisition. In such geological settings long-offset acquisition schemes need to be tested. Their optimal parameters may be defined by a feasibility study that includes a combination of forward modelling and PreSDM (Colombo 1999). Such study starts with ray-tracing through a predefined velocity model that is aimed at the computation of actual fold at the target reflections, with special consideration to long-offset arrivals. It allows focusing effects and determination of low-coverage areas to be analysed. The next step is the generation of synthetic data via finite difference FD modelling. The migrated synthetic images support the choice of the appropriate offset range and are used to investigate the effects of the introduction of long offsets into the poor coverage regions.
Conventional spread length seismic data acquisition schemes often fail to image complex imbricate structures associated with thrust tectonics (Yilmaz et al. 2010). If the ground-roll energy dominates the recorded wavefield at the near offsets, the exploitation of wide-angle reflections which are present outside the surface-wave cone may be the only solution for the signal to noise ratio enhancement. In some situations, P-wave reflections from the deeply-buried horizons hit the overlying thrust with critical or over critical angles. In such cases, the incident rays are reflected down again and no energy from these layers is recorded. Additionally, the waves reflected from the thrust planes may migrate much further than the position of the last active receiver (SeyedAli & Zabihi 2012). For that reason, forward modelling remains the only method to estimate the maximum offset and correct value of the migration aperture.
Other techniques, such as the PSTM modelling, can be utilized to gain an understanding of the most common structures encountered in FTBs: fault-bend faults, fault-propagation folds, and trishear fault-propagation folds (Li & Mitra 2020). The first fold type is generally well-imaged on seismic sections, but in the presence of the latter, strong pull up effects of the underlying layers and poor reflections from the steep front limb of the fold can be expected. The pre-stack depth migration applied to synthetic representations of these structures produces better and more geologically-consistent images than pre-stack time migration PreSTM, assuming that the accurate velocity model is known.
There are many other recent inventions developed to improve seismic imaging in the discussed regions: waveform and traveltime inversion (Jaiswal et al. 2008), joint inversion of magnetotellurics and seismic data (Liu et al. 2018), diffraction imaging (Schmelzbach et al. 2008), and the use of gravity data (Singh & Pandey 2007). Last but not least are interpretation-type methods which are helpful to reduce conceptual uncertainty: restoration and balancing techniques (Malz et al. 2015) or the exploitation of well-exposed outcrop examples and analog models of compressional structures (Serra 2018).

THE GEOLOGICAL SETTINGS
The study area is located in south-eastern Poland, at the border of two major geological units: the Carpathian Foredeep and the Carpathian orogenic wedge (Żytko et al. 1989, Ślączka et al. 2006) (Fig. 1). Over the years, a significant number of analytical studies have been performed to fully understand both the evolution and the structure of the Carpathian thrust belt (for summary and further references see e.g. Ślączka et al. 2006). Until 1996, when the first 3D seismic survey was acquired in the discussed zone, the recognition of the frontal Carpathians was mostly based on well data and outcrops' analyses, supported by the interpretation of 2D seismic data. Despite the lack of three-dimensional seismic images of this segment of the frontal part of the Carpathians, various types of studies were completed, and their contribution towards an understanding of the geological structure and evolution is undisputable.
The Outer Carpathians are built of a stack of thrust sheets, showing different lithostratigraphy and tectonic structures (cf. Ślączka et al. 2006). All Outer Carpathian thrust sheets (nappes) are thrust over the Miocene fill of the Carpathian Foredeep that is regarded as a typical peripheral foreland basin (Krzywiec 1999. The Skole Unit (Skole Nappe) covers the majority of the study area ( Fig. 1). It occupies the outermost position in the Polish Outer Carpathians and widens towards the east (Oszczypko 2006). This nappe is mostly composed of the Lower Cretaceous to Lower Miocene deep-sea sediments accumulated in the Skole Basin. In the discussed part of the Skole Unit, the following formations are present: the Krosno Beds, the Menilite Beds, the Hieroglyphic Beds, the Variegated Shales, the Ropianka Formation, along with the Bełwin Mudstones, the Dołhe Formation and the Spas Shales (Kotlarczyk 1978, Malata 1996, Waśkowska et al. 2019). Among them, the one with the highest thickness, reaching 2400 m, is the Ropianka Formation. It is represented by Turonian siliceous marls (with their thickness of 100-150 m) and the overlying sequence of Santonian-Paleocene sandstones and shales. Within the latter, one can distinguish: the Cisowa Member, the Wiar Member, the Leszczyny Member, and the Wola Korzeniecka Member (Kotlarczyk 1978). The total thickness of the above sandy-shaly series is as high as 2250 m.  Żytko et al. 1989, modified) North East from the Skole Nappe there is narrow, tectonic unit built of deformed foredeep deposits i.e. the Stebnik Nappe (Ney 1969, Oszczypko 2006. It is composed of mostly terrigenous sediments, formed in the internal part of the Carpathian Foredeep. The thickness of the considered fragment of the orogenic wedge varies from 0 to 100 m and generally increases towards SE. The basal thrust of the Stebnik Nappe, known as Stebnik Detachment, is formed along the Lower Miocene salt-bearing Carpathian formation (Oszczypko et al. 2008). Further towards the west from the study area, the Stebnik Unit is thrust over the Zgłobice thrust sheet (Fig. 1, Krzywiec et al. 2014), however, in the study area, no indication of the latter unit was found in the outcrops, deep wells and in the seismic data.

THE ROLE OF MODERN SEISMIC DATA IN RECOGNITION OF THE CARPATHIAN OROGENIC WEDGE
Over the last two decades the geological structure of this part of the Carpathian thrust belt has been recognised by several 3D seismic surveys with various shot-receiver configurations, along with different source parameters. In 2019 these data sets were merged with the newly-acquired seismic volume and cohesively reprocessed. The total binning area of the considered survey is 886 km 2 , out of which approximately 650 km 2 were subject to detailed structural interpretation. This new data set is the largest onshore 3D seismic volume that has ever been processed and interpreted in Poland. The ultimate goal of the above project was to determine the internal structure of the Miocene sediments and evaluate the extent of the Miocene clastic reservoir beneath the several-kilometres thick frontal part of the orogenic wedge. Deposits of autochthonous Miocene in this region are prospective for finding gas accumulations.
The technological improvement of reflection seismic method and corresponding seismic imaging enhancement is clearly illustrated by comparing the legacy 2D data and newly-processed 3D data (Fig. 2). The line was acquired in 1994 by utilizing a dynamite source, and in terms of its overall quality, is an average one among the archive data. The other parameters are representative to 1990's methodology of seismic data acquisition: the nominal fold is 60, the maximum offset is 2400 m, the shot interval is 40 m, and the receiver interval equals 20 m.
A basic and somewhat subjective structural interpretation can be made on the basis of the legacy data ( Fig. 2A). One may notice an unclear frontal Carpathian thrust, separating two different parts of the image: one with the steeply-dipping reflections (upper left) and one where the dominating dip is not present (lower right). The events observed within the orogenic wedge, which represent the true dip of the allochthonous strata, are highly-mixed with various out-of-plane waveforms. These are the noise and seismic artifacts caused by the limitations of the 2D time-domain migration to account for the complexity of the existing geological structure. The individual thrusts can only be interpreted as indicative ones, due to low variability of seismic amplitudes within the Skole Unit and unclear angular unconformities. No specific seismic signature can be attributed to the Stebnik Nappe, hence it cannot be distinguished from the overlying formations.
Apart from the deepest infill of the paleovalley located beneath the orogenic wedge, the reflection alignment within the autochthonous Miocene deposits to a considerable degree departs from the correct one. Beneath the thick cover of the orogenic wedge, seismic response of the Precambrian basement is weak and disturbed. The seismic event associated with the top of the basement can only be interpreted locally. There is virtually no information regarding the true internal structure of the low-thickness autochthonous foredeep Miocene underlying the thickest part of the Carpathians. The ambiguity level of the geological interpretation, based upon 2D legacy data, is undoubtedly high.
In Figure 2B one can notice significantly higher reflectivity, higher S/N ratio, improved definition of the thrust sheets and minimum number of migration artefacts. Well-preserved morphology of the basement, along with the interpretable infill of the paleovalley, together with high continuity and correct dips of shallower Miocene horizons, are observed on the modern newly-processed 3D data set.
Compared to the 2D data, more conclusions can be drawn from the legacy 3D seismic data, recorded in 2007 (Fig. 3A). Some elements of the geological structure were even imaged more precisely than on the new volume (Fig. 3B). The past acquisition and processing methodology (e.g. the common use of Spectral Whitening applied to pre-stack data) allowed for obtaining high vertical resolution and well-defined thrusts within the Carpathian wedge.
In the case of the newly-processed data, the biggest asset is more regular offset-azimuthal distribution that is reflected by spatially-uniform data (Fig. 3B). The imaging of the frontal thrust, along with the geometry and amplitude of the Precambrian seismic reflection, is also unquestionably better on the current volume. The new data set is also free from the unprocessed noise and migration artifacts which dominate the NE part of the archive one, obscuring the true alignment of the Carpathian Foredeep horizons. In addition, the older vintage does not properly show the steepest reflections and faults, even in the presence of significant throws. Due to the above reasons, the interpretation of the 3D legacy data still poses the risk of exploration failure.
An example of structural interpretation performed with the newly-acquired 3D seismic data is presented in Figure 4. One can notice a reliable image of the frontal thrust, the internal Carpathian thrusts and associated sheets, the top  of the basement and the autochthonous Miocene sediments. Large tectonic zones, such as the Kniażyce fault, and many other, low-throw faults, are well-visible on the new data. By means of the new seismic volume one can directly determine the laterally-varying thickness of the Carpathian Foredeep autochthonous deposits that ranges from virtually 0 m (SW) to over 3500 m (NE). Due to the significant seismic wave attenuation and substantial burial depth, there is no evidence that the Carpathians rest directly on the Precambrian basement in the southernmost part of the area. The thickness of undeformed Miocene sediments can be considered minimal and below vertical resolution of the seismic data. Apart from the above general structural features, the modern image reveals plenty of details, a selection of which is analysed and described below (Fig. 5).

Fig. 4. Interpreted time-domain seismic section from the newly-acquired 3D seismic volume (PreSTM)
latter case, the wavelet extraction from current data set shows that the tuning thickness varies from 13-19 m in the Carpathian Foredeep area, up to 33-47 m within an infill of the paleovalleys, deeply buried (>3 km) below the orogenic wedge. At the same time, analysis of the seismic image of the Skole Unit indicates that the layers with the thickness of approx. 10 m are distinguishable and can be unmistakably correlated with the well data (Fig. 5B). Actual vertical resolution related to the frontal thrust cannot be determined because it is strongly depth-and dip-dependent.
The reflectivity corresponding to the Carpathian flysch sediments is not only thickness-dependent, but also lithology-dependent. Some strata, such as the Dołhe Formation, usually exhibit high reflectivity, caused by the presence of cluster of shales within the sequence dominated by the sandstones and mudstones. That feature, in turn, allows for an unambiguous interpretation to be performed (Fig. 5C).
In practice, even if the new seismic data are used, the frontal thrust can only be precisely inspected at the following depth interval: 0.5-2.5 km. At such depths, the thrust is clearly visible as a sharp, angular contact between the Stebnik Unit and undeformed Carpathian Foredeep infill. In the deeper part, the thrust surface becomes relatively flat and resembles flat-lying Miocene autochthonous strata. Among other factors which obstruct the exact tracking of the frontal Carpathian thrust, at least two should be listed: an identical lithological composition of the Stebnik Nappe and underlying sub-thrust foredeep infill that implies zero impedance contrast at their contact (Fig. 5A), and limited vertical resolution of the seismic data.
Considering the latter aspect, due to a joint effect of seismic wave attenuation and sediment burial depth, vertical seismic resolution within the flysch formations of the Outer Carpathians is higher than in the autochthonous strata. In the The seismic response of the individual thrusts depends on the internal composition of the surrounding sediments. In some relatively rare cases, the thrust surface imitates a typical seismic horizon, and therefore, can be picked along the specified high-amplitude event (Fig. 5D, bright green). Particularly essential is the fact that such surface can be further utilized in the process of PreSDM velocity model building. The more common and distinct features of the thrusts, such as the dip changes and angular unconformities, are clearly observable within the new data set (Fig. 5D, red rectangle).
Interpretation of the internal structure of individual sheets remains a nontrivial task, not only due to high degree of tectonic deformation, but also due to the low acoustic impedance difference at between particular lithological formations of the Outer Carpathians. In contrast to sandy-shaly interfaces, at the top and base of evaporites, such as anhydrite and salt, a large reflection coefficient is present. Analysis of the local well logging data revealed that acoustic impedance within the evaporites is considerably higher than in surrounding flysch. Consequently, the seismic signature of these rocks is a positive peak at their top and negative trough at their base (Fig. 5E). The existence of evaporites within the study area, both allochthonous and autochthonous ones, as well as their origin, is well-documented (Garlicki 1979, Czapowski 1994). In the analysed 3D seismic data these rocks are visible as high-amplitude events within the Stebnik Nappe (Fig. 5E, F). The stiff anhydrites, stuck within the ductile flysch deposits, and cut by a multitude of faults, show a largescale deformation that takes place under the burden of the Carpathian wedge. The presence of the allochthonous evaporites is a distinctive feature of the Stebnik Unit, and they should be regarded as a marker horizon used during the interpretation of the frontal Carpathian thrust.
Besides the above mentioned phenomena, the modern seismic volume allows for accurate mapping of the Precambrian basement, verifying known tectonic zones and delineating the correct geometry of already recognized structural traps. Additionally, some new structural elements, interesting from the exploration point of view, could be discerned on described data set.

DATA AND METHODS
In this project, the following types of data were utilized: 3D seismic data set (merged and processed with neighbouring legacy volumes), surface geologic maps, well stratigraphy, uphole data, VSP, formation evaluation results, sonic and density logs. Presented analysis is an extension of standard exploration scheme that includes: seismic data acquisition, processing and interpretation. Table 1 contains the main parameters of this new seismic survey. The processing consisted of the following key stages: signal matching between the newly-acquired and archive volumes, noise removal, surface-consistent amplitude scaling, three iterations of stacking velocity analysis and statics, External Model Trim Statics, 5D interpolation, followed by Kirchhoff VTI pre-stack time migration. The above standard processing route was further followed by the creation of PreSDM initial velocity field. During that stage, many questions arose regarding the correct approach to near-surface velocity model building. Initial PreSDM velocity model was built upon utilizing major seismic boundaries, interpreted on the time-migrated section: the Precambrian basement, the frontal thrust of Carpathians, intra-Carpathian thrusts and several horizons representing the autochthonous Miocene deposits. One of the main challenges was to assign correct velocities for the Skole Nappe. Because of its complex internal structure, it had to be subdivided into several smaller-scale sheets using both seismic data and sonic logs. Within each sheet, a corresponding fragment of the considered log was found in all wells located in the seismic survey area and generalized to establish one recurring (velocity versus depth) trend. These trends were further spread into the seismic volume by means of interpreted horizons (Fig. 6).
The main geological target of the seismic survey was a tight gas reservoir, buried below the Carpathians at an approximate depth of 3 km. The goal was achieved, but the processing did not ensure satisfactory visibility of the shallowest events (up to 700 m below the ground). In practice, imaging of this portion of the subsurface is related to largest minimum offset, equal to the length of the diagonal of the box, formed by two consecutive receiver lines and two neighbouring source lines.
This parameter reached 492 m, which admittedly, is a reasonable result for land seismic data, but still not as good as in modern onshore high-density seismic surveys, where it ranges from 160 to 300 m (Liu et al. 2009). The potential improvement of this parameter would result in an unjustifiable rise in acquisition costs. The shallow velocity model remained unrecognised also due to the lack of reliable sonic logs. Because of the deficiency of the well data, the velocity field from 3D tomographic inversion of the first breaks was taken into account in the process of PreSDM model building. The tomography method is an industry-standard one, commonly applied during the calculation of time-domain statics (Li 1996, Zhu et al. 1998, Zhu 1999. The tomography results were further merged along the pseudo-datum with sonic log velocities (Fig. 6). Due to the fact that these two types of velocities occurred to be markedly different, the analysis below, utilizing finite difference modelling and inversion of the synthetic data, was carried out.
Finite-difference (FD) methods are precise and robust for solving the wave equation that is widely utilized in reflection seismic (Ikelle & Amundsen 2018). FD equations play a prominent role in both industrial and academic computational seismology, especially with respect to their application in synthetic seismogram calculation and migration (Ristow & Ruhl 1994). Their advantage over other computational methods is the ability to entirely describe wave motion in geological media with

Fig. 6. Time domain comparison of the velocities derived from 3D tomographic inversion of the refracted waves (above the pseudo-datum) and well-log-based velocities (below the pseudo-datum)
any variation of the elastic properties. Apart from high costs and high hardware requirements, major drawback of the discussed techniques is that they may be considered a black-box (Lecomte et al. 2015). The user sets up the parameters and when the process is finished, plenty of various wave types appear, including: the head waves, surface waves, diffractions and free-surface effects. For that reason, analysis of the output wavefield may require a lot of expertise and effort.
The FD algorithm utilized in this study is an integral part of the Promax/SeisSpace processing software (Vidale 1990). The module computes the synthetic data in common shot domain. A 2D acoustic approach requires velocity and density attributes as input to generate the seismograms. Two 2D velocity models of the uppermost part of the Carpathian orogenic wedge were utilized to perform the simulations (Fig. 7A, B). Each of them contains 1400 traces, spaced out by 1 metre in the horizontal direction, with the same depth sampling rate. The choice of modelling parameters was made to obtain suitably high amplitude of the refracted wave, sufficient vertical resolution and good separation between the direct wave and the refracted wave (Tab. 2). In order to find out if refracted wave tomography provides the correct solution in the shallow section of the Carpathian fold and thrust belt, the modelled first arrivals (Figs. 8,9) were picked and inverted in the chosen ZondST2D commercial software (Kaminsky 2020). Refraction tomography method is believed to resolve both vertical and lateral velocity changes and may be effectively applied in complex geological settings, i.e. the areas of compaction, dipping bedrock and fault zones. Nonetheless, we should not forget that each commercially available algorithm has its own theoretical assumptions and limitations (Sheehan et al. 2005). In the method of refracted wave tomography, the subsurface medium is divided into a set of cells with known size, usually increasing with depth. After setting the background model, the picked first break times are used to iteratively update the model, simultaneously the difference between predicted and actual travel times, is minimized. The model refinement may be performed on the basis of various solutions, through: 2-D Waveform Eikonal Travel-time WET (Schuster & Quintus-Bosz 1993), Occam inversion (Constable et al. 1987) or many other algorithms. Among the available computation techniques, the Occam inversion was selected as an optimal one to be used for the studied synthetic data. This method was originally invented for electromagnetic sounding data (Constable et al. 1987), but can be used to solve refraction tomography as well. This algorithm can be considered the least squares process that uses a smoothing operator and additional contrast minimization. The initial model should be kept as simple as possible. Rather than matching the first arrivals accurately, a smooth solution that fits the data within an expected tolerance, is found. The inversion, run on the examined data set, turned out to be stable, efficient, and typically converged in no more than 10 iterations.

Near-surface velocity model
The challenge in recovering correct P-wave velocities of the shallow medium can be addressed by showing the statistical summary of the uphole velocity and the velocity derived from first-arrival refraction tomography on newly-processed 3D seismic data (Fig. 10). The same depth range, reaching 50-90 m from the ground elevation, is taken as an input for these statistics. The uphole data are represented by the maximum velocity encountered at the measurement point, and in the case of the tomography, each 5 m depth sample of the resulting volume, is taken into account. What is particularly important, the presented histograms are calculated separately for two major regions: the Carpathian orogenic wedge and the Carpathian Foredeep, both limited laterally by an extent of the current seismic data. As it can be noticed by comparing the upper and lower part of Figure 10, the tomography-based velocities in the Carpathian Foredeep region are comparable to the uphole ones. The mean values are respectively 2083 and 2043 m/s, and both are in excellent agreement with the expected values within the shallow undeformed Miocene sediments. Their consistency with the velocity trends from the sonic logs can also be concluded from upper right part of Figure 6 (above pseudo-datum). In case of the Carpathian wedge area, one can notice a broader range of values derived from tomography volume. A significant number of samples are far beyond 2000 m/s, with the mean value of 2475 m/s. It is definitely not in accordance with the uphole data, which represent the average value of 1844 m/s. As a consequence, the difference between the results obtained using these two methods, is equal to 34%. From the presented calculations an immediate conclusion can be drawn that under favourable conditions of a near-surface medium, various methods of assessing seismic wave velocity are in agreement with each other, but in the examined part of the Carpathians, they are definitely not.
The statistical analysis of existing uphole data was one of the most essential stages of this project. The total number of accessible uphole locations exceeds 1000. Instead of a standard practice that amounts to data interpolation between the measurement points, these points were projected onto the surface geologic map and assigned to individual formations (Fig. 11). In order to fully understand the complexity of existing near-surface medium one can focus on the selected fragment of the discussed map. It shows various formations underlying the weathering zone, which are distinguished in prolific outcrops and the shallow wells (Szotek et al. 2018). This map informs us of the highly-varying spatial distribution of the formations belonging to the Skole Unit. Subsequently, for each investigated unit, separate statistics of velocities and thicknesses were generated (Tab. 3). Fig. 11. Uphole locations superimposed on the fragment of the surface geological map (after Szotek et al. 2018, modified). The map show the concept of the uphole velocity assignment to individual formations, observed on the ground surface The maximum velocity registered in each uphole was attributed to the bedrock layer (i.e. the one that underlies the weathering zone) and averaged within individual unit. As it can be noticed from Table 3, there are significant differences between the observed formation velocities. The lowest value was found in the Krakowiec Clays (V = 1528 m/s), whereas the highest one, was computed in Cretaceous sediments (the Fucoid Marl, Wiar Member of the Ropianka Formation, V = 2290 m/s). Two main causes of the velocity variability were noticed: the age of the sediments (Miocene strata usually have lower velocity than Cretaceous ones) and lithology -sandstones and marls exhibit higher velocities than shales.
The presented analysis revealed that the P-wave velocity in the near-surface geological medium is varying and formation-dependent. The Krakowiec Clays represent a simple model with significant contrast at the base of the weathering zone (Fig. 12A). Within this formation, in more than 80% of the examined uphole locations, two-layered velocity model, consisting of the weathered and the bedrock layer, was interpreted, whereas only in 5% of the measurement points, 5 or more layers were distinguished (Tab. 3, 4 th column). As opposed to the Skole Nappe formations, the Krakowiec Clays bedrock is approximately flat, therefore, we are dealing with favourable conditions to generate the refracted waves. In such situation the output of the first-arrival inversion must be realistic. On the other hand, interval velocity computed in the shallow lithological formations of the Skole Unit (e.g. the Ropianka Formation -undivided), do not show a single refractor (Fig. 12B). Many of them (the Menilite Beds, the Spas Shales, the Upper Krosno Beds, the Wiar Member of the Ropianka Formation) can be regarded as a gradient medium, built of several distinct layers (Tab. 3, 4 th column). Their bedrock is definitely untypical, with the horizon dips exceeding 70 degrees, which cannot be considered an ideal medium for generating the refracted waves. The total depth of the upholes range from 50 to 90 m, and just below that depth, the prediction of the velocity field is imprecise. The sonic logs frequently do not cover the near-surface well section (several hundred metres from the ground), whereas check shot interval velocities are prone to errors, due to the requirement to correct the first-arrival times for non-vertical wave propagation. Despite the fact that a large number of wells have been drilled within the study area, the majority of them are of no little or no use, as they contain low-reliability velocity data. Among the investigated wells, a new one drilled in 2011 was used to determine the average velocity gradient in the depth range: LVL base -700 m. The resulting velocity gradient is equal to 0.71 m/s per metre, which seems to be considerably lower comparing to other geological regions in Poland. It could be regarded as evidence of untypical feature of the Carpathian wedge with respect to velocity field and seismic wave propagation.
On the basis of the upholes, topography, surface geologic map and interpretation of the 3D seismic data, several models which represent the shallow part of Carpathian fold and thrust belt, were built. During the model building process, it was necessary to extend the dipping seismic horizons and thrust surfaces visible at greater depths up to the ground and link them to the geologic map.
An example of 1.4 km long and 700 m deep, two-dimensional SW-NE model that crosses the Skole Unit, is shown in Figure 7. Note the dipping bedrock layers and variable, formation-dependent LVL thicknesses, derived directly from the uphole data. Model A presents a simplified situation, in which the velocity is constant within the bedrock formation, whereas in model B, the computed velocity gradient of 0.71 is assigned to each unit. It is evident that the gradient may vary to some extent between particular thrust sheets, however, one must take into account that there is not enough reliable data to compute representative gradients for each formation. The presented models should be considered a good starting point for any simulation of wave propagation in the examined medium.

Finite-difference modelling
An analysis of the modelling results revealed major features of the seismic wavefield generated in the shallow strata (Figs. 8,9). The refracted wave arrivals are evident, visible mainly as positive peaks preceded by the side lobes. Despite utilizing the theoretical Ricker wavelet, the shape of the refracted wave deviates from the zero-phase to a certain extent, and evidently depends on the geological model of the shallow zone. At the locations where the depth of the high-velocity bedrock formation is minimal, the energy of the upper side lobe tends to increase. The arrangement and position of the first arrival times are consistent with the geometry of the low-velocity zone. The cross-over point where the direct and the refracted wave meet is not a single point in space but rather an interference zone, where a sudden change of the first arrivals occurs and which might lead to some confusion during the picking process. As it can be noticed from Figure 7, at the base of the weathering zone and also at the contact of individual sheets, rapid velocity changes are introduced into the model. These changes are the main reason for the generation of diffracted waves. These waves interfere with the adjacent refracted waves, but in most cases do not hinder the interpretation of the first arrivals.
In industry refraction surveys some shot points could also be placed in the areas of low (2-3 m) weathering zone thickness. In such cases, the distinction of the refracted waves may not be possible on local fragments of the receiver array due to the presence of the shadow zones. The contacts of individual sheets also have their representation on the seismograms. Because of the assumed velocity differences between the formations, these contacts are equivalent to typical seismic horizons. At least two sheets' contacts are noticeable on the synthetic data as rough, dipping events, truncated at the first arrival line. These reflections can be regarded as a seismic signature of the thrusts from the right end of the model. The generated signal propagates through the LVL, and for this reason, the shape of the weathered strata is superimposed on the thrust signatures. Consequently, in order to obtain the correct image of individual sheets and their contacts, the true LVL model needs to be determined at first.
The overall difference between the synthetic responses, generated for constant (Fig. 8) and gradient velocity model (Fig. 9), is barely detectable. One can notice only minor changes in the amplitude of the refracted and upgoing waves, and variations in their time position. Hence, the differentiation between these two model types, and computation of the velocity gradient is a tough assignment.

Refracted wave tomography
The verification of the tomography solution, using the modelled and interpreted first arrivals, was performed in dedicated ZondST2D commercial software (Kaminsky 2020). The results obtained for the discussed velocity models are shown in Figure 13. In both cases the overall shape of LVL zone is recovered correctly and includes its features, such as local thickening of the weathered layers. It can be noticed that the capability of the tomography to recreate the formation velocity is considerably better where the topography is relatively uncomplicated and flat. In such conditions it may recover even a minor thickening of shallowest layer, as it can be observed at the right end of the line. Conversely, in the presence of the ground slopes, the inversion process tends to underestimate the velocity and introduce velocity gradient at sharp formation contacts.
A clear advantage of the Occam algorithm is that, despite using an imperfect initial model, the process converges to close-to-true velocities in the bedrock formations. The background model does not contain any details about the dips and velocities of the sub-LVL formations. As can be seen in Figure 13, after the inversion, thrust sheet with the highest velocity (i.e. Ropianka Fm., traces: 910-1140, V = 2060 m/s) could easily be distinguished from the surrounding ones (i.e. Variegated Shales, V = 1757 m/s). The computation output is fully in accordance with the hallmark of the Occam algorithm, because the final result seems to be a smooth version of real geological structure. The inversion also tends to overestimate the velocities in the central part of the high-velocity formations, however, the difference between the set and recovered values doesn't exceed 300 m/s. The difference in the obtained tomography image between the constant and gradient-velocity models is insignificant. In the latter case, the rays propagate to greater depths, therefore, especially in the sloping terrain part of the model, the inverted section is more reliable. The tomography did not recover the assumed velocity change with depth. It is caused by the limited number of rays penetrating the bedrock. In practical seismic surveys, performed in the areas where sub-LVL strata exhibit low velocity gradient, the recognition of the gradient seems to be on the verge of feasibility.

SUMMARY AND CONCLUSIONS
Newly-acquired 3D seismic data play a major role in the process of hydrocarbon exploration in the frontal Carpathians in SE Poland. These data show the complexity of the fold and thrust belt in much greater details than the legacy seismic data, both 2D and 3D. The current methodology of planning, acquisition and processing is suitable for imaging the frontal Carpathian thrust, together with its internal thrusts and sheets. One can notice a vast progress that took place in seismic data technology during last two decades. At the end of the processing, these data seems to be ready to be exploited by geologists, however, we should not forget about the limitations related to acquisition costs and processing algorithms, together with the common deficiency of the well data constraints.
Despite the high technological improvement of the seismic reflection method, there are still unsolved issues. The major problem encountered in the study area is the shallow part, with sparse reflections and difficult-to-predict velocity model.
In the presented work, the near-surface velocity model was examined by means of various geophysical methods. The analysis revealed a considerable difference between the results of: 3D refraction tomography solution, the sonic logs and the upholes. Such a discrepancy is, however, limited to an extent of investigated part of the Carpathian wedge.
In the attempt to solve this problem, a different approach to the process of velocity model building is proposed. Instead of using the 3D tomography-based volume, the shallow model is composed of interpreted thrusts and sheets, with assigned uphole-based velocities. Such model cannot be used during standard time-domain processing, since it utilizes a detailed interpretation of the pre-stack migrated horizons, together with the thrusts. Nevertheless, it is a good starting point for other methods, such as: seismic survey design and forward modelling.
In the presented work, in order to verify the industry-standard refraction tomography solution,seismic modelling together with 2D tomographic inversion in selected software, was performed. Synthetic seismograms show the high degree of complexity of the wavefield in the near-surface fragment of the Carpathian orogenic wedge. The inversion, applied to the synthetic data, converges to close-to-assumed velocities in the steeply-dipping formations underlying the weathering zone. The presented algorithm needs to be further validated on the basis of the measured seismic data.
Some state-of-the-art methods of the modelling, processing, and interpretation, developed with the aim of improving seismic imaging in the fold and thrust areas, are also demonstrated. These are the potential directions of the development of the seismic technology for the Polish reflection seismic exploration in the Carpathian region.
The author wishes to acknowledge the help of the reviewers of the presented study, especially Piotr Krzywiec from the Polish Academy of Sciences, for their kind and willing cooperation, useful suggestions, valuable feedback and improving the quality of the publication. The author would like to thank POGC, Warsaw, Poland for supplying access to their data.