Investigation of Pore Pressure on Seismic Amplitude Response : A Well-Based Modeling Case Study of Sojuko Field , Niger Delta Difference

The Sojuko field was discovered in 2001 in the eastern shallow offshore area of the Niger Delta, Nigeria. Three (3) exploration wells have so far been drilled in the field, two (2) of which are reasonably vertical and the third highly deviated. Three (3) key reservoirs which are laterally continuous across the wells have been identified with proven oil and gas reserves. Pore pressure data from repeat formation test (RFT) measurements acquired in the deviated well show that the wells are entirely hydrostatic to true depth (TD). This research focuses on investigating how seismic amplitudes change with offset/angle of incidence in relation to varying pore pressure regimes at the shale-hydrocarbon sand and shale-brine sand interfaces using well data. The aim is to aid quantitative interpretation in an on-going field-wide exploration drive to de-risk hydrocarbon exploration in the deeper plays in the area which are below TD, and are expected to be overpressured. The study is hinged on end-member shale elastic parameter substitution in which the shales are subjected to varying overpressure regimes while keeping the reservoirs (sands) at in situ (hydrostatic) condition. The end-member shale property substitution simulated shale compaction dis-equilibrium as the main overpressure generation mechanism in this study. The results show that top gas sands, top oil sands and top brine sands would be visible on seismic in the deeper plays where pore pressures are expected to be very high, but with distinctive seismic amplitude with offset/angle behavior. The top gas sands are visible as blue loop with small positive reflection coefficients at the near offsets/angles, but with polarity reversal to red loop with negative reflection coefficients which become more and more negative at the far angles at hard overpressure regimes. Top oil sands are recognized as blue loop with large positive reflection coefficients at the near angles; the coefficients becoming less and less positive at the far angles/offsets. The top oil sands may not be detected on seismic at the far angles/offsets unless at very hard overpressures. Brine sands have similar seismic response as oil sands at hard overpressures, but can be distinguished from oil sands based on their much higher amplitudes over the entire offset/angle range. The study is also aimed at removing uncertainty in seismic-based pore pressure quantification at the deeper targets where there is absence of well data for calibrating pore pressure effects at varying conditions.


