Acoustic finite difference method (AFDM)
The AFDM is one of the first numerical methods used in the solution of differential wave equations and it has been widely used by many researchers in the calculation of seismic shot data. The method is fully applicable to numerical solutions of acoustic and elastic wave equations (Moczo et al. 2007; Margrave 2003; Carcione et al. 2002; Youzwishen and Margrave 1999; Manning and Margrave 1998; Krebes and Lee 1994; Marfurt 1984).
The AFDM, provides the direct solution of the partial differential acoustic wave equations given by Eq. (1) under certain initial and boundary conditions. The method provides the possibility to obtain the propagation of acoustic waves from simple to complex geological models. The 2D scalar (numerical) wave equation (x, z) of the method can be written in Cartesian coordinates as follows (Lines et al. 1999).
$$\frac{{{\partial ^2}\phi (x,z,t)}}{{\partial {t^2}}}={v^2}(x,z){\nabla ^2}\phi (x,z,t)$$
1
In Eq. 1, x and z are the horizontal and vertical distances (metres-m) on the grid network, respectively. T is the two-way travel time (seconds-s), v is the velocity of the medium in which the wave propagates (metres/second-m/s) and \({\nabla ^2}\)\(\phi (x,z,t)\)represents the wave potential or acoustic wave field. The Laplacian operator \({\nabla ^2}\)is given by Eq. (2).
$${\nabla ^2}\phi =\frac{{{\partial ^2}\phi }}{{\partial {x^2}}}+\frac{{{\partial ^2}\phi }}{{\partial {z^2}}}$$
2
The Laplacian operator can be approximated with central difference operators. The two approximations used within the AFD software are a second and a fourth order approximation. The approximations use five and nine points of the grid respectively. In this study, Eq. (3), which is a fourth order Laplace operator approximation with nine grid points, is used because it provides high accuracy and broadband solution, although it increases the computation time (Youzwishen and Margrave 1999).
$$\begin{gathered} {\nabla ^2}\phi \left( {x,z} \right) \approx \frac{{ - \phi \left( {x+2,z} \right)+16\phi \left( {x+1,z} \right) - 30\phi \left( {x,z} \right)+16\phi \left( {x - 1,z} \right) - \phi \left( {x - 2,z} \right)}}{{12\Delta {x^2}}}+ \hfill \\ \frac{{ - \phi \left( {x,z+2} \right)+16\phi \left( {x,z+1} \right) - 30\phi \left( {x,z} \right)+16\phi \left( {x,z - 1} \right) - \phi \left( {x,z - 2} \right)}}{{12\Delta {z^2}}} \hfill \\ \end{gathered}$$
3
The seismic wave field is calculated by iteratively solving with FDM at each mesh point shown in Fig. 1. However, there are few analytical solutions of the wave equations and approximate numerical solutions are widely preferred for wave field modelling.
Finite-difference approximation is a numerical process that is formally composed of constants and variables (Fig. 1). As shown in Fig. 1, any point (i, j) is called a grid point. The value of the dependent variable at one point in the grid region can be related to the values at neighbouring points. For example, the points (i + 1, j) and (i-1, j) are located to the right and left, while the points (i, j + 1) and (i, j-1) represent the points located above and below. For a grid mesh as shown in Fig. 1, the grid points are given by equations (4) and (5) below. If the mesh used is a regular grid, grid points can be easily expressed. For example, for a grid such as in Fig. 1, the mesh points are as follows,
\({x_i}=i.\Delta {x_i}\) \(i=0,1,2,...,M\) (4)
\({z_i}=j.\Delta {z_i}\) \(j=0,1,2,...,N\) (5)
Here, Δt: temporal sampling interval (sec-s), Δx: spatial sampling interval (meter-m) (Δx = Δz). The calculation of the time derivative on the left hand side of Eq. (1) using the second order central difference approach is given in Eq. 5.
$$\frac{{{\partial ^2}\phi }}{{\partial {t^2}}}=\frac{{\phi (t+\Delta t) - 2\phi \left( t \right)+\phi (t - \Delta t)}}{{\Delta {t^2}}}$$
6
If equations (3) and (6) are substituted into the numerical wave equation in Eq. (1), the wave field at time t + Δt can be solved iteratively (Eq. 7).
$$\phi (x,z,t+\Delta t) \approx (2+\Delta {t^2}\,V{\left( {x,z} \right)^2}{\nabla ^2})\phi (x,z,t) - \phi (x,z,t - \Delta t)$$
7
In Eq. 7, \(\phi\) is the P-wave potential. Eq. 7 demonstrates that if the wave field is known at time t and t-Δt, it is possible to calculate the wave field at time t + Δt. This process is called a snapshot or time stepping of the wave propagation at each time. The wavefield is extracted at time t-∆t and the Laplacian operator is applied to the wavefield at time t.
Applications
In this study, geometrical (stratified structure, discontinuity) and physical properties (velocity and density information) of simple and complex trap models were determined from the literature (Kearey and Brooks 1984; Mavko et al. 1998). The source function, spatial and temporal calculation parameters are given in Table 1. The variation in the volumetric density values of the layers was considered constant (ρ = 2.0 g/cm3), as it was significantly small compared to the seismic wave velocity. For modelling, the reflective interfaces were digitised and the depth, distance and velocity information were introduced into the modelling software.
Table 1
Parameters used for modelling
Modelling Parameters | Channel | Granit Wash | Normal Fault | Bioherm | Anticlinal |
Profile length (m) | 1000 | 2000 | 2000 | 2000 | 1000 |
Maximum depth (m) | 1000 | 1000 | 1000 | 1000 | 700 |
Receiver interval (m) | 10 | 10 | 10 | 10 | 10 |
Shot interval (m) | 80 | 40 | 40 | 40 | 40 |
Number of shots | 25 | 46 | 25 | 25 | 25 |
Number of receiver | 101 | 201 | 201 | 201 | 101 |
Maximum velocity (m/s) | 4000 | 4000 | 4000 | 4000 | 4000 |
Minimum velocity (m/s) | 2000 | 2000 | 2000 | 2000 | 2000 |
Maximum offset (m) | 2000 | 2000 | 2000 | 2000 | 1000 |
Minimum ofset (m) | 100 | 100 | 100 | 100 | 100 |
Calculation time step (ms) | 0.02 | 0.02 | 0.02 | 0.02 | 0.02 |
Sampling time (ms) | 4 | 4 | 4 | 4 | 4 |
Record lenght (ms) | 1000 | 1000 | 1000 | 1000 | 1000 |
Minimum Phase Ricker Source Wavelet (Hz) | 30 | 30 | 30 | 30 | 30 |
Trap models are constructed with layers by defining layers with different P-wave velocities as polygonal regions. Initially, the entire area is gridded with horizontal and vertical distances of 10 meters each. This grid mesh has a total of 1000 nodes of 10 x 10 m and is evenly distributed. Accordingly, the number of nodes in x and z directions of this grid mesh is 100 and the grid distance in x and z directions is 10 m. The total number of nodes is 10,000. After the gridding process, the velocity matrix is determined. The minimum and maximum velocity are assigned. Then, the corner coordinates (x, z) of the polygonal region are defined. Velocity models are created by superimposing several polygonal regions. To calculate the shot records of the obtained velocity models, a series of shot records are generated by defining the seismogram sampling interval (dt), time calculation step (dtstep) and recording time (tmax) values. Later, the modeled shot data are written to a SEGY file (Margrave 2000). An example grid mesh illustrating the synthetic model containing source-receiver locations for modeling is shown in Fig. 2.
Channel Trap Model
Channel structures are considered to be one of the best stratigraphic traps for hydrocarbon accumulation (El-Nikhely et al. 2022). It has been determined that the majority of oil and gas fields in the world are associated with channel-type structures (Mohebian et al. 2018). Channel models are generally 100–900 meters wide and characterised as less meandering structures. Therefore, channel systems are considered as traps with a high probability of containing hydrocarbons (Shuxin et al. 2017). Moreover, although the channel structures are interrupted by unconformities since they are porous and permeable, these channels can be filled with hydrocarbons and form high-quality reservoirs. The channel reservoir model developed in a regularly stacked horizontal stratified environment is shown in Fig. 3a. The channel structure in the model is designed to be 500 meters deep, 30 meters thick, 50 meters wide, porous and permeable sandstone containing natural gas with a velocity of 1400 meters per second. The channel is surrounded by imporous and impermeable marl stone with a velocity of 3000 metres/second. The model consists of alluvium at Z = 0-200 meters, Vp = 1500 m/sec, clay at Z = 200–250 meters, Vp = 1700 m/sec, porous sandstone at Z = 250–350 meters, Vp = 2000 m/sec, saturated shale at Z = 350–450 m, Vp = 2100 m/sec, limestone at Z = 450–500 meters, Vp = 2500 m/sec and dolomite at Z = 700–1000 meters, Vp = 3500 m/sec at the bottom.
Granit Wash Trap Model
The granite wash trap model, which is a stratigraphic trap type, has been the target of petroleum research and development studies with the discovery of petroleum (Sproule 1956). A significant amount of hydrocarbons in granite wash traps accumulate in low permeability reservoirs, unlike conventional medium high-permeability oil traps. These hydrocarbon resources are referred to as unconventional resources and are known as stratigraphically tight gas-sand reservoirs according to the classification of Hyne (1984). These type of traps shown in Fig. 3b are sandstones associated with decomposed granite basement rock and are deep structures (Dec et al. 1996). Underground geological
structures containing such hydrocarbon traps are quite complex and do not have similar characteristics, so it can be quite difficult to identify and define them by seismic surveys.
Figure 3b shows a model illustrating a granite wash trap within a multilayered medium. The model is designed to encompass a range of geological formations, including porous and permeable sandstone, limestone, dolomite, or fractured rock, with a velocity unit of 2700 m/s beneath the granite rock cover. The area surrounding the peak point, extending 1200 m/s on each side from the apex, is depicted in a turquoise blue colour. However, the reservoir is completely covered by rocks with velocities of 2650 m/s (brown unit between 70–850 m), 2500 m/s (light green unit between 500–600 m) and 2580 m/s (yellow units between 400–500 m and 600–750 m; impermeable marl, shale, salt or micritic limestone). The other units were stacked in harmony with the cover rock (the yellow-coloured unit between 400–500 m) which has an anticlinal appearance towards the surface. The velocity of the surface layer is 1000 m/s and its thickness is 150 m on the left side, 50 m in the centre and 125 m on the right side.
Normal Fault Trap Model
Normal fault type traps are structural traps that are curved, formed as a result of the intersection of two faults or the intersection of many faults. Normal fault traps occur when the movement of fault blocks inhibits the migration of oil. For instance, on one side of the fault, an impermeable and sealed formation may shift in the opposite direction compared to the oil-bearing formation on the other side. In this case, the impermeable layer prevents the flow of oil and an oil or natural gas reservoir is formed against the fault (Biddle and Wielchowsky 1994). An example normal fault trap is shown in Fig. 3c. In this model, natural gas (1500 m/s - turquoise blue colour) is placed inside the trap structures and covered with overburden rock (2580 m/s - light green colour).
Bioherm Model
The bioherm trap model, a type of stratigraphic trap, consists of the skeletons of various marine organisms such as corals, spines and molluscs (Cuming and Shrock 1928). They are defined as mound and lens-shaped biological structures that are unconformably aligned within a stratigraphic series of different lithologies (Riding 2002). Bioherms, known as unconventional stratigraphic traps according to Hyne's (1984) classification, are also recognised as anticlinal-type non-structural traps (Sosnin et al. 2021). Bioherms constitute an important part of stratigraphic traps in hydrocarbon exploration and have great potential for oil and natural gas reservoir formation. Therefore, bioherms have become important targets for hydrocarbon exploration. They are characterized on seismic sections by chaotic bedding, intermittent internal reflections, chaotic or blank reflections, mounded reflections, and apparent amplitude anomalies (Sattler et al. 2004; Yu et al. 2005).
An example bioherm trap model is shown in Fig. 3d. In this model, a reservoir containing natural gas with an average seismic P-wave velocity of 1200 m/s is placed in the 380–510 m depth range. From the surface to the depth, it is covered by porous water-saturated sandstone with Vp = 2300 m/s (the green-coloured unit between 110–430 m), marl with Vp = 2600 m/s (the orange unit between 250–690 m) and saturated shale with Vp = 2700 m/s (the yellow coloured unit between 350–800 m). The velocity of the surface layer is 1800 m/s and its thickness is 225 m on the right side, 110 m in the middle and 150 m on the left side. The other units are stacked downwards as sandy limestone with Vp = 2700 m/s (the grey unit between 510–860 m), porous water-saturated sandstone with Vp = 2800 m/s (the green unit between 690–1000 m) and basalt with Vp = 2900 m/s (the dark red unit between 840–1000 m).
Anticline Model
The anticlinal trap model, a type of structural trap, has an arch-like shape and is a type of fold with the oldest rocks at its centre. They usually develop above thrust faults. Therefore, any small compression within the inner crust causes large impacts on the upper rock layer (Earle 2015).
Anticlines are one of the first known examples of oil traps. About 80 per cent of the world's oil is found in anticlinal traps. Within the anticlinal trap, the low density of the oil causes it to migrate towards the surface, floating in the source rock until it is deposited and stored in the reservoir rock, such as sandstones or porous limestones. Oil, along with water and natural gas, is held by a overburden rock consisting of an impermeable layer or an impermeable layer such as a fault zone (Riva 2015).
The anticlinal trap model in a partially fluid-saturated porous medium is shown in Fig. 3e. The model includes oil sands at Z = 620–720 m, Vp = 2580 m/s and gas sands at Z = 540–620 m, Vp = 2540 m/s above it. The oil sands reservoir is sealed by rocks with Z = 0-490 m, Vp = 2700 m/s shale, Z = 190–720 m, Vp = 2620 m/s water-sand and Z = 410–820 m, Vp = 2750 m/s. However, the trap model was designed with source rock with Z = 690–1000 m, Vp = 2800 m/s (dolomite or fractured rock, limestone, permeable and porous sandstone) and slurry-sand with Vp = 2620 m/s on both sides of the source rock (left side Z = 720–1000 m and right side Z = 720–910 m) (Quintal 2012; Kuteynikova et al. 2014).
Data Processing
The processing of the calculated synthetic shot records was carried out with ProMax software. In relation to the different trap models, the seismic data contain noise (such as scattering, edge reflections, multiple reflections) due to wave propagation. In general, a standard workflow was determined, including pre-processes, pre-stack and post-stack processes. (Fig. 4). However, considering the reflection quality of the data, a workflow update was performed, including appropriate parameter selection and process sequencing aimed at enhancing the primary reflections. In this context, the data processing processes applied to the seismic shot data calculated for each trap model are briefly summarised below.
In the granite wash modelling, bandpass filtering was applied to the artificial shot data at frequencies [10 16 60 80] to enhance the primary reflection, followed by predictive deconvolution with an prediction lag of 6 ms and an operator length of 80 ms to attenuate the multiple reflections. As a normal consequence of this process, since the noise amplitudes in the high frequency region outside the useful band of the data increase, bandpass filtering was applied to the data at cut-off frequencies of [8 10 70 90] Hz after deconvolution (Fig. 5a, c and Fig. 5b, d). Thus, with the application of deconvolution (Fig. 9c and 9d), the multiples were eliminated and the vertical resolution was increased. In Fig. 6 is given velocity spectrums plots which predictive deconvolution (right) applied and without predictive deconvolution applied (left). In the velocity spectrum calculation, since deep reflections become especially more prominent after deconvolution and their amplitude becomes stronger, the corresponding amplitudes of the velocity spectrum are also strengthened at that rate, and therefore the velocity estimation is also made more confidently.
When Fig. 7a, b and c are compared, it is observed that the reflections from the bottom of the fault indicated by the red arrows are distorted and the amplitudes are attenuated. It has been determined that this effect increases with reflector depth (Fig. 7b and c). Therefore, predictive deconvolution with an operator length of 80 ms and an prediction lag of 7 ms was applied to the shot data after the pre-data processing stages were applied. In this way, both the multiple reflections that are likely to occur are attenuated and the low amplitude primary reflection hyperbolas from each interface in the shadow zone become more expressible. The 330th CMP trace group is seen at approximately 1750 m of the normal fault trap model (Fig. 8). The resulting image is a velocity spectrum process based on the input-output energy ratio, and the most accurate velocity selection is made interactively on this spectrum.
In bioherm modelling, the coherence filtering proposed by many researchers (Neidell and Taner 1971; McMechan 1983; Leven and Roy-Chowdhury 1984; Kong et al. 1985; Milkereit and Spencer 1989; Li et al. 1997) was applied in the frequency-wavenumber domain for the filtering of edge reflections from the model edges. The coherent filter is used for direct waves, refracted waves and for the damping of all linear undesired events (noises). The velocity value and frequency range of the event to be filtered are required. Therefore, in the design of the filter, the velocity value was set to 2000 m/s and the frequency range was set to 20–60 Hz (Fig. 9b). As seen in Fig. 9, as a result of the application of the adaptation filter to the data, edge reflections were filtered and the signal-to-noise ratio of the shot group was increased.