Application of pre-stack seismic waveform inversion and empirical relationships for the estimation of geomechanical properties in Ruby field, central swamp depobelt, Onshore Niger Delta, Nigeria

Pre-stack seismic inversion, well log analysis approach and empirical relations were adopted in this study to better estimate geomechanical properties of Ruby field with minimum error. The use of conventional well log empirical method alone to evaluate geomechanical properties in oil/gas fields sometimes becomes problematic. Geomechanical properties were divided into: elastic moduli [Young’s modulus, shear modulus, bulk modulus and Poisson ratio (PR)] and rock mechanical strength properties (closure stress ratio (CSR), brittleness (BRI) and compressibility). Four geomechanical earth models (CSR, BRI, Young’s modulus and PR) were generated from the inversion analysis to understand the distribution of rock strength properties across the field. The results deciphered high Young’s, shear and bulk modulus in the reservoir zone compared to the cap/seal rocks and a decrease in PR. This implies that, the cap/seal are more ductile and less compressible than the reservoir rocks, indicating that the reservoirs are highly brittle. CSR result reveals high in cap/seal indicating that the cap/seal rock are harder to fracture and has a greater chance to withstand higher compressive stress before failing as opposed to reservoir rocks. The inverted earth model shows that, Young’s modulus and brittleness increase toward the northeastern part of the field, while CSR and PR increase toward the southwestern part of the field. These results suggest that harder, stiffer, highly compressible and easily fractured rocks are found in the northern and eastern part of the field as opposed to the southern to western part of the field that is ductile.


Introduction
Over the years, the evaluation of geomechanical properties in the Niger Delta oil and gas field has been neglected. This is because the field rarely suffers from most geomechanicalrelated problems such as reservoir compaction, subsidence and wellbore stability problems. However, the estimation of elastic moduli and rock mechanical properties plays a salient role in reservoir analysis from developmental stage down to production, completion and abandonment stage. Understanding the strength of reservoirs, cap/seal rocks helps us to better predict the competence of the formation and hence aid in planning drilling location and wellbore stability programme. However, two main hurdles are in borehole stability studies often due to: (a) data regarding rock mechanical properties, insitu stress state and geologic structures that cannot be accurately define and are often provided with uncertainties, (b) well engineers, petrophysicist and drillers who are directly involved in well delivery may not have adequate knowledge in identifying and managing borehole stability without looking at the geomechanical properties of the reservoir (Xiaohui et al. 2016). Therefore, the estimation, evaluation and establishment of geomechanical earth model enhance understanding of the field in solving most geomechanical issues. The estimation of elastic moduli (Young's Modulus) and rock mechanical strength properties (Closure Stress ratio, Brittleness and Poisson ratio) from seismic data has over the years posed great challenge to geophysicist and petroleum engineers in the Niger Delta Basin. Advancement in seismic technology has led to the ease in estimation of geomechanical properties. Advances on the use of Prestack time-migrated (PreSTM) data have enormously helped in evaluating geomechanical properties with minimum error, thereby reducing the problem posed by inadequate knowledge of these properties (Okeugo et al. 2018). The issue of band-limit vertical resolution and uniqueness of the deterministic inversion results has led to the use of Prestack seismic inversion for the evaluation of geomechanical properties (Trudeng et al. 2014). One fundamental advantage of using Prestack simultaneous waveform inversion is the ability of the inversion process to remove the effect of wavelet from the seismic data (removal of side lobe and tuning effect introduced by wavelet) and produce high-resolution image of the subsurface (Margrave et al. 1998;Veeken and Silva 2004;Moosavi and Mokhtari 2016). Clearly, another advantage of the seismic-driven geomechanical model is that Prestack seismic inversion provides additional means to constrain the uncertainties of mechanical properties in the 3D space and better account for local variation introduced by complex geological structures (Xiaohui et al. 2016). As a necessary and integral part of field development, a full-field geomechanical model is essential to understand reservoir and the overburden behavior. Estimation of geomechanical properties in the Niger Delta Basin using prestack timemigrated seismic data and well log requires great skill and immense understanding of seismic inversion technique as this is the key to unraveling geomechanical properties from seismic data. Seismic inversion method provides the most detailed view of the earth's subsurface and its applications in the Niger Delta over the years have yielded increase in oil production (Dagogo et al. 2016;Akpan et al. 2020a, b).
Some studies have been done to evaluate geomechanical/rock mechanical properties in the Niger Delta Basin (Eyila and Oladunjoye 2014; Abijah and Tse 2016;Osaki et al. 2018;Abe et al. 2018), but none employed pre-stack seismic inversion technique to evaluate these properties. Most of the study focussed on the use of well log data to evaluate these properties and developed a one-dimensional (1D) and 3D seismic stochastic inversion which is basically geostatistical and are filled with uncertainties in the development of mechanical earth models. Abijah and Tse (2016) evaluated geomechanical properties of an onshore oil and gas field in the Niger Delta. In their findings, significant variation of elastic moduli was observed between cap and reservoir rocks. The cap rock, which is shale, has high average Poisson's ratio, elastic, bulk and rigidity modulus of 0.4,23,466.1 MPa,253,312.8 MPa and 8,717 MPa, respectively. The reservoir rocks, which were inferred to be sandstones revealed relatively lower Poisson's ratio, elastic, bulk and rigidity modulus, but higher compressibility and rock strength making them more brittle with higher potential of tensile fracture. Eyila and Oladunjoye (2014) estimated geomechanical strength of reservoir rocks with well log data for safety limit in sand free production. The study considered bulk modulus (K) and shear-to-compressibility ratio (S/c)) for solving sanding problem in the reservoir. It was observed that the average value of K is lower than the threshold value, which indicates the minimum value at which fluid may be produced safely at any rate, and S/c is lower but tends toward acceptable range. Results from geomechanics help to address production-related challenges such as seal integrity, reservoir competence, fault reactivation, compaction, subsidence, injection, depletion, borehole stability and sand control.
In view of these, several researches were carried out within the study area. This study focuses on the estimation of geomechanical properties (elastic moduli and rock mechanical) from well data using empirical equations and seismic data using prestack time-migrated seismic waveform inversion in Ruby field reservoir competent assessment.