Introduction
Seismic amplitude variation with shot-receiver offset or angle of incidence is an essential tool for oil and gas exploration.Quantitative interpretation of seismic amplitude response has been used as a good lithology indicator, and when properly calibrated, a good hydrocarbon indicator in a variety of basins around the world, including the Niger Delta.Seismic amplitude anomalies, however, would not provide direct hydrocarbon indicator (DHI) support at all times, as there have been instances in the Niger Delta, for example, when false positives, as well as hydrocarbon-bearing reservoirs without clear DHI support have been identified (Wojcik et al., 2016).The seismic response of a reservoir package is an interface property which depends on the elastic properties of the reservoir sandstones and the bounding shales across the shale-sandstone boundary.Well-derived rock property trends provide the basis for predicting seismic amplitude response for different combinations of lithology and pore fill.
Abnormally high pore pressures have been reported to be prevalent in the Niger Delta sedimentary basin, especially in the more prolific deeper plays which are currently the main exploration frontiers in the Niger Delta.Exploration and exploitation of hydrocarbon reserves in these overpressured zones could be quite expensive with a number of safety challenges.This does not preclude the fact that pore pressure may be higher at shallower depths in overpressured formations.Adequate knowledge of the pore pressure, especially pre-drill, is therefore very important to de-risk exploration and exploitation of these reserves and make the operations more cost effective.Direct, in situ measurements of pore pressure are made at discrete points in reservoir sections of a well.Shales are fine-grained and often impermeable.Pore fluids expelled from the pore spaces of the rocks arising from compaction due to geostatic load are retained within the impervious shale and as such, shales are often overpressured.Thus, overpressures in sand-shale sequences, especially in deltaic environments such as is the case in the current study, are mainly associated with under-compacted shales, although other factors such as tectonic activities, hydrothermal expansion and diagenesis may also play some roles.Characterizing overpressure in the shales would therefore quantify pore pressure effects associated with subsurface rocks in an exploration project.
Standard abnormal geopressure quantification methods using empirical relationships between seismic velocity and effective stress have been reported in the literature (Eberhart-Phillips et al., 1989), and have been extensively used for pre-drill pore pressure predictions in sedimentary basins.These empirical relationships are saddled with explicit and implicit assumptions; typical workflows determine elastic properties which are then used in simple shalebased poro-elastic formulations.The velocity-based pore pressure prediction methods are affected by factors such as lithology, fluid content and stress history (Ciz et al., 2005).
The amplitude of a P-wave reflected from an interface is dependent on the contrast in acoustic impedance across the interface.Considering a sand-shale lithologic sequence, for example, the reflection depends on the contrast in the acoustic impedance of a bounding shale lithology and the underlying reservoir sands.Acoustic impedance is a layer property which depends on the P-wave velocity and bulk density of rocks in the individual layers.In addition to the P-wave velocity and bulk density of the rocks, knowledge of the shear properties of the rock, particularly the shear velocity is also important in the determination of seismic amplitude variation with offset/angle of incidence.This arises from the physical condition that a P-wave incident obliquely at an interface is partitioned into a reflected and transmitted shear wave in addition to the reflected and transmitted P-waves which are also generated.
Seismic amplitude variation with offset has been used to study the expected behavior of hydrocarbon sands to aid quantitative interpretation of seismic data.In the Niger Delta, such studies have led to the discovery of giant oil and gas fields in the shallow marine/deltaic Nembe Creek Field, discovered in 1973 and the deepwater turbidites Bonga Field, discovered in 1996 (Nelson, 1973;Chapin et al., 2002;Wojcik et al., 2016).Changes in pore pressure can affect the properties of the rock frame, porosity and compaction which may be large enough to influence the seismic response (Lumley et al., 1987;Dolan, 2000).Prasad (2002) carried out laboratory measurement of amplitude spectra with increasing load, but did not draw explicit conclusions from the observations (Ciz et al., 2005).In this paper, we study the effect of pore pressure on seismic response in an oil and gas field in the Niger Delta shallow marine/deltaic system using well log and synthetic seismograms generated for different pore pressure regimes.The study aims to provide understanding of how seismic amplitudes change with offset/angle of incidence in relation to elevated pore pressures in the bid to understand the resultant changes in the rock's elastic properties.The study is important for removing uncertainty in pore pressure quantification at the deeper targets below well TD, where there is absence of well data for calibrating pore pressure effects at varying conditions.
The Zeoppritz equations (Equation 1) completely determine the amplitudes of seismic waves reflected from, or transmitted through a planar surface bounding two elastic media in welded contact with one another for all incidence angles (Haase, 2004).A number of linearized approximations have been made to the original Zeoppritz equations to simply calculations and offer more understanding of the factors which control seismic amplitude variation with offset/incidence angles, such as Bortfeld, 1961;Aki & Richards, 1980;Shuey, 1985;Smith & Gidlow, 1987;Fatti et al., 1994.In particular, the Bortfeld approximation which was used for this study, made it easier to understand how reflection amplitudes depend on incidence angle and physical parameters (Feng & Bancroft, 2006).In this approximation, the reflection coefficients are separated into three (terms) which include a normal incidence term, fluid factor term and a rigidity term on the basis of the Poisson's ratio (VerWest, 2004).The Bortfeld approximation is given by: where  is a ray parameter, defined as:  =   and  and  are acoustic impedance of the overlying and underlying layer, defined as   and  * , where  ,  and  ,  are the P-wave velocity and bulk density of the upper and lower layer, respectively. and  are the shear velocity of the upper and lower layer.

