Pore pressure prediction in Niger Delta high pressure, high temperature (HP/HT) domains using well logs and 3D seismic data: a case study of X-field, onshore Niger Delta

Few wells targeting high temperature, high pressure intervals in most tertiary sedimentary basins have achieved their objective in terms of technicalities and cost. Since most shallow targets have been drilled, exploration focus is drifting into deeper plays both onshore and in deep offshore areas. To ensure safe and economic drilling campaigns, pore pressure prediction methodologies used in the region needs to be improved. The research aims at generating and testing a modification of Eaton’s equation fit for high temperature, high pressure intervals on a field. The evolution of pore pressure in the field was established from offset well data by making several crossplots, and fracture gradient was computed using Mathew and Kelly’s equation. Eaton’s equation parameters were then calibrated using several wells until a desired field scale result was achieved when compared with information from already drilled intervals i.e., kicks and RFT data. Seismic velocity data resulting from high density, high resolution velocity analysis done to target deep overpressured intervals were then used to predict 1D pore pressure models at six selected prospect locations. Analyses reveal depths shallower than 3800 m TVD/MSL with geothermal gradient 3.0 °C/100 m and pressure gradient less than 1.50sg EMW are affected mainly by undercompaction; depths greater than 3800 m TVD/MSL with geothermal gradient of 4.1 °C/10 m and pressure gradients reaching 1.82–2.12sg EMW are affected by unloading with a narrow drilling margin for the deep highly pressured prospect intervals. Eaton’s n-exponent was modified to 6, and it proved accurate in predicting high overpressure in the first prospect wells drilled.


Introduction
Rapid sedimentation associated with young tertiary sedimentary basins [e.g., Niger Delta, Gulf of Mexico, South East Asia, Trinidad] have led to entrapment of pore fluids in thick sequences of low permeability shales [Holland et. al. 1990;Nehring, 1991;Grauls and Baleix 1994;Heppard et. al. 1998]. The reduction in porosity and permeability associated with this burial coupled with increase in geothermal gradient with depth sometimes lead to the development of deep high pressured, high temperature intervals. These intervals create challenges for both exploration and drilling activities as high overpressures might affect the hydrocarbon retention capacity of traps due to hydrofracturing of seals as well as pose additional risk and uncertainties to drilling teams and facilities.
If overpressures are not accurately predicted prior to drilling, it can greatly increase drilling non-productive time, result in reduced drilling speed and may cause catastrophic incidents such as well blowouts, pressure kicks and well complexities (Zhang 2011;Zhang 2013;Sen and Ganguli 2019;Baouche et al. 2020aBaouche et al. , b, 2020cGanguli and Sen 2020;Radwan and Sen 2020). Other minor incidents that might occur include unintended fracturing of the formation and consequent mud losses. While it is rare to find a well that did not meet at least a part of its objectives, only few wells drilled in Niger Delta to target high pressure, high 1 3 temperature intervals achieved their stated objectives. Some drilling campaigns have either failed, been terminated or resulted in actual drilling cost very high compared to well cost estimates (Hope 2017).
In most explored intervals of the Niger Delta basin, the dominant primary overpressure mechanism is undercompaction due to high sedimentation rate which is more rapid than fluid expansion rate. Based on the latest downhole measurements, shallow sandstone reservoirs are found to have sub-hydrostatic pore pressure, while the lower reservoir unit is presently at hydrostatic pressure regime (Agbasi et al. 2020). As drilling activities are now focused into deeper plays, secondary overpressure mechanisms become increasingly important which necessitates a proper understanding of these mechanisms and modified approach to prediction techniques.
This study integrates well logs, seismic velocity data, 3D seismic data, reservoir formation pressures (RFT data), leak-off test (LOT data), temperature data from several wells in the field to generate a modified Eaton's equation which would efficiently predict pore pressure evolution Onshore, Niger Delta. The objective of the project includes (a) to identify the different overpressure mechanisms in the Niger Delta (b) to identify regions impacted by secondary overpressure mechanisms (c) to set regional shale normal compaction trends and compute overburden gradients (d) to predict pore pressure in disequilibrium compactions and secondary overpressure domains (e) to compute fracture gradient (FIP) profiles, (f) to build cube models of pore pressure gradient, overburden gradient and fracture gradient for the field.