Geology of the study area
The Niger Delta Basin is situated at the top of the Gulf of Guinea on the west coast of Africa (Fig. 1a). The Basin lies between latitudes 3°N and 6°N and longitudes 4°E and 8°E (Ejedawe et al. 2007). The base map of the study area showing locations of wells in Ruby field is shown in Fig. 1b. Niger Delta Basin is one of the most prolific deltaic hydrocarbon provinces in the world. The reservoirs of the field have been dated to Middle Miocene and identified to be deposited in a para-sequence of shallow marine and deltaic plain deposits (Reijers 2011). The Delta has one petroleum system called the Akata-Agbada system of the Tertiary Niger Delta (Chukwuemeka et al. 2018). It is made up of three lithostratigraphic units, ranging from the oldest Akata Formation, Agbada Formation and the youngest Benin Formation (Fig. 2), which are all diachronous. The Akata Formation is the oldest lithostratigraphic unit in the Niger Delta. It is Palaeocene to Recent in age and occurs as thick marine sedimentary succession that is laid in front of an advancing delta. It is estimated to be about 21,000 ft thick in the central part of the clastic wedge. It consists of thick shale sequence and turbidite sands in deep waters with little amount of clay and silt. The shale serves as the main source rock in the basin while the turbidites serve as reservoirs in the deep offshore Niger Delta (Doust and Omatsola 1990).
The Agbada Formation (paralic clastic) is the second of the three strongly diachronous Niger Delta multifaceted Formations (Short and Stauble 1967;Frankl and Cordry 1967). It is characterized by paralic interbedded sandstone and shale (Reijers et al. 1996). In the lower Agbada Formation, sand and shale were deposited in equal proportion while the upper portion is mostly sand with minor shale interbedded. The continental Benin Formation comprises the top part of the Niger Delta clastic wedge, from the Benin-Onitsha area in the north to beyond the coastline (Short and Stauble 1967). The top of the formation is a recent, sub-aerially exposed delta top surface and its base extends to a depth of 4600 ft. The base of the Benin Formation is defined by the youngest marine shale. The shallow parts of this Formation are composed entirely of non-marine sands deposited in alluvial or upper coastal plain environments during a southward progradation of the delta (Doust and Omatsola 1990).