Location and Geological Setting
The study area is located in the eastern part of the Niger Delta shallow offshore, Nigeria (Figure 1).The Niger Delta occupies an area of about 75,000km 2 at the eastern part of the Gulf of Guinea, and is one of the world's most prolific hydrocarbon provinces.It is sub-divided into six megastructural units, also known as depobelts, which include the Northern Delta, Greater Ughelli, Central Swamp, Coastal Swamp, Shallow Offshore and the Deepwater, all of which have distinct sedimentary and structural styles, stratigraphy, as well as unique hydrocarbon generation, migration and distribution patterns.The structural styles, stratigraphy and petroleum system of the Niger Delta have been extensively studied and is well documented in the literature (Short & Stauble, 1967;Burke, 1972;Avbovbo, 1978;Evamy et al., 1978;Ejedawe, 1981;Whiteman, 1982;Doust & Omatsola, 1990;Reijers et al., 1997;Tuttle et al., 1999).The depositional system of the reservoir sands in the shallow offshore study area is of Late Miocene age, and is made up mainly of shallow sediments associated with the continued Niger Delta development.
Figure 1.Map of the Niger Delta showing the study area and the exploratory well locations.Wells XP2 and XP3 are vertical while XP1 is highly deviated

Materials
Of the three (3) exploratory wells in the field, only the deviated well has the complete suite of recorded well log data required for this study, including measured pore pressure (repeat formation tester, RFT) data.The suite of logs comprised GR, resistivity, density, neutron and compressional and shear sonic logs.The deviation data for the well was also provided.An arbitrary 3D pre-stack depth migration data traversing the three wells was provided for accurate wavelet extraction for the generation of synthetic seismograms.

Log Editing and Data Conditioning
The log data provided for this study were already edited by the petrophysicists in the team.It was noticed, however, that the editing was only restricted to the reservoirs and did not extend to non-reservoir layers.Minor mechanical editing was therefore carried out to edit obvious errors or spikes in the curves.The final log curves, especially the acoustic and density logs which are very important for this study were of good quality and did not require substitution by the most likely formation response interpreted from other logs.The GR log was used for sand-shale identification and volume of shale calculation and the resistivity logs helped in hydrocarbon identification for this study.The neutron log was used in combination with bulk density for gas-oil differentiation.The compressional and shear sonic slownesses were converted to compressional and shear velocities along hole which, when combined with bulk density, would provide seismic amplitude information for synthetic seismogram generation.
Figure 2 shows the conditioned well logs, converted to true vertical depth sub-sea.

Estimation of Optimum Seismic Wavelet
The seismic wavelet is a pulse of short duration (Sherrif, 2002); it is the link between the seismic data and subsurface rock properties, including stratigraphy (Cui & Margrave, 2014).Numerous seismic data analysis procedures such as seismic-to-well tie, seismic inversion and amplitude variation with offset/angle studies require knowledge of the embedded source wavelet which, basically, is the time-domain representation of the shape of reflection from a single positive reflector for a seismic wave incident normally at the reflection interface.Special techniques can be used in seismic data acquisition to measure the wavelet (Ikelle et al., 1997), but it is more commonly estimated from seismic data using statistical technique, as well as from the reflection coefficients or reflectivity series calculated from the sonic and bulk density logs.Whereas, well log data is able to provide information on both the amplitude and phase spectra of the wavelet, the statistical methods provide only amplitude spectrum and not the phase.The phase is assumed to be known in statistical wavelet estimation.Both the amplitude and phase of the wavelet are important for accurate interpretation of seismic data.Even though wavelet estimation using well logs can give the amplitude and phase information, a good tie between the well and seismic data must be achieved first to avoid errors in the resulting phase spectrum.Unfortunately, to obtain a good seismic-to-well tie, the optimum wavelet must have been known first.
In most practical eases of wavelet estimation using seismic and well data, a statistical wavelet is usually estimated from the seismic data as the first step and thereafter, the logs are matched to the seismic data for a new wavelet to be extracted using the well logs.This process may have to be repeated several times to obtain a reasonable wavelet from combination of the data.We used a different approach for wavelet extraction in this study.In the first instance, we created a 35Hz Ricker wavelet, sampled at 2ms for the Vp-Vs-RHOB log set for the well.A Ricker wavelet is often used as a zero-phase embedded wavelet in modeling and synthetic seismogram generation (Sherrif, 2002).
Figure 3 shows the amplitude and phase of the Ricker wavelet indicating positive polarity for acoustic impedance increase (blue loop) and a trough for acoustic impedance decrease (red loop).By comparing the statistical wavelet to the well-based wavelet, it was observed that the 35Hz well-based wavelet introduced so much of high frequencies into the data, although the amplitudes were reasonably comparable (Figure 5a).To derive an optimum wavelet, we used an approach which involves amplitude matching whereby keeping the sampling rate constant at 2ms, we varied the frequency in the Ricker wavelet estimation to match the amplitude spectrum of the statistical wavelet that was derived from well data (Figure 5b).In this approach, the frequency that best matches the amplitudes of the different wavelets is the optimum wavelet of the datasets.After several attempts, a 15.5Hz Ricker wavelet was obtained as the optimum wavelet and was used for subsequent modeling.