Theoretical background
Overpressure generating mechanisms has been broadly classified into five categories; namely disequilibrium compaction, fluid expansion, diagenesis, tectonic compression and pressure transfer (Zhao et. al 2018). The cause of overpressure in many sedimentary basins is undercompaction i.e., disequilibrium compaction which makes it the primary overpressure mechanism. Secondary overpressure mechanisms include fluid expansion which is related to thermal stress (hydrocarbon generation/kerogen cracking, oil cracking to gas, hydrothermal expansion), diagenesis (illite-smectite transformation, quartz dissolution, cementation) and lateral and vertical overpressure transfer (dynamic transfer) (Osborne and Swarbrick 1997). An example of fluid transfer is that of a sealed interval having pore fluid pumped into it from another, higher-pressure zone. Sometimes this can be caused by charging along faults or along dipping sand enclosed in shale. The sand transmits pore fluid and pressure from deeper shales up dip (Hall 1993;Swarbrick et al. 2002;Yardley and Swarbrick 2000). Causes of overpressure differ in lithology; for mudstones, the overpressured formation in source rocks is usually different from those of non-source rocks, the former which is frequently related to hydrocarbon generation and sometimes also affected by diagenesis, while the latter which is commonly related to disequilibrium compaction, diagenesis and pressure transfer. Phenomena such as quartz dissolution and cementation is also important in some Cenozoic basins especially in Nigeria where fine grained, poorly sorted texturally mature, platykurtic siltstones are heavily affected by cementation and compaction (Alabere et al. 2020). Dehydration reactions associated with mineral diagenesis is another mechanism of overpressure development. Smectite dehydration is a complex process but can lead to overall volume increases of both the rock matrix and the pore water system (Hall 1993). Fluid expulsion from montmorillonite lattices during phase transition to illite typically occur at a temperature of 100 °C in the Gulf of Mexico and this correlate to the depth at which overpressure has been observed (Bruce 1984). The transition of anhydrite to gypsum is another dehydration reaction that can lead to overpressure development, but only at relatively shallow depths as the temperature at which this dehydration occurs is only about half that of the smectite-illite transition.
Several direct (empirical analysis) and indirect (theoretical analysis) have been proposed for the study of overpressure. The empirical analysis refers to the analysis of the geologic-geophysical logging response to and the experimental test of overpressures, while the latter refers to the analysis of geologic conditions for and numerical simulation of overpressured formation (Zhao et al. 2018). With more and more application of empirical methods in the study of overpressured formations, most of the overpressure cases that are traditionally thought to be caused by disequilibrium compaction are denied totally or partially. Instead hydrocarbon generation is demonstrated to be of increasing significance, the clay diagenesis (especially illite-smectite transformation) as well as tectonic compression and pressure transfer also started becoming important in recent literatures. In addition, the overpressured formation in many basins is thought to be influenced by the combination of two or more overpressuring mechanism (Zhao et al. 2018).
Six methods have been proposed for the study of overpressure, and they include; Bower's method, velocity-density crossplotting method, porosity correlation method, log curves combination analysis, pressure calculation and correlation method and the comprehensive analysis method (Zhao et al. 2018). In this research, a combination of velocity-density crossplotting method, pressure calculation and correlation method was used.
The Bowers approach proposed by Bowers (1994) which involves the use of porosity-vertical effective stress chart which is also known as the Bower's chart or loading-unloading curves pointed out that overpressures caused by disequilibrium compaction would reflect on the loading curve while those caused by fluid expansion would reflect on the unloading curve. He proposed two key indicators for the identification of unloading origin; which are that the velocity of sediments should be very small when compared with the noticeable change of effective stress and that density is not sensitive to unloading (Bowers 1994(Bowers , 2002. Therefore, to determine overpressures caused by unloading it's necessary to make both the effective stress-velocity relation diagram and the effective stress-density relation diagram when compiling loading-unloading charts. Sediment unloading can be elastic, inelastic or transitional. For elastic unloading, the acoustic velocity decreases while the density remains constant in the overpressured interval (Bowers 1994), or porosity slightly increases due to elastic rebound while the velocity decreases significantly (Bowers 2011). In the case of inelastic unloading, the density increases while the acoustic velocity slightly decreases or remains constant (Lahann and Swarbrick 2011).
The proponents of the velocity-density crossplotting method (Bowers 1995(Bowers , 2011Lahann et al. 2001;Hoesni 2004;O'Conner et al. 2011;Tingay et al. 2013;Lahann and Swarbrick 2011) have showed the relationship between different overpressuring mechanism and variation of density and velocity. Fluid expansion mechanisms such as gas generation results in a decrease in acoustic velocity as overpressure increases with little or no effect on the density (Hoesni 2004). Clay diagenesis or chemical compaction results in increasing density as overpressure increases with little or no effect on the acoustic velocity (Lahan et al. 2001;Lahann and Swarbrick 2011;O'Conner et al. 2011). On the contrary load transfer or combination of causal mechanisms usually leads to decreasing acoustic velocity, increasing density as overpressure increases but the magnitudes of the changes are between those described above (Zhao et al. 2018). Tingray et al. (2013) suggested that overpressure caused by normal hydrostatic pressure and disequilibrium compaction are located on the loading curve while those caused by other mechanism are on the unloading curve. For the use of acoustic velocity and density crossplots, care as to be taken to ensure correction for organic matter content as they have obvious effects on the density and acoustic velocity of formations (Li et al. 2016). Figure 1 shows the evolution of pore pressure with depth and the effective stress response. A couple of researchers (Hertmanrud et al. 1988;Tiege et al. 1999;Tingray et al. 2009;Dasgupta et al. 2016) used the porosity correlation method which depends on comparison between porosity of overpressured mudstones with normally pressured intervals at similar depth to determine the presence of a porosity anomaly which would indicate the activity of an overpressure generating mechanism. The log curves combination analysis which is the most basic method and one of the most reliable in investigating overpressure origin in use since 1970s (Fertl and Timko 1972;Fertl 1976). Some researchers such as Bowers (2002) believe that when multiple logs are used for overpressure prediction, a synchronous reversal of acoustic, resistivity and density logging curves should occur. On the contrary if the reversal is not synchronous i.e., if density reversal lags behind, that indicates that the overpressure present is not caused by undercompaction but by other mechanisms such fluid expansion, pressure transfer, tectonic loading or smectiteillite transformation.
The pressure calculation and correlation method involves calculating overpressure using different methods (Eaton's method, equivalent depth method, the Bower's method, etc.) and correlating it with measured formation pressures. The mechanism by which the calculated formation pressures are most approximate to the measured pressure is most likely the cause of the overpressure. The equivalent depth method is based on porosity and thus effectively predicts disequilibrium compaction while the Eaton's equation prediction depends on the value of the Eaton's exponent. An exponent (n) of 3.0 suggests disequilibrium compaction while exponent (n) greater than 6 could be predictive of fluid expansion or pressure transfer mechanisms. Whichever of these methods that matches the measured pressure in the wells ready suggests the possible overpressure mechanism at play. The comprehensive analysis combines the empirical methods described above with knowledge of geological setting and distribution of overpressures as well as simulation studies (Zhao et al. 2018). Knowledge of the basin would involve understanding the geothermal and compaction setting and depositional environment.

Materials and methods
A dataset that comprised 3D seismic, reservoir formation pressure (RFT), leak-off test (LOT), Suite of well logs (wireline and LWD logs consisting of sonic, density, gamma ray and resistivity logs), mudlogging data and end of well reports were employed in this study. Our analysis was carried out in five stages: (1) Offset wells post-mortem pore pressure analyses and establishing the cause of overpressure.
To investigate the presence and understand control on distribution of high temperature and pressure intervals in the field, temperature measurements from 17 wells across the field were plotted with depth, and the geothermal gradient was calculated using Eq. 1. Also, plots of reservoir pressure (RFT) and leak-off pressure (LOT data) versus depth were made to see the evolution of reservoir pressure and fracture gradient with depth, and the regional minimum stress was estimated from the leak-off pressure. The cause of the observed high temperature and high overpressure in the field was investigated using acoustic velocity vs density crossplots and comparison with well-known standard (Fig. 2) modified by Swarbrick (2012). This is then followed by high density, high resolution velocity analysis. The velocity picking process is done using Delta Stack, an automatic 3D velocity picking tool. The 3D HDHR velocity fields were obtained by carrying out a Normal Moveout Correction (NMO) correction in a time domain at every CMP gather. We then define a propagation constraint on a seed line in 3D which drove the picking Fig. 2 Schematic diagram of velocity-density showing the normal compaction trend and the possible trends for the different overpressure generating mechanism (Swarbrick 2012after Hoesni 2004 and is automatically propagated and adapted according to neighboring values already available. The interval velocity cube is then computed from this RMS velocity field at two offset (post-mortem) well locations and Dix converted (Eq. 2), then pseudo-sonic was computed using Eq. 3. Anisotropy corrections were applied to obtain a pseudo-sonic consistent with actual sonic acquired in the wells. This was done to calibrate the pseudo-sonic and determine what corrections needed to be applied at any prospect well location. 1D interval velocities were picked at six random locations to cover the area where the possible deep prospects are located and extended to reach deep prospect intervals between 4200 and 4500 m TVD/MSL.
The LWD and Wireline logs had gamma ray log, density and sonic log though the density and GR log were missing in some wells. Where there was no density log, a pseudo density log was created using Gardner's equation (Eq. 4) and applying the constant A & B calibrated from wells where both sonic and density logs were present.
The gamma ray logs were first interpreted to delineate sandstones from shaly rocks, then a line group was drawn using a cut-off of 75 API which was increased in some areas to 80 API for some wells to ensure that only the shales were included in the line groups used to create shale points. The shale points were then transferred to the sonic log so as to have interval transit time for only the shaly intervals as reservoir pressure would be an input from the RFT data. Then logs were then filtered using shrink boxcar (equal weight method) to reduce the density of points. In wells were there was no gamma ray logs, the whole sonic logs were filtered as shale points cannot be picked.
Since the density logs are usually acquired in the reservoir sections of interest, thus lacking at shallow intervals, there is a need to extrapolate the density to the surface. This is important for overburden gradient computation as it involves integrating the density log from the surface till total depth (TD) of the well as all the overlying weight of rocks has to be considered. The extrapolation was done using Miller's equation (Eq. 5) and then all the depth-corrected density logs for the 11 wells as well as the Rhob_miller (density from miller's equation) was averaged to get the regional Rhob_average log which is filtered using the shrink boxcar (4) = A * 10 6 DT sonic B (equal weight) approach and then displayed in all wells after the depth was corrected for ground level elevation due to differences in topography. The correction was necessary because the reference for overburden gradient calculations onshore is the ground level. The overburden gradient was then computed by integrating the depth-corrected regional Rhob_average in each well from the top to TD (Eq. 6).
The normal compaction trend was then computed for all eleven wells using NCT-3P equation (Eq. 7), and calibrated using values from the hydrostatic shale points on sonic log by adjusting the DT max , DT min and K parameters of the NCT-3P until the resulting normal compaction trend best fits the filtered sonic log. Any sharp points of departure that cannot be fitted already gives an indication of overpressure.
The shale pore pressure was computed using Eaton's sonic equation (Eq. 8) and the computed normal compaction trend, overburden gradient, shale point filtered sonic log and Eaton exponent, n = 3 for regions shallower than 3800 m TVD/MSL as the input parameters. For regions deeper than 3800 m TVD/MSL which are impacted by unloading, the normal Eaton's equation could not adequately account for the pressures observed in shales assuming pressure equilibrium between shale and reservoir sands in these depth intervals. Pressure estimates came from kicks recorded in thin isolated sand layer within the shales. Modified Eaton's equation was used for pore pressure prediction below this depth by fitting shale pore pressure predicted to the observed events in these wells (kicks and RFT) by adjustment of the Eaton exponent. The Eaton's n-parameter was tested in the range of 3-6, with n = 6 best predicting the shale pore pressure in the deep pressured intervals. Measured reservoir pressure from RFT measurements was integrated into the shale pore pressure prediction to get the final pore pressure profile.
The fracture gradient was then computed using Matthews & Kelley's equations (Eq. 9), and the effective stress ratio (K) was calibrated by using the leak-off test results as control to see the minimum and maximum case.

3D cube modeling
The 1D pore pressure models created using well data and seismic velocities give a good representation of the pressure regime at the well scale which is sparsely located across the field. As these wells are a few kilometers apart, the entire process of pore pressure prediction using seismic velocities would have to be repeated for each future exploration well location. The purpose of 3D modeling is to create cubes of the geopressure data that would give a more cohesive evolution of pore pressure components across the field such that 1D pressure profiles could be extracted from the cube at any prospect or infill well location. Drillworks 3D application was used to create 3D models of geopressure data from 1D post-drilling and pre-drilling models. Drillworks 3D software converts the 1D interpretation of geopressure data from all wells into 3D models taking into account the spatial distance between these one-point datasets and interpolating using kriging.

Results and discussion
The results from well log analysis, seismic interpretation and crossplots revealed the origin and causes of overpressure. Normally, geothermal gradient in tertiary sedimentary basin like the Niger Delta is around 3 °C per 100 m of descent (9) F g = K * (S − PP) + PP through the crust but temperature versus depth plot shown in Fig. 3 reveals two geothermal gradient trends. A trend in the shallow depths (< 3800 m TVD/MSL) with the normal geothermal gradient and a second geothermal gradient of 4.1 °C/100 m in the deeper intervals (> 3800 m TVD/ MSL) suggestive of temperatures higher than 110 °C beyond this depth. Crossplots of the reservoir pressure and leak-off pressure (LOT) as shown in Fig. 4 revealed that pressures shallower than 3700 m TVD/MSL are mostly hydrostatic or lower due to ongoing production from the field while virgin reservoir pressures are observed below this depth. Also, much higher reservoir pressure trends are observed below 4000 m TVD/MSL depth. This proves that high temperature, high pressure zones exist at depth deeper than 3800 m TVD/ MSL. Based on the pressure and LOT points on this plot, a synthetic regional minimum principal stress (S3) was estimated (Fig. 5) that gives an idea of the expected hydrocarbon column retention limit which if exceeded, the hydrocarbon trap would be considered blown or leaked. Most overpressure mechanisms result in an unexpected decrease in density and acoustic velocity with increase in depth. Thus, a cross-plot of these two properties reveals that the overpressure mechanism in the shallower intervals (< 3800 m TVD/MSL) is undercompaction/disequilibrium compaction. For deeper intervals (> 3800 m TVD/MSL) where higher temperatures have been recorded, overpressure appears to be due to a combination of unloading mechanisms (fluid expansion and fluid transfer) and clay dehydration (Fig. 6). From the surface, the pore pressure is hydrostatic, to the top of overpressure where pressure increases gradually to a maximum of about 1.5sg EMW within the intervals < 3800 m TVD MSL (Fig. 9). The top of overpressure is shallower in the northern part and deepest in the southeastern part, a direct effect of complex fault system on the field. Onset of overpressure is observed to be deeper for the exploration wells, which suggests that pseudo-sonic extraction from interval velocity tends to underestimate the onset of overpressure in the field. A general drop in reservoir pressure was observed between 2000 and 3600mTVD/ MSL interval due to depletion by ongoing production from these reservoirs. This zone with depleted reservoirs might be susceptible to differential sticking risks during drilling operations. The top of hard overpressure coincides with the onset of secondary overpressure mechanism which causes the sharp pressure ramp observed between 380 and 4000 m TVD/MSL interval with pressure gradient reaching 1.80 to 2.12sg EMW. This corresponds to zones in which unloading (fluid transfer/expansion and clay dehydration) is at play in addition to undercompaction.
The range of fracture gradient values is between 1.40 and 1.6sg EMW at shallow depths and at deeper intervals that there is a sharp increase to between 1.80 and 2.20 sg EMW. This suggest that drilling with mud weight heavier than 2.2sg would leak to fracking of the formation while drilling below 1.4sg might result in kicks. The earliest onset of observed reduction in drilling window for the field is at 3700 m TVD/MSL as most well would already encounter a narrow window (lesser than 0.54). Figures 10 and 11 below show that at depths shallower than 3700 m TVD/MSL, the drilling window is between 0.54 and 0.7sg EMW. Below 3900 m TVD/MSL which correspond to overpressured  interval, the drilling window starts to become narrow and drops to a range of between 0.4 and 0.17sg EMW at around 4300 m TVD/MSL. Accurate understanding and observance of mud weight within the drilling window are required to prevent kicks and wellbore instability issues.

Uncertainty and sensitivity analysis
Uncertainty in shale pore pressure prediction comes from uncertainty in the parameters used in its estimation such as the normal compaction trend, overburden gradient and the Eaton's exponent. Uncertainty in the overburden gradient calculation is in the range of (± 0.13-0.21) of the actual value and in the Eaton's exponent is between 3 and 6. Normal compaction trend changes based on the choice of DTmax, DTmin and compaction factor. All this put together meant that the final shale pore pressure is within 0.9sg of the actual shale pore pressure in equivalent mud weight as seen in Fig. 12.

Conclusion
(1) There is two main overpressure regimes onshore Niger Delta. The first regime is caused primarily by disequilibrium compaction and occurs from the top of undercompaction (ranging from 1000 to 2200 m TVD/MSL) to about 3800 m TVD/MSL where pressure gradient is less than 1.50sg EMW.
(2) Below 3800 m, TVD/MSL is the second regime caused by secondary overpressure mechanism (unloading) in addition to disequilibrium compaction. This regime is marked by a sharp pressure ramp with pressure gradients reaching 1.82-2.2sg EMW below 3900 m TVD/ MSL. (3) Kerogen transformation, thermal cracking of oil and clay dehydration reactions which all occur at temperatures above 100 °C are probable causes of fluid volume increase which coupled with expected low permeability at such dept resulted in the observed secondary overpressure. (4) Deeper than 3700 m, TVD/MSL depth narrow drilling margin due to high overpressures associated with unloading pose serious operational challenges if conventional drilling methods are used. (5) The 3D cubes of the most likely pore pressure and fracture gradient prediction will be a cost saving tool as the companies can plan their base case architecture of future infill and exploration wells on the most likely case and potentially save money on unnecessary casings. Also, the cost of exploration wells will reduce if such wells are not planned on the max case (commitment case) prediction but on the base case (most likely case). This will increase ability to drill more wells within the exploration budget. However, there is still a potential risk which necessitates follow up on pore pressure while drilling so as to quickly react to any unpredicted events.