The conceptual framework
Pre stack seismic simultaneous inversion is a stratigraphic deconvolution technique used for reservoir characterization due to its high resolution and ability to track thin beds (Okeugo et al. 2018). One fundamental advantage of using this method of seismic inversion technique for constructing an earth model of various geomechanical properties is its ability to extract additional information such as shear velocity (V S ) and density (ρ) from seismic data. This gives the method an edge over other techniques where these properties (shear velocity and density) are determined using empirical formulae. The Zoeppritz first-order approximation to reflectivity is given by Aki and Richard (1980) in Eq. 1 as: It was observed that for the equation to be correct, the relative changes of the property ( ΔV P∕ V P , ΔVs ∕ V P and Δ ∕ ) should be small so that the second-order terms can be neglected and that is much less than 90 0 . With these assumptions, Eq. 1 can be rewritten in terms of P-wave and S-wave impedances as: where R(θ) is the P-wave reflection coefficient, θ is the angle of incidence, Z P = V P is the average acoustic impedance, Z S = V S i s t h e ave r a g e s h e a r i m p e d a n c e , the zero-offset S-wave reflection coefficient. The P-wave reflection coefficient was later simplified by Fatti et al. (1994) as given in Eq. 3 as: But Ma (2002) replaced V S V P by Z S Z S , such that the reflection coefficients R( ) are only a function of three parameters: Z P , Z S and . With this, Eq. 3 then reduces to: Equation 4 serves as the bedrock of prestack seismic simultaneous inversion employed in estimating acoustic impedance and shear impedance (Z P and Z S ) for prestack seismic gather. Once acoustic impedance (Z P ), shear impedance (Z S ) and density (ρ) volume have been generated, geomechanical properties such as Young's modulus, brittleness, closure stress ratio and Poisson's ratio can also be generated (Hampson et al. 2005).
(2) Mavko et al. (2009) defined dynamic Young's modulus as the extensional stress to the extensional strain ratio in a uniaxial stress state. Assuming there is no change in the orthogonal normal stresses, the ratio of incremental normal stress to its resultant strain gives the Young's modulus for any mechanical body (Bjørlykke et al. 2015). Dynamic Young's modulus can be obtained from elastic wave velocity from well logs (Archer and Rosouli 2012;Eyinla and Oladunji 2014;Jamshidian et al. 2017) as given in Eq. 5: where ρ is density, V P is the compressional wave velocity and V S is the shear wave velocity.
Shear or rigidity modulus is the degree to which a material can be subjected to a shearing stress and still resist shearing. Dynamic shear modulus can be estimated from density measurements and sonic derived P-wave using the relation of Mavko et al. (2009) given in Eq. 6 as: Bulk modulus (Kb) defines the extent or measure of a substance resistance to isotropic squeezing. A material subjected to an isotropic pressure from P 1 to pressure P 2 , will compress from its initial volume (V 1 ) to a new smaller volume V 2 . Dynamic bulk modulus can be computed from velocity and density data using relationship proposed by (Jizba 1992;Brinson and Brinson 2008;Archer and Rasouli 2012) given in Eq. 7 as: Poisson's ratio (PR) is the measure of longitudinal contraction due to lateral elongation caused by the application of a uniaxial stress (Fjar et al. 2008). Poisson's ratio is estimated using the relation given Eq. 8 as: Closure stress ratio (CSR) is the degree of pressure at which fracture closes after the fracturing pressure is relaxed. Rocks with high closure stress are basically harder to fracture than rocks with lower closure stress. This parameter is computed empirically using Hampson et al. (2005) relation given in Eq. 9 as: where the parameters are as earlier defined.
Brittleness (BRI) defines the ability of a material to break with little elastic deformation and significant plastic deformation when subjected to some degree of stress. Brittle rocks are harder and rigid but with very little tensile strength (Hampson et al. 2005). This geomechanical property can be estimated empirically using Eq. 10: where E index is the Young modulus index given as: and PR index is the Poisson ratio index given as: E is the Young's modulus, E min and E max are the minimum and maximum values of Young's modulus, respectively. Also, PR is the Poisson's ratio, PR min and PR max are the minimum and maximum values of Poisson's ratio, respectively.