NMO Stretch and Full Waveform Synthetic
Normal moveout (NMO) is the travel time difference between a seismic wave at normal incidence and at a given offset.NMO correction is an important step in seismic data processing as it attempts to place all reflections at the zero offset.Dunkin & Levin (1973) were among the first researchers to discover that this correction causes "stretching" of the seismic waveform, resulting in increase in wavelength and loss of high frequencies, thereby distorting the amplitude spectrum of the wavelet.Wavelet stretching causes reduction in seismic resolution and adversely affects amplitude variation with offset/angle analysis (Zhang et al., 2011).The effect is more pronounced on the far offsets than the near data (Xu & Chopra, 2007).This problem is resolved at the data processing stage by the application of an offset mute, or an automatic mute designed on the basis of a pre-defined limit, known as stretch factor or tolerance , to allow only a certain range of useful offsets thereby removing the effect of NMO stretch from the data.Knowledge of the NMO stretch tolerance is very important in seismic response modeling to determine the maximum useful offset or angle that can be allowed in the data to prevent bias in the model.In this study, we carried out a full waveform AVO modeling on synthetic seismogram generated from the well data and previously derived 15.5Hz Ricker wavelet to determine the effect of NMO stretch in the data (Figure 6).Two (2) stretch factors,  = 1.2 and  = 1.3 were tested which defined the far angle at 35 0 and 45 0 , respectively.NMO stretch tolerance of 1.3 was chosen for the dataset.1.3 or 30% NMO stretch tolerance increases the wavelength of the data by 30% and decreases the frequency by 30%.The test shows that it is safe to model up to 45 0 far angle of incidence or its equivalent in offset without any adverse effect from NMO stretch in the area of the study (Figure 6).

Seismic Response Modeling
Setting the far angle to the maximum allowable stretch (45 0 ), we tested three (3) different ray-tracing-based AVO modeling methods to predict AVO response at the reservoir tops to determine the AVO modeling method most suitable for predicting seismic response at the reservoir tops in the study area.The methods tested were the exact Zeoppritz and it appromations such as Bortfeld and Bortfeld Aki-Richards Wiggins.The Bortfeld's AVO modeling method was chosen from this test and used in our subsequent modeling to study the effect of pore pressure on seismic response.
We hinged our seismic response modeling in relation to pore pressure on end-member shale substitution in which the shales were progressively overpressured while the reservoirs were kept at in situ pore pressure condition (hydrostatic).In the first step in the modeling, we modeled the lithostatic pressure  in the area using bulk density logs from the three (3) wells (Figure 7a).Next, we performed a regression analysis on the measured pore pressure (RFT) data which were made in the reservoirs, to obtain pore pressure  values in the non-reservoir sections (Figure 7b).The regression yielded the relation: () = 0.433 * () + 14.7 ( 6 ) Figure 7. Modeled lithostatic pressure (a) and measured pore pressure regression Thereafter, we calculated vertical effective stress in the shales using the Terzaghi's principle (Equation 7). where,
The vertical effective stress calculated in Equation 7is the vertical effective stress in the shales at in situ (hydrostatic) condition.We transformed this vertical effective stress in the shales to compressional velocity in the shales at hydrostatic condition, in accordance with the Bower's velocity-to-vertical effective stress model (Bowers, 1995).The velocity was derived using Equation 8: where  is compressional velocity at the sediment top which equals the water velocity, and  and  are Bower's parameters determined by calibrating compressional sonic velocity picks in shales with the measured pore pressure data.These parameters were respectively 1,524m/s, 0.3059 and 1.0124.
Having obtained compressional velocity in the shale at in situ condition, we carried out a shale property substitution whereby we raised the pore pressure in the shales to higher pore pressure regimes and modeled the resulting vertical effective stress for each case.We used two (2) overpressure scenarios, a mild overpressure regime with pore pressure gradient of 0.55psi/ft, and a hard overpressure regime with pore pressure gradient of 0.70psi/ft.For each regime, we computed the compressional shale velocity using Equation 3.This enabled us to have, at this point, compressional velocity in the shales for three pore pressure regimes: a hydrostatic (in situ) regime, mild overpressure regime and hard overpressure regime, with pore pressure gradients of 0.433psi/ft, 0.55psi/ft and 0.70psi/ft, respectively.
In addition to the compressional velocity, shear velocity and bulk density are the other elastic parameters which are important for seismic response modeling.We cross plotted the original compressional and shear velocity logs delimited by GR cut-off of 80 API (Figure 8) to obtain Vs versus Vp trend in shale (Equation 9), and then compressional and bulk density logs also delimited by the 80 API GR cut-off (Figure 9) to obtain shale bulk density versus compressional velocity trend (Equation 10).Using these trends, we derived the shear velocity and bulk density in the shales for the hydrostatic, mild and hard overpressure regimes by replacing Vp in Equations 9 and 10 with the output from Equation 8for the hydrostatic pressure, mild and hard overpressure scenarios respectively.Finally, we carried out the seismic response modeling with Bortfeld AVO modeling technique, having obtained the complete elastic parameters in the shales for the hydrostatic pressure, mild and hard overpressure cases respectively, to determine the effect of changing pore pressure on the seismic response.

Results and Discussion
The shale end-member substitution technique used in this study ensured that only the elastic properties of compressional and shear velocity, and bulk density in the non-reservoir sections of the well are subjected to different pore pressure regimes while the reservoirs are kept at in situ condition, which is hydrostatic.The modeling is focused on the shale-reservoir boundary only, and does not include fluid contacts in the reservoirs.Vertical effective stress computed from the lithostatic stress and measured pore pressure was the pivot of the modeling, since the elastic parameters were derived from regressions involving the compressional velocity computed from vertical effective stress for the respective pore pressure regimes.Figure 10 shows the result of elastic parameters derived for the hydrostatic pressure, mild and hard overpressure regimes.The result clearly shows that increase in pore fluid pressure results in decrease in the compressional and shear velocity, as well as bulk density of subsurface rocks, suggesting an inverse relation between the elastic parameters and pore fluid pressure.
Figure 11a shows the seismic response modeling for the hydrostatic pressure, mild overpressure and hard overpressure regimes for the shallowest reservoir, R100.This reservoir has a gross thickness of about 32 ft of gas sand, 41 ft of oil sand and 91 ft of brine sand.The respective reservoir units are indicated by the dashed horizontal lines in the first panel of Figure 11a.The reservoir exhibits a class III AVO response at the shale-gas sand boundary at in situ, hydrostatic pressure condition.The top gas sand is visible on seismic as a strong red loop with negative reflection coefficient; the reflection coefficients get more negative with angle of incidence/offset.The top reservoir is clearly visible on all angle ranges at this pressure regime without reversal in polarity.With increase in pore pressure, the reflection coefficients become less and less negative at the near angles (that is, the AVO intercept becomes less and less negative) and eventually, the gas sand changes to Class II and then Class I sands at hard overpressure, with evidence of polarity reversal at 32 0 incident angle from positive to negative amplitudes at the hard overpressure.The amplitude modeling shows that at hydrostatic pressure and mild overpressure regimes, top gas sand is visible as a red loop on seismic, but with very hard overpressures, top gas sand is recognized on seismic as a blue loop with positive reflection coefficient at the near angles, but the amplitudes reverse polarity to red loop with negative reflection coefficients at the far angles.Figure 11b shows the gas sand event for the different pore pressure regimes for R100, confirming Class III AVO response for the in situ case at pressure gradient of 0.433psi/ft, Class II response at mild overpressure of pressure gradient 0.55psi/ft and Class I at hard overpressure of pressure gradient of 0.70psi/ft.Similar results were obtained for the seismic response modeling at reservoir R200 which also contains aggregation of gas, oil and brine sands.With increase in pressure, polarity reversal occurred, and the reflection coefficients changed from negative to positive, even at mild over pressures.At hard overpressures, the top oil sand is visible on seismic as a blue loop with positive reflection coefficients, exhibiting a Class I seismic amplitude response with large positive amplitudes at the near angles which become less and less positive with increasing incidence angle.This amplitude response is clearly evident in the top oil sand event plot in Figure 12b which shows plots of the reflection coefficients versus angle for the respective pressure regimes.The seismic response modeling with respect to varying pore pressure regimes carried out in this study was done to examine the behavior of known hydrocarbon accumulations in the prospect area to provide better risk assessment for the deeper plays in the prospect area which are below the true depth of the well used for the modeling.In addition to providing understanding of the behavior of the hydrocarbon and brine bearing sands by understanding the resultant changes in elastic parameters due to changes in pore pressure, the study has provided additional information on the likely far angle ideal for amplitude modeling and inversion to aid quantitative interpretation in the area.Shale compaction disequilibrium is the major cause of overpressure encountered in most young, Tertiary sedimentary basins such as the Niger Delta, but other factors such as gas generation, diagenesis, tectonics and salt deposition can also generate overpressure.However, we have predicated our study on shale compaction disequilibrium as the main overpressure generation mechanism in the study area by the use of shale end-member elastic parameter substitution technique for this study.The modeling results have shown that changes in pore pressure in the prospect area can be detected on seismic.This result is consistent with those of Dunne (2013) who converted pore pressure measurements from Australia into logs of overpressure and cross-plotted with compressional velocity to study the effect of pressure on velocity and seismic response.