Well control
Well log method of geophysical exploration depends on the ability of the interpreter to make geologic inference away from the well bore. The extrapolation of well data outside the area of interest is often necessary as it aid in the identification of faults, sequence boundaries and pressure compartments. The presence of well control in any field plays a vital role in seismic inversion method of geophysical exploration. As stated by Simm and Becon (2014), the problem of determining wavelet phase becomes much more difficult when there is no well control. Wavelet phase controls the quality of the inversion output. The presence of more well control in any field will aid in validating the output of the inversion results based on the inverted amplitudes around the well. With this, deductions can be made on the inverted amplitudes away from the well. In general, having more well control will improve the significance of the inverted seismic results especially away from the well which can lead to drilling more production wells in the field.

Materials and methodology
The available dataset employed in this study include 3D prestack time-migrated (PreSTM) seismic volume and suites of wire line log labelled Well_001, Well_025 and Well_029. The log comprises gamma ray, density-neutron, resistivity and compressional wave velocity as shown in Fig. 3. The PreSTM seismic volume and wells were quality-checked (QC) and prepared for seismic simultaneous inversion. This is because the quality of seismic volume is critical in building a reliable earth models. Seismic data are often populated by artefact which has strong impact on processing and interpretation. Geostatistical data enhancement was used in order to identify artefact or anomalies in seismic data. This   process also helped in coherence of event related to seismic velocities. However, check shot correction was also carried out to properly scale the compressional velocity (Vp) log. P-wave velocity was determined from sonic log using the Ogagure (2008) localized relationship model for Niger Delta Basin given in Eq. 13 as: where Δt P is the interval transit time for compressional sonic log (Table 1). Shear wave velocity (V S ) was determined from the compressional wave velocity using Greenberg and Castagna (1992) relation given in Eq. 14 as: For a successful evaluation of geomechanical properties from seismic data, wavelet was extracted (Fig. 3) and tie to wells to properly calibrate the time event on both seismic and well log. The extracted wavelet was observed to be symmetrical with its maximum amplitude value at zero time. This symmetrical shape aided the picking of reflection peak and trough for the generation of earth models of various geomechanical properties. Synthetic seismogram was also generated in the wells to aid in the inversion. It was observed that the synthetically generated information from the well data (blue traces) was not properly matched with the original seismic data (red traces). Time shift and stretching were applied to the synthetic traces to fit the original seismic data. After this operation was done, a good to excellent correlation coefficient of 92.7% and 94.1% were obtained for Well_001 and Well_029 wells, respectively, as shown in Figs. 4 and 5.
Four elastic moduli (Young's modulus, Shear modulus, Bulk modulus and Poisson's ratio) were estimated (13) V P = 0.305 Δt P × 10 6 (14) V P = 1.16V S + 1.36  dynamically using V P from depth of 11,047-11,936 ft in Well_001 and 10,643-11,803 ft in Well_029. Statistical values were estimated and plotted against depth to determine the trend for both wells (Well_001 and Well_029). These parameters were estimated using the relations given in Eq. 5-12 accordingly. Earth models were estimated for Young modulus, closure stress ratio, brittleness and Poisson's ratio from the inversion analysis after successful seismic well to tie. The methodology workflow employed in achieving the aim of this study is shown in Fig. 6.