Conclusion
In the foregoing, we carried out end-member shale elastic parameter substitution to study the effect of changing pore pressure on seismic response at the shale-reservoir sand boundary for three proven hydrocarbon reservoirs and a brine reservoir, in an on-going effort to de-risk hydrocarbon exploration in the deeper plays in a shallow offshore prospect in the Niger Delta.The deeper plays are below the true depth of the existing wells, and are expected to be overpressured.The procedure involved subjecting the shales to progressively higher pore pressure regimes while keeping the reservoirs at in situ condition to study the change in seismic amplitudes resulting from changes in pore pressure.The shale property substitution technique used in this study simulated shale compaction disequilibrium as the main cause of overpressure in the area, this mechanism being the main cause of overpressure in young Tertiary sediments.The results show that pore pressure changes can be detected on seismic in the deeper plays in the field albeit the top gas sands, top oil sands and top brine sands have distinctive seismic amplitude behavior with offset/angle of incidence.The study also provides information on the maximum far angle/offset for future quantitative interpretation workflows and a robust seismic wavelet estimation method in addition to the expected behavior of gas, oil and brine sands at elevated pore pressure.

Figure 4 .
Figure 4. Statistical wavelet derivation across target reservoirs (a) and amplitude and phase spectra of the wavelet (b)

Figure 5 .
Figure 5. Derivation of optimum wavelet from seismic and well log data: (a) comparison of initial Ricker and statistical wavelets (b) amplitude matching of well and statistical wavelets

VsFigure 8 .
Figure 8. Vs versus Vp cross plot in shale: The Vs and Vp values for this cross plot were delimited using GR cutoff of 80 API to restrict the regression to the non-reservoir sections

Figure 10 .
Figure10.Elastic parameters resulting from shale-end member substitution for hydrostatic pressure (blue), mild overpressure (dark green) and hard overpressure (red) regimes.The elastic parameters decrease as a result of increase in pore pressure

Figure 11b .
Figure 11b.R100 seismic amplitudes versus angle of incidence plot for different pressure regimes.This reservoir contains gas, oil and brine.The reservoir top is the shale-gas sand boundary

Figure 12b .
Figure 12b.R300 seismic amplitudes versus angle of incidence plot for different pressure regimes.This reservoir contains oil and brine.The reservoir top is the shale-oil sand boundary