Elastic moduli
The results of four elastic moduli: Young's modulus, Shear modulus, Bulk modulus and Poisson's ratio estimated dynamically using V P for Well_001 and Well_029 are shown in Tables 2 and 3. The tables show statistical  Hakiki et al. (2018). This means that the cap/seal are more ductile and less compressible than the reservoir rocks, indicating that the reservoir are stiffer, harder and more compressible and hence competent in hydrocarbon exploitation.

Rock mechanical properties
Rock strength properties were estimated to evaluate the competence of the formation between the hydrocarbon and non-hydrocarbon unit. Three rock strength properties were estimated: closure stress ratio (CSR), brittleness (BRI) and compressibility (λ). Among the estimated parameters, CSR and BRI were seen to play a great role in proper formation strength evaluation. Tables 4 and 5 show values for estimated properties used for the evaluation. These properties were estimated from V P and density logs and were plotted against depth (Figs. 8 a  and b). It was observed that, the reservoirs are highly brittle ranging from 9.59 to 20.24% than the cap/seal rocks meaning the aforementioned is more ductile than the reservoir rocks. If both formations are subjected to high stress, the reservoir will fracture first before the cap/seal rocks, thus suffers less elastic deformation. The reservoir is inferred to be prone to high tensile fracture as compared to the cap/seal rocks. CSR was observed to be low in the reservoirs ranging from 3.6 to 4.6 Gpae −1 and high in cap/seal indicating that the cap/seal rocks are harder to fracture than the reservoir rocks. The observed fracture will tend to close rapidly in the cap/seal as compared to the reservoir rocks when the fracturing pressure is relaxed.

Prestack seismic simultaneous inversion estimated Geomechanical properties
The evaluation of geomechanical properties from seismic data requires a good estimation of amplitude and phase spectra. The inversion analysis results for the two wells are presented in Figs. 9 (a and b). The inverted trace for Z P , Z S , density and V P -V S ratio are displayed in track 1, track 2, track 3 and track 4, respectively. The inverted results (in red) overlay the original logs from well (in blue) and initial model in (black). Analysis of the inverted parameters shows some slight mismatch in the initial, original and inverted model of Z P , Z S and density, respectively, at some intervals. Although the inverted result was not affected, the error estimated between the density, synthetic and seismic trace remains approximately zero. However, the error gather, which is the difference between the two results, shows an approximately zero error of 0.091. The small error observed indicates that the inversion is mathematically correct and subsequent synthetic trace which matches the real seismic trace was generated. The observed zero error implies that the inversion algorithm has created impedance trace that tends to be consistent with the extracted wavelet and input seismic trace.

Young's modulus (E)
Result from the Prestack seismic inverse modeling shows high to low values of Young's modulus values ranging from 13.28 to 19.46 GPa (Fig. 10a). The horizon time slice generated from sand body-2 divides the horizon into three zones: high, intermediate and low Young's modulus zones. These zones were identified by the values of Young's modulus within the horizon with high zone above 17.66 GPa, intermediate zone ranging from 15.08 to 17.66 GPa and low zone below 15.08 GPa. The 3D volume generated using a time slice of 2875 ms identified harder and stiffer rocks across the study area (Fig. 10b). From the zones, the degree of rock strength and competence of the horizon slice and across the study area were established. Since, elastic moduli and rock strength are linearly proportional; it indicates that an increase in elastic moduli results in an increase in rock strength. This implies that, high Young's modulus zones signify harder, stiffer and incompressible rock with good degree of competence. It is also inferred from the 3D volume result (Fig. 10b), that harder, stiffer and high degree of competent rocks are found at northeastern part of the field.

Closure stress ratio (CSR)
The result of inverse modeling of the CSR shows high to low values of CSR ranging from 3.7 × 10 −4 to 6.8 × 10 −4 (Fig. 11a). A horizon time slice generated from sand body-2 reveals high, intermediate and low CSR zones.
These zones were marked out through the values of CSR as deciphered in the colour key with high zone ranging above The 3D volume which was generated to evaluate the response of the rock body to fracture when stress is applied across the study area is shown in Fig. 15b. From the result, three fracture response zones of the rock body were determined. These zones are grouped as: highly prone area to fracture and less fracture movement to closure (low CSR), intermediately prone area to fracture (intermediate CSR) and less prone area to fracture with high fracture movement to closure (high CSR). However, CSR was observed to be greater at the south to western part of the field, making rocks in this zone harder to fracture than rocks in the northeastern part of the field, as increase in CSR depicts hardness to fracture. When fracturing pressure is relaxed in the field (if any), the fracture will tend to close rapidly in the southwestern part of the field as compared to the northeastern part of the field (Fig. 11b).

Brittleness (BRI)
The brittleness result ranges from 16 to 40% across the field as deciphered from the colour key. The horizon time slice generated from sand body-2 divides the horizon into three zones: high, intermediate and low brittleness zones as shown in Fig. 12a. These zones were identified by the values of BRI within the horizon with high zone above 32.5%, intermediate zone ranging from 23.4 to 32.5% and low zone below 23.4%. The 3D volume generated using horizon time slice of 2875 ms to evaluate the level of brittleness (easiness to break) across the study area reveals three levels of brittleness (Fig. 12b)  to western part of the field. Conversely, the southwestern part is seen to be more ductile than the northeastern part. This implies that if fracturing stress is applied, the northeastern part of the field will fracture first before the southwestern part but will suffer more elastic deformation than the northern part of the field. However, the northeastern part of the field is prone to tensile failure and the southwestern part is prone to compressive shear failure as shown in Fig. 12b.

Poisson ratio (PR)
The inverse modeling result shows high value to low values of PR ranging from 0.25 to 0.36 as deciphered in the colour of Fig. 13a. From the horizon time slice generated using sand body-2, three PR zones: high, intermediate and low are observed. These zones were identified by the values of PR within the horizon with high zone above 3.3 × 10 −4 , intermediate zone ranging from 3.0 × 10 −4 to 3.3 × 10 −4 and low zone below 3.0 × 10 −4 GPa. However, a 3D volume generated with a time slice of 2875 ms identified harder and stiffer rocks across the study area (Fig. 13b). From the zones, the degree of rock strength and competence of the horizon slice across the study area was established. It was observed from the volume that harder, stiffer and high degree of competent rocks is found at the northeastern part of the field.

Conclusion
An integrated prestack simultaneous seismic inversion and empirical relation determination approach from both seismic and well logs data were employed to elucidate the geomechanical characteristics of the ''Ruby'' field, within the Central Swamp depobelt of the Niger Delta for reservoir competence assessment. These datasets provided a robust interplay that enabled the delineation of elastic moduli and rock mechanical strength properties within the well locations from empirical relations across the field and from prestack seismic inversion analysis. Evaluation of geomechanical properties within the zone of interest (reservoir and cap/seal) reveals increase in elastic moduli, brittleness and decrease in closure stress ratio and Poisson's ratio. The reservoirs are generally inferred to be brittle while the cap/seal are ductile in nature. Geomechanical properties evaluated from the inversion analysis of Young's modulus volume reveal that harder rocks are found at the northern part of the field. Closure stress ratio decreases toward the north and increase in the south indicating that rocks in the south are harder to fracture. Brittleness increases toward the north indicating that the rocks in the north are brittle and prone to easy break if exposed to stress, but will suffer less elastic deformation as opposed to the rocks in the southwestern part of the field, which will suffer great elastic deformation when fractured and hence prone to compressive shear failure as compared to rocks in the north that may be prone to tensile failure. Utilizing 3D prestack seismic inversion constrain by local geology has yielded a high-resolution predictive model for reservoir competence assessment and improve significantly wellbore stability prediction. This will help to reduce development well cost by having a more precise and predictable geomechanical model.